diff --git a/docs/paper/reductions.typ b/docs/paper/reductions.typ index 15c431e52..376712b9e 100644 --- a/docs/paper/reductions.typ +++ b/docs/paper/reductions.typ @@ -154,6 +154,7 @@ "ConsecutiveBlockMinimization": [Consecutive Block Minimization], "ConsecutiveOnesMatrixAugmentation": [Consecutive Ones Matrix Augmentation], "ConsecutiveOnesSubmatrix": [Consecutive Ones Submatrix], + "FeasibleBasisExtension": [Feasible Basis Extension], "SparseMatrixCompression": [Sparse Matrix Compression], "DirectedTwoCommodityIntegralFlow": [Directed Two-Commodity Integral Flow], "IntegralFlowHomologousArcs": [Integral Flow with Homologous Arcs], @@ -6607,6 +6608,101 @@ A classical NP-complete problem from Garey and Johnson @garey1979[Ch.~3, p.~76], ] } +#{ + let x = load-model-example("FeasibleBasisExtension") + let A = x.instance.matrix + let m = A.len() + let n = A.at(0).len() + let rhs = x.instance.rhs + let S = x.instance.required_columns + let cfg = x.optimal_config + // Free column indices (those not in S) + let free-cols = range(n).filter(j => j not in S) + // Selected free columns from config + let selected = cfg.enumerate().filter(((i, v)) => v == 1).map(((i, v)) => free-cols.at(i)) + // Full basis: required + selected + let basis = S + selected + [ + #problem-def("FeasibleBasisExtension")[ + Given an $m times n$ integer matrix $A$ with $m < n$, a column vector $overline(a) in bb(Z)^m$, and a subset $S$ of column indices with $|S| < m$, determine whether there exists a _feasible basis_ $B$ --- a set of $m$ column indices including $S$ --- such that the $m times m$ submatrix $A_B$ is nonsingular and $A_B^(-1) overline(a) >= 0$ (componentwise). + ][ + The Feasible Basis Extension problem arises in linear programming theory and the study of simplex method pivoting rules. It was shown NP-complete by Murty @Murty1972 via a reduction from Hamiltonian Circuit, establishing that determining whether a partial basis can be extended to a feasible one is computationally intractable in general. The problem is closely related to the question of whether a given linear program has a feasible basic solution containing specified variables. The best known exact algorithm is brute-force enumeration of all $binom(n - |S|, m - |S|)$ candidate extensions, testing each for nonsingularity and non-negativity of the solution in $O(m^3)$ time.#footnote[No algorithm improving on brute-force enumeration is known for the general Feasible Basis Extension problem.] + + *Example.* Consider the $#m times #n$ matrix $A = mat(#A.map(row => row.map(v => str(v)).join(", ")).join("; "))$ with $overline(a) = (#rhs.map(str).join(", "))^top$ and required columns $S = \{#S.map(str).join(", ")\}$. We need $#(m - S.len())$ additional column from the free set $\{#free-cols.map(str).join(", ")\}$. Selecting column #selected.at(0) gives basis $B = \{#basis.map(str).join(", ")\}$, which yields $A_B^(-1) overline(a) = (4, 5, 3)^top >= 0$. Column 4 makes $A_B$ singular, and column 5 produces a negative component. + + #pred-commands( + "pred create --example " + problem-spec(x) + " -o feasible-basis-extension.json", + "pred solve feasible-basis-extension.json --solver brute-force", + "pred evaluate feasible-basis-extension.json --config " + x.optimal_config.map(str).join(","), + ) + + #figure( + canvas(length: 0.7cm, { + import draw: * + let cell-size = 0.9 + let gap = 0.15 + // Draw the matrix + for i in range(m) { + for j in range(n) { + let val = A.at(i).at(j) + let is-basis = j in basis + let is-required = j in S + let f = if is-required { + graph-colors.at(1).transparentize(50%) + } else if is-basis { + graph-colors.at(0).transparentize(30%) + } else { + white + } + rect( + (j * cell-size, -i * cell-size), + (j * cell-size + cell-size - gap, -i * cell-size - cell-size + gap), + fill: f, + stroke: 0.3pt + luma(180), + ) + content( + (j * cell-size + (cell-size - gap) / 2, -i * cell-size - (cell-size - gap) / 2), + text(8pt, str(val)), + ) + } + } + // Column labels + for j in range(n) { + let label-color = if j in S { graph-colors.at(1) } else if j in basis { graph-colors.at(0) } else { black } + content( + (j * cell-size + (cell-size - gap) / 2, 0.4), + text(7pt, fill: label-color)[$c_#j$], + ) + } + // Row labels + for i in range(m) { + content( + (-0.5, -i * cell-size - (cell-size - gap) / 2), + text(7pt)[$r_#(i + 1)$], + ) + } + // RHS vector + let rhs-x = n * cell-size + 0.6 + content((rhs-x + (cell-size - gap) / 2, 0.4), text(7pt, weight: "bold")[$overline(a)$]) + for i in range(m) { + rect( + (rhs-x, -i * cell-size), + (rhs-x + cell-size - gap, -i * cell-size - cell-size + gap), + fill: luma(240), + stroke: 0.3pt + luma(180), + ) + content( + (rhs-x + (cell-size - gap) / 2, -i * cell-size - (cell-size - gap) / 2), + text(8pt, str(rhs.at(i))), + ) + } + }), + caption: [Feasible Basis Extension instance ($#m times #n$). Orange columns are required ($S = \{#S.map(str).join(", ")\}$), blue column is the selected extension. Together they form a nonsingular basis with non-negative solution.], + ) + ] + ] +} + // Completeness check: warn about problem types in JSON but missing from paper #{ let json-models = { diff --git a/docs/paper/references.bib b/docs/paper/references.bib index 3b63d53a5..f9f3064f0 100644 --- a/docs/paper/references.bib +++ b/docs/paper/references.bib @@ -1454,3 +1454,13 @@ @techreport{plaisted1976 number = {STAN-CS-76-583}, year = {1976} } + +@article{Murty1972, + author = {Murty, Katta G.}, + title = {A fundamental problem in linear inequalities with applications to the travelling salesman problem}, + journal = {Mathematical Programming}, + year = {1972}, + volume = {3}, + pages = {326--370}, + doi = {10.1007/BF01584550} +} diff --git a/problemreductions-cli/src/cli.rs b/problemreductions-cli/src/cli.rs index 9c40cfdd9..b71d1b244 100644 --- a/problemreductions-cli/src/cli.rs +++ b/problemreductions-cli/src/cli.rs @@ -272,6 +272,7 @@ Flags by problem type: ConsecutiveOnesMatrixAugmentation --matrix (0/1), --bound ConsecutiveOnesSubmatrix --matrix (0/1), --k SparseMatrixCompression --matrix (0/1), --bound + FeasibleBasisExtension --matrix (JSON 2D i64), --rhs, --required-columns SteinerTree --graph, --edge-weights, --terminals MultipleCopyFileAllocation --graph, --usage, --storage AcyclicPartition --arcs [--weights] [--arc-costs] --weight-bound --cost-bound [--num-vertices] @@ -734,6 +735,12 @@ pub struct CreateArgs { /// Query attribute index for PrimeAttributeName #[arg(long)] pub query: Option, + /// Right-hand side vector for FeasibleBasisExtension (comma-separated, e.g., "7,5,3") + #[arg(long)] + pub rhs: Option, + /// Required column indices for FeasibleBasisExtension (comma-separated, e.g., "0,1") + #[arg(long)] + pub required_columns: Option, /// Number of groups for SumOfSquaresPartition #[arg(long)] pub num_groups: Option, diff --git a/problemreductions-cli/src/commands/create.rs b/problemreductions-cli/src/commands/create.rs index 0fdb3fe8e..a3e266d3a 100644 --- a/problemreductions-cli/src/commands/create.rs +++ b/problemreductions-cli/src/commands/create.rs @@ -9,7 +9,7 @@ use anyhow::{bail, Context, Result}; use problemreductions::export::{ModelExample, ProblemRef, ProblemSide, RuleExample}; use problemreductions::models::algebraic::{ ClosestVectorProblem, ConsecutiveBlockMinimization, ConsecutiveOnesMatrixAugmentation, - ConsecutiveOnesSubmatrix, SparseMatrixCompression, BMF, + ConsecutiveOnesSubmatrix, FeasibleBasisExtension, SparseMatrixCompression, BMF, }; use problemreductions::models::formula::Quantifier; use problemreductions::models::graph::{ @@ -193,6 +193,8 @@ fn all_data_flags_empty(args: &CreateArgs) -> bool { && args.conjuncts_spec.is_none() && args.deps.is_none() && args.query.is_none() + && args.rhs.is_none() + && args.required_columns.is_none() } fn emit_problem_output(output: &ProblemJsonOutput, out: &OutputConfig) -> Result<()> { @@ -766,6 +768,9 @@ fn example_for(canonical: &str, graph_type: Option<&str>) -> &'static str { "StringToStringCorrection" => { "--source-string \"0,1,2,3,1,0\" --target-string \"0,1,3,2,1\" --bound 2" } + "FeasibleBasisExtension" => { + "--matrix '[[1,0,1,2,-1,0],[0,1,0,1,1,2],[0,0,1,1,0,1]]' --rhs '7,5,3' --required-columns '0,1'" + } _ => "", } } @@ -893,6 +898,9 @@ fn help_flag_hint( } ("ConsecutiveOnesSubmatrix", "matrix") => "semicolon-separated 0/1 rows: \"1,0;0,1\"", ("SparseMatrixCompression", "matrix") => "semicolon-separated 0/1 rows: \"1,0;0,1\"", + ("FeasibleBasisExtension", "matrix") => "JSON 2D integer array: '[[1,0,1],[0,1,0]]'", + ("FeasibleBasisExtension", "rhs") => "comma-separated integers: \"7,5,3\"", + ("FeasibleBasisExtension", "required_columns") => "comma-separated column indices: \"0,1\"", ("TimetableDesign", "craftsman_avail") | ("TimetableDesign", "task_avail") => { "semicolon-separated 0/1 rows: \"1,1,0;0,1,1\"" } @@ -2831,6 +2839,53 @@ pub fn create(args: &CreateArgs, out: &OutputConfig) -> Result<()> { ) } + // FeasibleBasisExtension + "FeasibleBasisExtension" => { + let usage = "Usage: pred create FeasibleBasisExtension --matrix '[[1,0,1],[0,1,0]]' --rhs '7,5' --required-columns '0'"; + let matrix_str = args.matrix.as_deref().ok_or_else(|| { + anyhow::anyhow!( + "FeasibleBasisExtension requires --matrix (JSON 2D i64 array), --rhs, and --required-columns\n\n{usage}" + ) + })?; + let matrix: Vec> = serde_json::from_str(matrix_str).map_err(|err| { + anyhow::anyhow!( + "FeasibleBasisExtension requires --matrix as a JSON 2D integer array (e.g., '[[1,0,1],[0,1,0]]')\n\n{usage}\n\nFailed to parse --matrix: {err}" + ) + })?; + let rhs_str = args.rhs.as_deref().ok_or_else(|| { + anyhow::anyhow!( + "FeasibleBasisExtension requires --rhs (comma-separated integers)\n\n{usage}" + ) + })?; + let rhs: Vec = rhs_str + .split(',') + .map(|s| s.trim().parse::()) + .collect::, _>>() + .map_err(|err| { + anyhow::anyhow!( + "Failed to parse --rhs as comma-separated integers: {err}\n\n{usage}" + ) + })?; + let required_str = args.required_columns.as_deref().unwrap_or(""); + let required_columns: Vec = if required_str.is_empty() { + vec![] + } else { + required_str + .split(',') + .map(|s| s.trim().parse::()) + .collect::, _>>() + .map_err(|err| { + anyhow::anyhow!( + "Failed to parse --required-columns as comma-separated indices: {err}\n\n{usage}" + ) + })? + }; + ( + ser(FeasibleBasisExtension::new(matrix, rhs, required_columns))?, + resolved_variant.clone(), + ) + } + // LongestCommonSubsequence "LongestCommonSubsequence" => { let usage = @@ -7612,6 +7667,8 @@ mod tests { storage: None, quantifiers: None, homologous_pairs: None, + rhs: None, + required_columns: None, } } diff --git a/src/models/algebraic/feasible_basis_extension.rs b/src/models/algebraic/feasible_basis_extension.rs new file mode 100644 index 000000000..d80cc1ac9 --- /dev/null +++ b/src/models/algebraic/feasible_basis_extension.rs @@ -0,0 +1,349 @@ +//! Feasible Basis Extension problem implementation. +//! +//! Given an m x n integer matrix A (m < n), a column vector a_bar of length m, +//! and a subset S of column indices with |S| < m, determine whether there exists +//! a feasible basis B (a set of m column indices including S) such that the +//! m x m submatrix A_B is nonsingular and A_B^{-1} a_bar >= 0. +//! +//! NP-complete (Murty, 1972). + +use crate::registry::{FieldInfo, ProblemSchemaEntry}; +use crate::traits::Problem; +use serde::{Deserialize, Serialize}; + +inventory::submit! { + ProblemSchemaEntry { + name: "FeasibleBasisExtension", + display_name: "Feasible Basis Extension", + aliases: &[], + dimensions: &[], + module_path: module_path!(), + description: "Given matrix A, vector a_bar, and required columns S, find a feasible basis extending S", + fields: &[ + FieldInfo { name: "matrix", type_name: "Vec>", description: "m x n integer matrix A (row-major)" }, + FieldInfo { name: "rhs", type_name: "Vec", description: "Column vector a_bar of length m" }, + FieldInfo { name: "required_columns", type_name: "Vec", description: "Subset S of column indices that must be in the basis" }, + ], + } +} + +/// The Feasible Basis Extension problem. +/// +/// Given an m x n integer matrix A with m < n, a column vector a_bar of length m, +/// and a subset S of column indices with |S| < m, determine whether there exists +/// a feasible basis B of m columns (including all of S) such that the submatrix +/// A_B is nonsingular and A_B^{-1} a_bar >= 0. +/// +/// # Representation +/// +/// Each non-required column has a binary variable: `x_j = 1` if column j is +/// selected. A valid config must select exactly m - |S| additional columns. +/// The problem is satisfiable iff the resulting A_B is nonsingular and +/// A_B^{-1} a_bar >= 0. +/// +/// # Example +/// +/// ``` +/// use problemreductions::models::algebraic::FeasibleBasisExtension; +/// use problemreductions::{Problem, Solver, BruteForce}; +/// +/// let matrix = vec![ +/// vec![1, 0, 1, 2, -1, 0], +/// vec![0, 1, 0, 1, 1, 2], +/// vec![0, 0, 1, 1, 0, 1], +/// ]; +/// let rhs = vec![7, 5, 3]; +/// let required = vec![0, 1]; +/// let problem = FeasibleBasisExtension::new(matrix, rhs, required); +/// let solver = BruteForce::new(); +/// let solution = solver.find_witness(&problem); +/// assert!(solution.is_some()); +/// ``` +#[derive(Debug, Clone, Serialize, Deserialize)] +pub struct FeasibleBasisExtension { + matrix: Vec>, + rhs: Vec, + required_columns: Vec, +} + +impl FeasibleBasisExtension { + /// Create a new FeasibleBasisExtension instance. + /// + /// # Panics + /// + /// Panics if: + /// - The matrix is empty or has inconsistent row lengths + /// - m >= n (must have more columns than rows) + /// - rhs length does not equal m + /// - |S| >= m (must have room for at least one additional column) + /// - Any required column index is out of bounds + /// - Required columns contain duplicates + pub fn new(matrix: Vec>, rhs: Vec, required_columns: Vec) -> Self { + let m = matrix.len(); + assert!(m > 0, "Matrix must have at least one row"); + let n = matrix[0].len(); + for row in &matrix { + assert_eq!(row.len(), n, "All rows must have the same length"); + } + assert!( + m < n, + "Number of rows ({m}) must be less than number of columns ({n})" + ); + assert_eq!( + rhs.len(), + m, + "rhs length ({}) must equal number of rows ({m})", + rhs.len() + ); + assert!( + required_columns.len() < m, + "|S| ({}) must be less than m ({m})", + required_columns.len() + ); + for &col in &required_columns { + assert!(col < n, "Required column index {col} out of bounds (n={n})"); + } + // Check for duplicates + let mut sorted = required_columns.clone(); + sorted.sort_unstable(); + for i in 1..sorted.len() { + assert_ne!( + sorted[i - 1], + sorted[i], + "Duplicate required column index {}", + sorted[i] + ); + } + Self { + matrix, + rhs, + required_columns, + } + } + + /// Returns the matrix A. + pub fn matrix(&self) -> &[Vec] { + &self.matrix + } + + /// Returns the right-hand side vector a_bar. + pub fn rhs(&self) -> &[i64] { + &self.rhs + } + + /// Returns the required column indices S. + pub fn required_columns(&self) -> &[usize] { + &self.required_columns + } + + /// Returns the number of rows (m). + pub fn num_rows(&self) -> usize { + self.matrix.len() + } + + /// Returns the number of columns (n). + pub fn num_columns(&self) -> usize { + self.matrix[0].len() + } + + /// Returns the number of required columns (|S|). + pub fn num_required(&self) -> usize { + self.required_columns.len() + } + + /// Returns the indices of non-required columns (the "free" columns). + fn free_columns(&self) -> Vec { + let required_set: std::collections::HashSet = + self.required_columns.iter().copied().collect(); + (0..self.num_columns()) + .filter(|c| !required_set.contains(c)) + .collect() + } + + /// Check if basis columns form a nonsingular system and the solution is non-negative. + /// + /// Uses exact rational arithmetic via integer Gaussian elimination with + /// numerator/denominator tracking to avoid floating-point errors. + #[allow(clippy::needless_range_loop)] + fn check_feasible_basis(&self, basis_cols: &[usize]) -> bool { + let m = self.num_rows(); + assert_eq!(basis_cols.len(), m); + + // Build augmented matrix [A_B | a_bar] in i128 to avoid overflow + // during Bareiss fraction-free Gaussian elimination. + let mut aug128: Vec> = Vec::with_capacity(m); + for i in 0..m { + let mut row = Vec::with_capacity(m + 1); + for &col in basis_cols { + row.push(self.matrix[i][col] as i128); + } + row.push(self.rhs[i] as i128); + aug128.push(row); + } + + // Bareiss algorithm: fraction-free Gaussian elimination. + // After elimination, the system is upper-triangular. + let mut prev_pivot = 1i128; + + for k in 0..m { + // Partial pivoting + let mut max_row = k; + let mut max_val = aug128[k][k].abs(); + for i in (k + 1)..m { + if aug128[i][k].abs() > max_val { + max_val = aug128[i][k].abs(); + max_row = i; + } + } + if max_val == 0 { + return false; // singular + } + if max_row != k { + aug128.swap(k, max_row); + } + + for i in (k + 1)..m { + for j in (k + 1)..=m { + aug128[i][j] = + (aug128[k][k] * aug128[i][j] - aug128[i][k] * aug128[k][j]) / prev_pivot; + } + aug128[i][k] = 0; + } + prev_pivot = aug128[k][k]; + } + + // Back-substitution to solve. We solve in rational form: x_i = num_i / det. + // The solution x = A_B^{-1} a_bar must satisfy x >= 0, which means + // num_i / det >= 0 for all i, i.e., num_i and det have the same sign (or num_i = 0). + // Back-substitution using rational arithmetic to check x >= 0. + // Simple rational back-substitution: + // x[i] = (aug128[i][m] - sum_{j>i} aug128[i][j] * x[j]) / aug128[i][i] + // We track x[i] as (numerator, denominator) pairs. + + let mut x_nums = vec![0i128; m]; + let mut x_dens = vec![1i128; m]; + + for i in (0..m).rev() { + // numerator of (aug128[i][m] - sum_{j>i} aug128[i][j] * x[j]) + let mut num = aug128[i][m]; + let mut den = 1i128; + + for j in (i + 1)..m { + // subtract aug128[i][j] * (x_nums[j] / x_dens[j]) + // num/den - aug128[i][j] * x_nums[j] / x_dens[j] + // = (num * x_dens[j] - den * aug128[i][j] * x_nums[j]) / (den * x_dens[j]) + let a = aug128[i][j]; + num = num * x_dens[j] - den * a * x_nums[j]; + den *= x_dens[j]; + // Simplify to avoid overflow + let g = gcd_i128(num.abs(), den.abs()); + if g > 1 { + num /= g; + den /= g; + } + } + // x[i] = (num/den) / aug128[i][i] = num / (den * aug128[i][i]) + let diag = aug128[i][i]; + x_nums[i] = num; + x_dens[i] = den * diag; + // Normalize sign: make denominator positive + if x_dens[i] < 0 { + x_nums[i] = -x_nums[i]; + x_dens[i] = -x_dens[i]; + } + let g = gcd_i128(x_nums[i].abs(), x_dens[i].abs()); + if g > 1 { + x_nums[i] /= g; + x_dens[i] /= g; + } + } + + // Check x >= 0: each x_nums[i] / x_dens[i] >= 0 + // Since x_dens[i] > 0 (normalized), we need x_nums[i] >= 0 + x_nums.iter().take(m).all(|&num| num >= 0) + } +} + +/// Compute GCD of two i128 values. +fn gcd_i128(mut a: i128, mut b: i128) -> i128 { + while b != 0 { + let t = b; + b = a % b; + a = t; + } + a +} + +impl Problem for FeasibleBasisExtension { + const NAME: &'static str = "FeasibleBasisExtension"; + type Value = crate::types::Or; + + fn variant() -> Vec<(&'static str, &'static str)> { + crate::variant_params![] + } + + fn dims(&self) -> Vec { + vec![2; self.num_columns() - self.num_required()] + } + + fn evaluate(&self, config: &[usize]) -> crate::types::Or { + let free_cols = self.free_columns(); + let num_free = free_cols.len(); + + if config.len() != num_free { + return crate::types::Or(false); + } + if config.iter().any(|&v| v >= 2) { + return crate::types::Or(false); + } + + let m = self.num_rows(); + let s = self.num_required(); + let needed = m - s; + + // Count selected free columns + let selected_free: Vec = config + .iter() + .enumerate() + .filter(|(_, &v)| v == 1) + .map(|(i, _)| free_cols[i]) + .collect(); + + if selected_free.len() != needed { + return crate::types::Or(false); + } + + // Form basis: required columns + selected free columns + let mut basis_cols: Vec = self.required_columns.clone(); + basis_cols.extend_from_slice(&selected_free); + + crate::types::Or(self.check_feasible_basis(&basis_cols)) + } +} + +crate::declare_variants! { + default FeasibleBasisExtension => "2^num_columns * num_rows^3", +} + +#[cfg(feature = "example-db")] +pub(crate) fn canonical_model_example_specs() -> Vec { + vec![crate::example_db::specs::ModelExampleSpec { + id: "feasible_basis_extension", + // 3x6 matrix, rhs=[7,5,3], required={0,1}, select col 2 -> B={0,1,2}, x=(4,5,3)>=0 + instance: Box::new(FeasibleBasisExtension::new( + vec![ + vec![1, 0, 1, 2, -1, 0], + vec![0, 1, 0, 1, 1, 2], + vec![0, 0, 1, 1, 0, 1], + ], + vec![7, 5, 3], + vec![0, 1], + )), + optimal_config: vec![1, 0, 0, 0], // select col 2 (first free column) + optimal_value: serde_json::json!(true), + }] +} + +#[cfg(test)] +#[path = "../../unit_tests/models/algebraic/feasible_basis_extension.rs"] +mod tests; diff --git a/src/models/algebraic/mod.rs b/src/models/algebraic/mod.rs index d341270c4..134662fd4 100644 --- a/src/models/algebraic/mod.rs +++ b/src/models/algebraic/mod.rs @@ -15,6 +15,7 @@ pub(crate) mod closest_vector_problem; pub(crate) mod consecutive_block_minimization; pub(crate) mod consecutive_ones_matrix_augmentation; pub(crate) mod consecutive_ones_submatrix; +pub(crate) mod feasible_basis_extension; pub(crate) mod ilp; pub(crate) mod quadratic_assignment; pub(crate) mod qubo; @@ -25,6 +26,7 @@ pub use closest_vector_problem::{ClosestVectorProblem, VarBounds}; pub use consecutive_block_minimization::ConsecutiveBlockMinimization; pub use consecutive_ones_matrix_augmentation::ConsecutiveOnesMatrixAugmentation; pub use consecutive_ones_submatrix::ConsecutiveOnesSubmatrix; +pub use feasible_basis_extension::FeasibleBasisExtension; pub use ilp::{Comparison, LinearConstraint, ObjectiveSense, VariableDomain, ILP}; pub use quadratic_assignment::QuadraticAssignment; pub use qubo::QUBO; @@ -40,6 +42,7 @@ pub(crate) fn canonical_model_example_specs() -> Vec FeasibleBasisExtension { + FeasibleBasisExtension::new( + vec![ + vec![1, 0, 1, 2, -1, 0], + vec![0, 1, 0, 1, 1, 2], + vec![0, 0, 1, 1, 0, 1], + ], + vec![7, 5, 3], + vec![0, 1], + ) +} + +#[test] +fn test_feasible_basis_extension_creation() { + let problem = issue_example(); + assert_eq!(problem.num_rows(), 3); + assert_eq!(problem.num_columns(), 6); + assert_eq!(problem.num_required(), 2); + assert_eq!(problem.dims(), vec![2; 4]); // 6 - 2 = 4 free columns + assert_eq!( + ::NAME, + "FeasibleBasisExtension" + ); + assert_eq!(::variant(), vec![]); +} + +#[test] +fn test_feasible_basis_extension_evaluate_satisfying() { + let problem = issue_example(); + // Free columns are [2, 3, 4, 5]. Select col 2 (index 0 in free list). + // B = {0, 1, 2}. A_B = I_3. x = (7, 5, 3) but actually A_B = [[1,0,1],[0,1,0],[0,0,1]] + // A_B^{-1} a_bar: solve [[1,0,1],[0,1,0],[0,0,1]] x = [7,5,3] => x = (4, 5, 3) >= 0 + assert!(problem.evaluate(&[1, 0, 0, 0])); +} + +#[test] +fn test_feasible_basis_extension_evaluate_satisfying_col3() { + let problem = issue_example(); + // Select col 3 (index 1 in free list). B = {0, 1, 3}. + // A_B = [[1,0,2],[0,1,1],[0,0,1]]. Solve: x = (1, 2, 3) >= 0 + assert!(problem.evaluate(&[0, 1, 0, 0])); +} + +#[test] +fn test_feasible_basis_extension_evaluate_singular() { + let problem = issue_example(); + // Select col 4 (index 2 in free list). B = {0, 1, 4}. + // A_B = [[1,0,-1],[0,1,1],[0,0,0]] => singular + assert!(!problem.evaluate(&[0, 0, 1, 0])); +} + +#[test] +fn test_feasible_basis_extension_evaluate_infeasible_negative() { + let problem = issue_example(); + // Select col 5 (index 3 in free list). B = {0, 1, 5}. + // A_B = [[1,0,0],[0,1,2],[0,0,1]]. Solve: x = (7, -1, 3). x_1 = -1 < 0 + assert!(!problem.evaluate(&[0, 0, 0, 1])); +} + +#[test] +fn test_feasible_basis_extension_evaluate_wrong_count() { + let problem = issue_example(); + // Need exactly 1 column selected (m - |S| = 3 - 2 = 1) + assert!(!problem.evaluate(&[1, 1, 0, 0])); // too many + assert!(!problem.evaluate(&[0, 0, 0, 0])); // too few +} + +#[test] +fn test_feasible_basis_extension_evaluate_wrong_config_length() { + let problem = issue_example(); + assert!(!problem.evaluate(&[1, 0])); // too short + assert!(!problem.evaluate(&[1, 0, 0, 0, 0])); // too long +} + +#[test] +fn test_feasible_basis_extension_evaluate_invalid_variable_value() { + let problem = issue_example(); + assert!(!problem.evaluate(&[2, 0, 0, 0])); +} + +#[test] +fn test_feasible_basis_extension_brute_force() { + let problem = issue_example(); + let solver = BruteForce::new(); + let solution = solver + .find_witness(&problem) + .expect("should find a solution"); + assert!(problem.evaluate(&solution)); +} + +#[test] +fn test_feasible_basis_extension_brute_force_all() { + let problem = issue_example(); + let solver = BruteForce::new(); + let solutions = solver.find_all_witnesses(&problem); + // From the issue: B={0,1,2} and B={0,1,3} are feasible, so 2 solutions + assert_eq!(solutions.len(), 2); + for sol in &solutions { + assert!(problem.evaluate(sol)); + } +} + +#[test] +fn test_feasible_basis_extension_unsatisfiable() { + // Construct an instance where no feasible basis exists. + // 2x3 matrix, rhs=[1,-1], S={} (need 2 columns from 3). + // A = [[1,0,1],[0,1,-1]], rhs = [1,-1] + // B={0,1}: solve [[1,0],[0,1]]x=[1,-1] => x=(1,-1), x_1<0 + // B={0,2}: solve [[1,1],[0,-1]]x=[1,-1] => x=(0,1), x>=0 => feasible! + // Let's try: A = [[1,1,1],[1,1,1]], rhs = [1,1]. All 2x2 submatrices are singular. + let problem = + FeasibleBasisExtension::new(vec![vec![1, 1, 1], vec![1, 1, 1]], vec![1, 1], vec![]); + let solver = BruteForce::new(); + assert!(solver.find_witness(&problem).is_none()); +} + +#[test] +fn test_feasible_basis_extension_serialization() { + let problem = issue_example(); + let json = serde_json::to_value(&problem).unwrap(); + assert_eq!( + json, + serde_json::json!({ + "matrix": [ + [1, 0, 1, 2, -1, 0], + [0, 1, 0, 1, 1, 2], + [0, 0, 1, 1, 0, 1], + ], + "rhs": [7, 5, 3], + "required_columns": [0, 1], + }) + ); + let restored: FeasibleBasisExtension = serde_json::from_value(json).unwrap(); + assert_eq!(restored.num_rows(), 3); + assert_eq!(restored.num_columns(), 6); + assert_eq!(restored.num_required(), 2); +} + +#[test] +fn test_feasible_basis_extension_paper_example() { + let problem = issue_example(); + // Verify B={0,1,2} is satisfying (config [1,0,0,0]) + assert!(problem.evaluate(&[1, 0, 0, 0])); + // Verify B={0,1,3} is satisfying (config [0,1,0,0]) + assert!(problem.evaluate(&[0, 1, 0, 0])); + + let solver = BruteForce::new(); + let solutions = solver.find_all_witnesses(&problem); + assert_eq!(solutions.len(), 2); +} + +#[test] +fn test_feasible_basis_extension_complexity_metadata() { + use crate::registry::VariantEntry; + + let entry = inventory::iter::() + .find(|entry| entry.name == "FeasibleBasisExtension") + .expect("FeasibleBasisExtension variant entry should exist"); + + assert_eq!(entry.complexity, "2^num_columns * num_rows^3"); +} + +#[test] +#[should_panic(expected = "must be less than")] +fn test_feasible_basis_extension_m_ge_n() { + // 3x3 matrix: m not < n + FeasibleBasisExtension::new( + vec![vec![1, 0, 0], vec![0, 1, 0], vec![0, 0, 1]], + vec![1, 1, 1], + vec![], + ); +} + +#[test] +#[should_panic(expected = "rhs length")] +fn test_feasible_basis_extension_rhs_length_mismatch() { + FeasibleBasisExtension::new( + vec![vec![1, 0, 1], vec![0, 1, 0]], + vec![1, 2, 3], // length 3, but m=2 + vec![], + ); +} + +#[test] +#[should_panic(expected = "|S|")] +fn test_feasible_basis_extension_too_many_required() { + // m=2, |S|=2 is not < m + FeasibleBasisExtension::new(vec![vec![1, 0, 1], vec![0, 1, 0]], vec![1, 2], vec![0, 1]); +} + +#[test] +#[should_panic(expected = "out of bounds")] +fn test_feasible_basis_extension_required_out_of_bounds() { + FeasibleBasisExtension::new( + vec![vec![1, 0, 1], vec![0, 1, 0]], + vec![1, 2], + vec![5], // out of bounds + ); +} + +#[test] +#[should_panic(expected = "Duplicate")] +fn test_feasible_basis_extension_duplicate_required() { + // 3x5 matrix so |S|=2 < m=3, but S has duplicates + FeasibleBasisExtension::new( + vec![ + vec![1, 0, 1, 0, 1], + vec![0, 1, 0, 1, 0], + vec![1, 1, 0, 0, 1], + ], + vec![1, 2, 3], + vec![0, 0], + ); +}