problemreductions/rules/

ksatisfiability_qubo.rs

//! Reduction from KSatisfiability to QUBO (Max-K-SAT).
//!
//! For K=2 (quadratic penalty), each clause contributes to Q based on literal signs:
//! - (x_i ∨ x_j): penalty (1-x_i)(1-x_j) → Q[i][i]-=1, Q[j][j]-=1, Q[i][j]+=1, const+=1
//! - (¬x_i ∨ x_j): penalty x_i(1-x_j) → Q[i][i]+=1, Q[i][j]-=1
//! - (x_i ∨ ¬x_j): penalty (1-x_i)x_j → Q[j][j]+=1, Q[i][j]-=1
//! - (¬x_i ∨ ¬x_j): penalty x_i·x_j → Q[i][j]+=1
//!
//! For K≥3, we use the Rosenberg quadratization to reduce degree-K penalty terms
//! to quadratic form by introducing auxiliary variables. Each clause of K literals
//! requires K−2 auxiliary variables.
//!
//! CNFClause uses 1-indexed signed integers: positive = variable, negative = negated.

use crate::models::algebraic::QUBO;
use crate::models::decision::Decision;
use crate::models::formula::KSatisfiability;
use crate::reduction;
use crate::rules::traits::{ReduceTo, ReductionResult};
use crate::variant::{K2, K3};
/// Result of reducing KSatisfiability to QUBO.
#[derive(Debug, Clone)]
pub struct ReductionKSatToQUBO {
    target: Decision<QUBO<i64>>,
    source_num_vars: usize,
}

impl ReductionResult for ReductionKSatToQUBO {
    type Source = KSatisfiability<K2>;
    type Target = Decision<QUBO<i64>>;

    fn target_problem(&self) -> &Self::Target {
        &self.target
    }

    fn extract_solution(
        &self,
        target_solution: &<Self::Target as crate::traits::Problem>::Solution,
    ) -> crate::rules::ExtractionResult<<Self::Source as crate::traits::Problem>::Solution> {
        crate::rules::traits::validate_target_witness(
            self.target_problem(),
            target_solution,
            |value| value.0,
            "QUBO energy does not meet the SAT zero-penalty threshold",
        )?;
        Ok(target_solution[..self.source_num_vars].to_vec())
    }
}

/// Result of reducing `KSatisfiability<K3>` to QUBO.
#[derive(Debug, Clone)]
pub struct Reduction3SATToQUBO {
    target: Decision<QUBO<i64>>,
    source_num_vars: usize,
}

impl ReductionResult for Reduction3SATToQUBO {
    type Source = KSatisfiability<K3>;
    type Target = Decision<QUBO<i64>>;

    fn target_problem(&self) -> &Self::Target {
        &self.target
    }

    fn extract_solution(
        &self,
        target_solution: &<Self::Target as crate::traits::Problem>::Solution,
    ) -> crate::rules::ExtractionResult<<Self::Source as crate::traits::Problem>::Solution> {
        crate::rules::traits::validate_target_witness(
            self.target_problem(),
            target_solution,
            |value| value.0,
            "QUBO energy does not meet the SAT zero-penalty threshold",
        )?;
        Ok(target_solution[..self.source_num_vars].to_vec())
    }
}

/// Add the quadratic penalty term for a 2-literal clause to the QUBO matrix.
///
/// For clause (l_i ∨ l_j), the penalty for the clause being unsatisfied is
/// the product of the complemented literals.
fn add_coefficient(
    matrix: &mut [Vec<i64>],
    row: usize,
    column: usize,
    coefficient: i64,
) -> Result<(), &'static str> {
    matrix[row][column] = matrix[row][column]
        .checked_add(coefficient)
        .ok_or("adding a SAT QUBO coefficient")?;
    Ok(())
}

