problemreductions/models/algebraic/

feasible_basis_extension.rs

//! 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::{CreateSpec, ProblemSchemaEntry};
use crate::traits::Problem;
use serde::{Deserialize, Serialize};

inventory::submit! {
    ProblemSchemaEntry {
        name: "FeasibleBasisExtension",
        display_name: "Feasible Basis Extension",
        aliases: &[],
        dimensions: &[],
        category: crate::registry::ProblemCategory::Algebraic,
        module_path: module_path!(),
        description: "Given matrix A, vector a_bar, and required columns S, find a feasible basis extending S",
        fields: FeasibleBasisExtensionCreateSpec::FIELDS,
    }
}

/// 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, 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.solve(&problem).unwrap();
/// assert!(solution.is_some());
/// ```
#[derive(Debug, Clone, Serialize, Deserialize)]
#[serde(try_from = "FeasibleBasisExtensionCreateSpec")]
pub struct FeasibleBasisExtension {
    matrix: Vec<Vec<i64>>,
    rhs: Vec<i64>,
    required_columns: Vec<usize>,
}

#[derive(Debug, Deserialize, crate::CreateSpec)]
struct FeasibleBasisExtensionCreateSpec {
    /// Integer matrix as JSON.
    #[create(codec = "json")]
    matrix: Vec<Vec<i64>>,
    /// Right-hand side vector.
    #[create(codec = "comma-separated")]
    rhs: Vec<i64>,
    /// Required column indices.
    #[create(codec = "comma-separated")]
    required_columns: Vec<usize>,
}

impl TryFrom<FeasibleBasisExtensionCreateSpec> for FeasibleBasisExtension {
    type Error = crate::registry::ConstructionError;

    fn try_from(spec: FeasibleBasisExtensionCreateSpec) -> Result<Self, Self::Error> {
        Self::try_new(spec.matrix, spec.rhs, spec.required_columns)
    }
}

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<Vec<i64>>, rhs: Vec<i64>, required_columns: Vec<usize>) -> Self {
        Self::try_new(matrix, rhs, required_columns).unwrap_or_else(|error| panic!("{error}"))
    }

    fn try_new(
        matrix: Vec<Vec<i64>>,
        rhs: Vec<i64>,
        required_columns: Vec<usize>,
    ) -> Result<Self, crate::registry::ConstructionError> {
        let m = matrix.len();
        let first = matrix.first().ok_or("matrix must have at least one row")?;
        let n = first.len();
        if matrix.iter().any(|row| row.len() != n) {
            return Err("all matrix rows must have the same length".into());
        }
        if m >= n {
            return Err("number of rows must be less than number of columns".into());
        }
        if rhs.len() != m {
            return Err("rhs length must equal number of rows".into());
        }
        if required_columns.len() >= m {
            return Err(
                format!("|S| ({}) must be less than m ({m})", required_columns.len()).into(),
            );
        }
        let mut seen = std::collections::HashSet::new();
        for &column in &required_columns {
            if column >= n {
                return Err(format!("required column {column} is out of bounds").into());
            }
            if !seen.insert(column) {
                return Err(format!("Duplicate required column index {column}").into());
            }
        }
        Ok(Self {
            matrix,
            rhs,
            required_columns,
        })
    }

    /// Returns the matrix A.
    pub fn matrix(&self) -> &[Vec<i64>] {
        &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<usize> {
        let required_set: std::collections::HashSet<usize> =
            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],
    ) -> Result<bool, crate::traits::EvaluationError> {
        let m = self.num_rows();
        assert_eq!(basis_cols.len(), m);

        // Build augmented matrix [A_B | a_bar] for Bareiss elimination.
        let mut augmented: Vec<Vec<i64>> = 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]);
            }
            row.push(self.rhs[i]);
            augmented.push(row);
        }

        // Bareiss algorithm: fraction-free Gaussian elimination.
        // After elimination, the system is upper-triangular.
        let mut prev_pivot = 1_i64;

        for k in 0..m {
            // Partial pivoting
            let mut max_row = k;
            let mut max_val = augmented[k][k].checked_abs().ok_or_else(|| {
                crate::traits::EvaluationError::IntegerOverflow(
                    "taking an elimination pivot magnitude".into(),
                )
            })?;
            for i in (k + 1)..m {
                let candidate = augmented[i][k].checked_abs().ok_or_else(|| {
                    crate::traits::EvaluationError::IntegerOverflow(
                        "taking an elimination candidate magnitude".into(),
                    )
                })?;
                if candidate > max_val {
                    max_val = candidate;
                    max_row = i;
                }
            }
            if max_val == 0 {
                return Ok(false); // singular
            }
            if max_row != k {
                augmented.swap(k, max_row);
            }

            for i in (k + 1)..m {
                for j in (k + 1)..=m {
                    let left = augmented[k][k]
                        .checked_mul(augmented[i][j])
                        .ok_or_else(|| {
                            crate::traits::EvaluationError::IntegerOverflow(
                                "multiplying Bareiss pivot and row entry".into(),
                            )
                        })?;
                    let right = augmented[i][k]
                        .checked_mul(augmented[k][j])
                        .ok_or_else(|| {
                            crate::traits::EvaluationError::IntegerOverflow(
                                "multiplying Bareiss elimination entries".into(),
                            )
                        })?;
                    let numerator = left.checked_sub(right).ok_or_else(|| {
                        crate::traits::EvaluationError::IntegerOverflow(
                            "subtracting Bareiss products".into(),
                        )
                    })?;
                    augmented[i][j] = numerator.checked_div(prev_pivot).ok_or_else(|| {
                        crate::traits::EvaluationError::IntegerOverflow(
                            "dividing by the previous Bareiss pivot".into(),
                        )
                    })?;
                }
                augmented[i][k] = 0;
            }
            prev_pivot = augmented[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![0_i64; m];
        let mut x_dens = vec![1_i64; m];

        for i in (0..m).rev() {
            // numerator of (aug128[i][m] - sum_{j>i} aug128[i][j] * x[j])
            let mut num = augmented[i][m];
            let mut den = 1_i64;

            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 left = num.checked_mul(x_dens[j]).ok_or_else(|| {
                    crate::traits::EvaluationError::IntegerOverflow(
                        "multiplying basis-solution numerator".into(),
                    )
                })?;
                let right = den
                    .checked_mul(augmented[i][j])
                    .and_then(|value| value.checked_mul(x_nums[j]))
                    .ok_or_else(|| {
                        crate::traits::EvaluationError::IntegerOverflow(
                            "multiplying basis-solution subtraction term".into(),
                        )
                    })?;
                num = left.checked_sub(right).ok_or_else(|| {
                    crate::traits::EvaluationError::IntegerOverflow(
                        "subtracting basis-solution terms".into(),
                    )
                })?;
                den = den.checked_mul(x_dens[j]).ok_or_else(|| {
                    crate::traits::EvaluationError::IntegerOverflow(
                        "multiplying basis-solution denominators".into(),
                    )
                })?;
                // Simplify to avoid overflow
                let g = gcd_i64(
                    num.checked_abs().ok_or_else(|| {
                        crate::traits::EvaluationError::IntegerOverflow(
                            "taking basis-solution numerator magnitude".into(),
                        )
                    })?,
                    den.checked_abs().ok_or_else(|| {
                        crate::traits::EvaluationError::IntegerOverflow(
                            "taking basis-solution denominator magnitude".into(),
                        )
                    })?,
                );
                if g > 1 {
                    num /= g;
                    den /= g;
                }
            }
            // x[i] = (num/den) / aug128[i][i] = num / (den * aug128[i][i])
            let diag = augmented[i][i];
            x_nums[i] = num;
            x_dens[i] = den.checked_mul(diag).ok_or_else(|| {
                crate::traits::EvaluationError::IntegerOverflow(
                    "multiplying basis-solution denominator by diagonal".into(),
                )
            })?;
            // Normalize sign: make denominator positive
            if x_dens[i] < 0 {
                x_nums[i] = x_nums[i].checked_neg().ok_or_else(|| {
                    crate::traits::EvaluationError::IntegerOverflow(
                        "normalizing basis-solution numerator sign".into(),
                    )
                })?;
                x_dens[i] = x_dens[i].checked_neg().ok_or_else(|| {
                    crate::traits::EvaluationError::IntegerOverflow(
                        "normalizing basis-solution denominator sign".into(),
                    )
                })?;
            }
            let g = gcd_i64(
                x_nums[i].checked_abs().ok_or_else(|| {
                    crate::traits::EvaluationError::IntegerOverflow(
                        "taking normalized numerator magnitude".into(),
                    )
                })?,
                x_dens[i].checked_abs().ok_or_else(|| {
                    crate::traits::EvaluationError::IntegerOverflow(
                        "taking normalized denominator magnitude".into(),
                    )
                })?,
            );
            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
        Ok(x_nums.iter().take(m).all(|&num| num >= 0))
    }
}