fn add_2sat_clause_penalty(
    matrix: &mut [Vec<i64>],
    lits: &[(usize, bool)],
) -> Result<(), &'static str> {
    assert_eq!(lits.len(), 2, "Expected 2-literal clause");

    let (var_i, neg_i) = lits[0];
    let (var_j, neg_j) = lits[1];

    // Ensure i <= j for upper-triangular form
    let (i, j, ni, nj) = if var_i <= var_j {
        (var_i, var_j, neg_i, neg_j)
    } else {
        (var_j, var_i, neg_j, neg_i)
    };

    match (ni, nj) {
        (false, false) => {
            // (x_i ∨ x_j): penalty = (1-x_i)(1-x_j) = 1 - x_i - x_j + x_i·x_j
            add_coefficient(matrix, i, i, -1)?;
            add_coefficient(matrix, j, j, -1)?;
            add_coefficient(matrix, i, j, 1)?;
        }
        (true, false) => {
            // (¬x_i ∨ x_j): penalty = x_i(1-x_j) = x_i - x_i·x_j
            add_coefficient(matrix, i, i, 1)?;
            add_coefficient(matrix, i, j, -1)?;
        }
        (false, true) => {
            // (x_i ∨ ¬x_j): penalty = (1-x_i)x_j = x_j - x_i·x_j
            add_coefficient(matrix, j, j, 1)?;
            add_coefficient(matrix, i, j, -1)?;
        }
        (true, true) => {
            // (¬x_i ∨ ¬x_j): penalty = x_i·x_j
            add_coefficient(matrix, i, j, 1)?;
        }
    }
    Ok(())
}

/// Add the QUBO terms for a 3-literal clause using Rosenberg quadratization.
///
/// For clause (l1 ∨ l2 ∨ l3), the penalty for not satisfying the clause is:
///   P = (1-l1)(1-l2)(1-l3) = y1·y2·y3
/// where yi = 1 - li (complement of literal).
///
/// We introduce one auxiliary variable `a` and quadratize using the substitution
/// a = y1·y2, adding penalty M·(y1·y2 - 2·y1·a - 2·y2·a + 3·a) where M is a
/// sufficiently large penalty (M = 2 suffices for Max-SAT).
///
/// The resulting quadratic form is:
///   H = a·y3 + M·(y1·y2 - 2·y1·a - 2·y2·a + 3·a)
///
/// `aux_var` is the 0-indexed auxiliary variable.
fn add_3sat_clause_penalty(
    matrix: &mut [Vec<i64>],
    lits: &[(usize, bool)],
    aux_var: usize,
) -> Result<(), &'static str> {
    assert_eq!(lits.len(), 3, "Expected 3-literal clause");
    let penalty = 2; // Rosenberg penalty weight

    let (v1, n1) = lits[0];
    let (v2, n2) = lits[1];
    let (v3, n3) = lits[2];
    let a = aux_var;

    // We need to express yi = (1 - li) in terms of xi:
    //   If literal is positive (li = xi): yi = 1 - xi
    //   If literal is negated (li = 1 - xi): yi = xi
    //
    // So yi = xi if negated, yi = 1 - xi if positive.
    //
    // We compute the QUBO terms for:
    //   H = a·y3 + M·(y1·y2 - 2·y1·a - 2·y2·a + 3·a)
    //
    // Each term is expanded using yi = xi (if negated) or yi = 1-xi (if positive).

    // Helper: add coefficient * yi * yj to the matrix
    // where yi depends on variable vi and negation ni
    let add_yy = |matrix: &mut [Vec<i64>],
                  vi: usize,
                  ni: bool,
                  vj: usize,
                  nj: bool,
                  coeff: i64|
     -> Result<(), &'static str> {
        // yi = xi if ni (negated literal), yi = 1 - xi if !ni (positive literal)
        // yi * yj expansion:
        if vi == vj {
            // Same variable: yi * yj
            // Both complemented the same way means yi = yj, so yi*yj = yi (binary)
            // If ni == nj: yi*yj = yi^2 = yi (binary), add coeff * yi
            // If ni != nj: yi * yj = xi * (1-xi) = 0 (always), add nothing
            if ni == nj {
                // yi * yi = yi (binary)
                if ni {
                    // yi = xi, add coeff * xi
                    add_coefficient(matrix, vi, vi, coeff)?;
                } else {
                    // yi = 1 - xi, add coeff * (1 - xi) = coeff - coeff * xi
                    // constant term ignored in QUBO (offset), diagonal:
                    add_coefficient(matrix, vi, vi, -coeff)?;
                }
            }
            // else: xi * (1-xi) = 0, nothing to add
            return Ok(());
        }
        // Different variables: yi * yj
        let (lo, hi, lo_neg, hi_neg) = if vi < vj {
            (vi, vj, ni, nj)
        } else {
            (vj, vi, nj, ni)
        };
        // yi = xi if neg, else 1-xi
        // yj = xj if neg, else 1-xj
        // yi*yj = (xi if neg_i else 1-xi) * (xj if neg_j else 1-xj)
        match (lo_neg, hi_neg) {
            (true, true) => {
                // xi * xj
                add_coefficient(matrix, lo, hi, coeff)?;
            }
            (true, false) => {
                // xi * (1 - xj) = xi - xi*xj
                add_coefficient(matrix, lo, lo, coeff)?;
                add_coefficient(matrix, lo, hi, -coeff)?;
            }
            (false, true) => {
                // (1 - xi) * xj = xj - xi*xj
                add_coefficient(matrix, hi, hi, coeff)?;
                add_coefficient(matrix, lo, hi, -coeff)?;
            }
            (false, false) => {
                // (1-xi)(1-xj) = 1 - xi - xj + xi*xj
                // constant 1 ignored (offset)
                add_coefficient(matrix, lo, lo, -coeff)?;
                add_coefficient(matrix, hi, hi, -coeff)?;
                add_coefficient(matrix, lo, hi, coeff)?;
            }
        }
        Ok(())
    };

    // Helper: add coefficient * yi * a to the matrix
    // where yi depends on variable vi and negation ni, a is aux variable
    let add_ya = |matrix: &mut [Vec<i64>],
                  vi: usize,
                  ni: bool,
                  a: usize,
                  coeff: i64|
     -> Result<(), &'static str> {
        // yi = xi if ni (negated literal), yi = 1-xi if !ni (positive literal)
        // yi * a:
        let (lo, hi) = if vi < a { (vi, a) } else { (a, vi) };
        if ni {
            // yi = xi, so yi * a = xi * a
            add_coefficient(matrix, lo, hi, coeff)?;
        } else {
            // yi = 1 - xi, so yi * a = a - xi * a
            add_coefficient(matrix, a, a, coeff)?;
            add_coefficient(matrix, lo, hi, -coeff)?;
        }
        Ok(())
    };

    // Term 1: a * y3 (coefficient = 1)
    add_ya(matrix, v3, n3, a, 1)?;

    // Term 2: M * y1 * y2
    add_yy(matrix, v1, n1, v2, n2, penalty)?;

    // Term 3: -2M * y1 * a
    add_ya(matrix, v1, n1, a, -2 * penalty)?;

    // Term 4: -2M * y2 * a
    add_ya(matrix, v2, n2, a, -2 * penalty)?;

    // Term 5: 3M * a (linear)
    // a is a binary variable, a^2 = a, so linear a → diagonal
    add_coefficient(matrix, a, a, 3 * penalty)?;

    Ok(())
}

/// Expand clause penalties and retain the constant omitted by QUBO.
/// K3 reserves one auxiliary per clause, including free auxiliaries for short
/// clauses; K2 reserves none. The source constructor validates clause widths.
fn build_qubo_matrix(
    num_vars: usize,
    clauses: &[crate::models::formula::CNFClause],
    num_aux: usize,
) -> Result<(Vec<Vec<i64>>, i64), &'static str> {
    let total = num_vars
        .checked_add(num_aux)
        .ok_or("computing the number of SAT QUBO variables")?;
    total
        .checked_mul(total)
        .ok_or("computing the SAT QUBO dense matrix entry count")?;
    let mut matrix = vec![vec![0; total]; total];
    let mut constant = 0i64;
    for (idx, clause) in clauses.iter().enumerate() {
        let literals: Vec<_> = clause
            .variables()
            .into_iter()
            .zip(&clause.literals)
            .map(|(variable, &literal)| (variable, literal < 0))
            .collect();
        let offset = match literals.as_slice() {
            [] => 1,
            &[(v, neg)] => {
                add_coefficient(&mut matrix, v, v, if neg { 1 } else { -1 })?;
                i64::from(!neg)
            }
            &[(_, n1), (_, n2)] => {
                add_2sat_clause_penalty(&mut matrix, &literals)?;
                i64::from(!n1 && !n2)
            }
            &[(_, n1), (_, n2), _] => {
                add_3sat_clause_penalty(&mut matrix, &literals, num_vars + idx)?;
                2 * i64::from(!n1 && !n2)
            }
            _ => unreachable!("the source validates clause width at most three"),
        };
        constant = constant
            .checked_add(offset)
            .ok_or("accumulating the SAT QUBO constant")?;
    }
    Ok((matrix, constant))
}