/// Compute GCD of two nonnegative i64 values.
fn gcd_i64(mut a: i64, mut b: i64) -> i64 {
    while b != 0 {
        let t = b;
        b = a % b;
        a = t;
    }
    a
}

impl Problem for FeasibleBasisExtension {
    const NAME: &'static str = "FeasibleBasisExtension";
    type Solution = Vec<bool>;
    type Value = crate::types::Or;

    crate::problem_parameters![
        ("num_rows", num_rows),
        ("num_columns", num_columns),
        ("num_required", num_required),
    ];

    fn variant() -> Vec<(&'static str, &'static str)> {
        crate::variant_params![]
    }

    fn evaluate(
        &self,
        config: &Self::Solution,
    ) -> Result<crate::types::Or, crate::traits::EvaluationError> {
        Ok({
            let free_cols = self.free_columns();
            let num_free = free_cols.len();

            if config.len() != num_free {
                return Err(crate::traits::EvaluationError::InvalidConfiguration(
                    "free-column selection length does not match the matrix".into(),
                ));
            }
            let m = self.num_rows();
            let s = self.num_required();
            let needed = m - s;

            // Count selected free columns
            let selected_free: Vec<usize> = config
                .iter()
                .enumerate()
                .filter(|(_, &v)| v)
                .map(|(i, _)| free_cols[i])
                .collect();

            if selected_free.len() != needed {
                return Ok(crate::types::Or(false));
            }

            // Form basis: required columns + selected free columns
            let mut basis_cols: Vec<usize> = self.required_columns.clone();
            basis_cols.extend_from_slice(&selected_free);

            crate::types::Or(self.check_feasible_basis(&basis_cols)?)
        })
    }
}

impl crate::solvers::BruteForceProblem for FeasibleBasisExtension {
    fn dimensions(&self) -> Vec<usize> {
        vec![2; self.num_columns() - self.num_required()]
    }
}

crate::declare_variants! {
    default FeasibleBasisExtension => "2^num_columns * num_rows^3" create FeasibleBasisExtensionCreateSpec,
}

crate::register_brute_force! {
    FeasibleBasisExtension decode |_, indices: Vec<usize>| crate::config::config_to_bits(&indices),
}

#[cfg(feature = "example-db")]
pub(crate) fn canonical_model_example_specs() -> Vec<crate::example_db::specs::ModelExampleSpec> {
    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: serde_json::json!(vec![true, false, false, false]), // 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;