#[crate::aggregate_reduction(identity)]
impl crate::rules::AggregateReductionResult for ReductionKSatToQUBO {}

#[crate::aggregate_reduction(identity)]
impl crate::rules::AggregateReductionResult for Reduction3SATToQUBO {}

#[reduction(
    transform = {
        exact { num_vars = "num_vars" },
        upper_bound { num_quadratic_terms = "num_vars * (num_vars - 1) / 2" },
    }
)]
impl ReduceTo<Decision<QUBO<i64>>> for KSatisfiability<K2> {
    type Result = ReductionKSatToQUBO;

    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
        let n = self.num_vars();
        let (matrix, constant) =
            build_qubo_matrix(n, self.clauses(), 0).map_err(|operation| {
                crate::rules::ReductionError::integer_overflow::<
                    KSatisfiability<K2>,
                    Decision<QUBO<i64>>,
                >(operation)
            })?;

        Ok(ReductionKSatToQUBO {
            target: Decision::new(
                QUBO::from_matrix(matrix).map_err(|message| {
                    crate::rules::ReductionError::construction::<
                        KSatisfiability<K2>,
                        Decision<QUBO<i64>>,
                    >(message)
                })?,
                -constant,
            ),
            source_num_vars: n,
        })
    }
}

#[reduction(transform = {
    exact {
        num_vars = "num_vars + num_clauses",
    },
    upper_bound {
        num_quadratic_terms = "(num_vars + num_clauses) * ((num_vars + num_clauses) - 1) / 2",
    },
})]
impl ReduceTo<Decision<QUBO<i64>>> for KSatisfiability<K3> {
    type Result = Reduction3SATToQUBO;

    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
        let n = self.num_vars();
        let (matrix, constant) =
            build_qubo_matrix(n, self.clauses(), self.num_clauses()).map_err(|operation| {
                crate::rules::ReductionError::integer_overflow::<
                    KSatisfiability<K3>,
                    Decision<QUBO<i64>>,
                >(operation)
            })?;

        Ok(Reduction3SATToQUBO {
            target: Decision::new(
                QUBO::from_matrix(matrix).map_err(|message| {
                    crate::rules::ReductionError::construction::<
                        KSatisfiability<K3>,
                        Decision<QUBO<i64>>,
                    >(message)
                })?,
                -constant,
            ),
            source_num_vars: n,
        })
    }
}

#[cfg(feature = "example-db")]
pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
    use crate::export::SolutionPair;
    use crate::models::algebraic::QUBO;
    use crate::models::formula::CNFClause;

    vec![
        crate::example_db::specs::RuleExampleSpec {
            id: "ksatisfiability_k2_to_qubo",
            build: || {
                let source = KSatisfiability::<K2>::new(
                    4,
                    vec![
                        CNFClause::new(vec![1, 2]),
                        CNFClause::new(vec![-1, 3]),
                        CNFClause::new(vec![-2, 4]),
                        CNFClause::new(vec![-3, -4]),
                    ],
                );
                crate::example_db::specs::rule_example_with_witness::<_, Decision<QUBO<i64>>>(
                    source,
                    SolutionPair {
                        source_config: serde_json::json!(vec![false, true, false, true]),
                        target_config: serde_json::json!(vec![false, true, false, true]),
                    },
                )
            },
        },
        crate::example_db::specs::RuleExampleSpec {
            id: "ksatisfiability_to_qubo",
            build: || {
                let source = KSatisfiability::<K3>::new(
                    5,
                    vec![
                        CNFClause::new(vec![1, 2, -3]),
                        CNFClause::new(vec![-1, 3, 4]),
                        CNFClause::new(vec![2, -4, 5]),
                        CNFClause::new(vec![-2, 3, -5]),
                        CNFClause::new(vec![1, -3, 5]),
                        CNFClause::new(vec![-1, -2, 4]),
                        CNFClause::new(vec![3, -4, -5]),
                    ],
                );
                crate::example_db::specs::rule_example_with_witness::<_, Decision<QUBO<i64>>>(
                    source,
                    SolutionPair {
                        source_config: serde_json::json!(vec![false, false, false, false, false]),
                        target_config: serde_json::json!(vec![
                            false, false, false, false, false, true, false, false, false, false,
                            false, false
                        ]),
                    },
                )
            },
        },
    ]
}

#[cfg(test)]
#[path = "../unit_tests/rules/ksatisfiability_qubo.rs"]
mod tests;