Add ClosestString -> ILP reduction (#1034)

Position-character ILP encoding for ClosestString: binary variables
x_{j,a} encode the consensus string character at position j (with
one-hot assignment constraints), and a single integer radius variable
R is bounded by per-string Hamming-distance constraints. The objective
minimizes R.

Closes #1034.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

Author: isPANN

docs/paper/reductions.typ

@@ -13869,6 +13869,30 @@ The following reductions to Integer Linear Programming are straightforward formu
   _Solution extraction._ Sort the selected vertices by position in $s_1$. Read off the characters to obtain the common subsequence, then pad to `max_length` with the padding symbol.
 ]
 
+#reduction-rule("ClosestString", "ILP")[
+  Binary variables select one alphabet symbol at each center position. An auxiliary radius variable upper-bounds the Hamming distance from the chosen center to every input string and is minimized.
+][
+  _Construction._ Given alphabet $Sigma$ of size $q$, $n$ input strings $s_1, dots, s_n in Sigma^m$ of common length $m$:
+
+  _Variables:_ (1) $x_(j, a) in {0, 1}$ for $j in {0, dots, m - 1}$ and $a in {0, dots, q - 1}$: $x_(j, a) = 1$ iff the center has symbol $a$ at position $j$. (2) Nonnegative integer $R$: an upper bound on the worst-case Hamming distance.
+
+  _Constraints:_ (1) Assignment: $sum_(a = 0)^(q - 1) x_(j, a) = 1$ for every position $j$. Combined with the nonnegativity built into the ILP, this also forces every $x_(j, a) in {0, 1}$. (2) Radius: $R + sum_(j = 0)^(m - 1) x_(j, s_i [j]) >= m$ for every input string $s_i$, which is equivalent to $R >= m - sum_j x_(j, s_i [j]) = d_H (c, s_i)$.
+
+  _Objective:_ Minimize $R$.
+
+  The ILP is:
+  $
+    "minimize" quad & R \
+    "subject to" quad & sum_(a = 0)^(q - 1) x_(j, a) = 1 quad forall j in {0, dots, m - 1} \
+    & R + sum_(j = 0)^(m - 1) x_(j, s_i [j]) >= m quad forall i in {1, dots, n} \
+    & x_(j, a) in {0, 1}, quad R in ZZ_(>= 0).
+  $
+
+  _Correctness._ ($arrow.r.double$) Given an optimal center $c^* in Sigma^m$, set $x_(j, c^*[j]) = 1$ for every $j$ and let $R = max_i d_H (c^*, s_i)$. Each assignment constraint is satisfied, and each radius constraint reduces to $R >= d_H (c^*, s_i)$, which holds with equality at the worst case. ($arrow.l.double$) The assignment constraints force each $x_(j, *)$ block to be a one-hot vector, hence encode a unique center string $c$. The radius constraint then gives $R >= d_H (c, s_i)$ for every $i$, so $R >= max_i d_H (c, s_i)$. Minimizing $R$ therefore minimizes the worst-case Hamming distance.
+
+  _Solution extraction._ For each position $j$, read the unique symbol $a$ with $x_(j, a) = 1$; the resulting length-$m$ vector is the source center.
+]
+
 #reduction-rule("LongestCommonSubsequence", "ILP")[
   An optimization ILP formulation maximizes the length of a common subsequence. Binary variables choose a symbol (or padding) at each witness position. Match variables link active positions to source string indices, and the objective maximizes the number of non-padding positions.
 ][

src/rules/closeststring_ilp.rs

@@ -0,0 +1,141 @@
+//! Reduction from ClosestString to ILP (Integer Linear Programming).
+//!
+//! Given an alphabet of size `q`, `n` input strings of common length `m`, and
+//! the goal of finding a center `c in Sigma^m` that minimizes the maximum
+//! Hamming distance to every input string, the natural encoding picks one
+//! alphabet symbol at every center position and constrains a radius variable
+//! to upper-bound every per-string Hamming distance:
+//!
+//! - Binary `x_{j, a} in {0, 1}` for `j in {0, ..., m - 1}` and `a in
+//!   {0, ..., q - 1}`: `x_{j, a} = 1` iff the chosen center has symbol `a` at
+//!   position `j`.
+//! - Nonnegative integer radius variable `R`.
+//! - Assignment constraint: `sum_a x_{j, a} = 1` for every position `j`.
+//!   Because every ILP variable is a nonnegative integer, this also forces
+//!   each `x_{j, a} in {0, 1}`.
+//! - Radius constraint per input string `s_i`:
+//!   `R + sum_j x_{j, s_i[j]} >= m`, which is equivalent to `R >= d_H(c, s_i)`.
+//! - Objective: minimize `R`.
+//!
+//! Reference: Ming Li, Bin Ma, and Lusheng Wang, "On the closest string and
+//! substring problems," Journal of the ACM 49(2):157-171, 2002.
+//! <https://doi.org/10.1145/506147.506150>
+
+use crate::models::algebraic::{LinearConstraint, ObjectiveSense, ILP};
+use crate::models::misc::ClosestString;
+use crate::reduction;
+use crate::rules::traits::{ReduceTo, ReductionResult};
+
+/// Result of reducing ClosestString to ILP.
+///
+/// Variable layout (`ILP<i32>`, all non-negative):
+/// - `x_{j, a}` at index `j * alphabet_size + a` for `j in [0, m)` and
+///   `a in [0, q)`, bounded to `{0, 1}`.
+/// - `R` (radius) at index `m * q`, an integer in `[0, m]`.
+#[derive(Debug, Clone)]
+pub struct ReductionClosestStringToILP {
+    target: ILP<i32>,
+    alphabet_size: usize,
+    string_length: usize,
+}
+
+impl ReductionResult for ReductionClosestStringToILP {
+    type Source = ClosestString;
+    type Target = ILP<i32>;
+
+    fn target_problem(&self) -> &ILP<i32> {
+        &self.target
+    }
+
+    /// Decode the integer ILP assignment into the source center config.
+    ///
+    /// For every position `j`, choose the unique alphabet symbol `a` with
+    /// `x_{j, a} = 1`. If the target assignment is missing or none of the
+    /// per-position `x_{j, *}` variables are set to 1, we fall back to symbol
+    /// `0` so the returned vector still has the expected length; partial /
+    /// infeasible ILP solutions are the caller's responsibility.
+    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
+        let q = self.alphabet_size;
+        (0..self.string_length)
+            .map(|j| {
+                (0..q)
+                    .find(|&a| target_solution.get(j * q + a).copied().unwrap_or(0) == 1)
+                    .unwrap_or(0)
+            })
+            .collect()
+    }
+}
+
+#[reduction(
+    overhead = {
+        num_vars = "alphabet_size * string_length + 1",
+        num_constraints = "string_length + num_strings",
+    }
+)]
+impl ReduceTo<ILP<i32>> for ClosestString {
+    type Result = ReductionClosestStringToILP;
+
+    fn reduce_to(&self) -> Self::Result {
+        let q = self.alphabet_size();
+        let m = self.string_length();
+        let strings = self.strings();
+        let n = strings.len();
+
+        let x_idx = |j: usize, a: usize| -> usize { j * q + a };
+        let r_idx = q * m;
+        let num_vars = q * m + 1;
+
+        let mut constraints: Vec<LinearConstraint> = Vec::with_capacity(m + n);
+
+        // Assignment constraints: exactly one symbol per center position.
+        // Together with the non-negativity built into ILP<i32>, this also
+        // forces every x_{j, a} to lie in {0, 1}.
+        for j in 0..m {
+            let terms: Vec<(usize, f64)> = (0..q).map(|a| (x_idx(j, a), 1.0)).collect();
+            constraints.push(LinearConstraint::eq(terms, 1.0));
+        }
+
+        // Radius constraints: R + sum_j x_{j, s_i[j]} >= m.
+        // Equivalently, R >= m - sum_j x_{j, s_i[j]} = d_H(c, s_i).
+        for s in strings.iter() {
+            let mut terms: Vec<(usize, f64)> = Vec::with_capacity(m + 1);
+            terms.push((r_idx, 1.0));
+            for (j, &symbol) in s.iter().enumerate() {
+                terms.push((x_idx(j, symbol), 1.0));
+            }
+            constraints.push(LinearConstraint::ge(terms, m as f64));
+        }
+
+        // Objective: minimize R.
+        let objective = vec![(r_idx, 1.0)];
+
+        let target = ILP::new(num_vars, constraints, objective, ObjectiveSense::Minimize);
+
+        ReductionClosestStringToILP {
+            target,
+            alphabet_size: q,
+            string_length: m,
+        }
+    }
+}
+
+#[cfg(feature = "example-db")]
+pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
+    vec![crate::example_db::specs::RuleExampleSpec {
+        id: "closeststring_to_ilp",
+        build: || {
+            // Canonical issue #1032 instance: binary alphabet, the four length-3
+            // strings 000, 011, 101, 110. Optimum radius is 2 (achieved by any
+            // binary length-3 center, e.g. 000).
+            let source = ClosestString::new(
+                2,
+                vec![vec![0, 0, 0], vec![0, 1, 1], vec![1, 0, 1], vec![1, 1, 0]],
+            );
+            crate::example_db::specs::rule_example_via_ilp::<_, i32>(source)
+        },
+    }]
+}
+
+#[cfg(test)]
+#[path = "../unit_tests/rules/closeststring_ilp.rs"]
+mod tests;

src/rules/mod.rs

@@ -164,6 +164,8 @@ pub(crate) mod capacityassignment_ilp;
 #[cfg(feature = "ilp-solver")]
 pub(crate) mod circuit_ilp;
 #[cfg(feature = "ilp-solver")]
+pub(crate) mod closeststring_ilp;
+#[cfg(feature = "ilp-solver")]
 pub(crate) mod clustering_ilp;
 #[cfg(feature = "ilp-solver")]
 pub(crate) mod coloring_ilp;
@@ -548,6 +550,7 @@ pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::Ru
         specs.extend(boundedcomponentspanningforest_ilp::canonical_rule_example_specs());
         specs.extend(capacityassignment_ilp::canonical_rule_example_specs());
         specs.extend(circuit_ilp::canonical_rule_example_specs());
+        specs.extend(closeststring_ilp::canonical_rule_example_specs());
         specs.extend(clustering_ilp::canonical_rule_example_specs());
         specs.extend(coloring_ilp::canonical_rule_example_specs());
         specs.extend(consecutiveblockminimization_ilp::canonical_rule_example_specs());

src/unit_tests/rules/closeststring_ilp.rs

@@ -0,0 +1,124 @@
+use super::*;
+use crate::models::algebraic::{ObjectiveSense, ILP};
+use crate::models::misc::ClosestString;
+use crate::rules::test_helpers::assert_bf_vs_ilp;
+use crate::solvers::{BruteForce, ILPSolver, Solver};
+use crate::traits::Problem;
+use crate::types::Min;
+
+/// Canonical issue #1032 instance: binary alphabet, four length-3 strings
+/// 000, 011, 101, 110. Every binary length-3 center attains radius exactly 2;
+/// no center attains radius 1.
+fn issue_instance() -> ClosestString {
+    ClosestString::new(
+        2,
+        vec![vec![0, 0, 0], vec![0, 1, 1], vec![1, 0, 1], vec![1, 1, 0]],
+    )
+}
+
+#[test]
+fn test_closeststring_to_ilp_structure() {
+    let source = issue_instance();
+    let reduction = ReduceTo::<ILP<i32>>::reduce_to(&source);
+    let ilp = reduction.target_problem();
+
+    // q = 2, m = 3 -> 2*3 + 1 = 7 variables.
+    assert_eq!(ilp.num_vars, 7);
+    // m = 3 assignment constraints + n = 4 radius constraints = 7.
+    assert_eq!(ilp.constraints.len(), 7);
+    assert_eq!(ilp.sense, ObjectiveSense::Minimize);
+
+    // The objective puts weight 1 on the radius variable only.
+    assert_eq!(ilp.objective.len(), 1);
+    let (r_idx, r_coeff) = ilp.objective[0];
+    assert_eq!(r_idx, 2 * 3);
+    assert!((r_coeff - 1.0).abs() < 1e-9);
+
+    // Each assignment constraint has q = 2 terms and rhs = 1.
+    for c in ilp.constraints.iter().take(3) {
+        assert_eq!(c.terms.len(), 2);
+        assert!((c.rhs - 1.0).abs() < 1e-9);
+    }
+
+    // Each radius constraint has m + 1 = 4 terms (one per position + R) and
+    // rhs = m = 3.
+    for c in ilp.constraints.iter().skip(3) {
+        assert_eq!(c.terms.len(), 4);
+        assert!((c.rhs - 3.0).abs() < 1e-9);
+    }
+}
+
+#[test]
+fn test_closeststring_to_ilp_closed_loop() {
+    let source = issue_instance();
+    let reduction = ReduceTo::<ILP<i32>>::reduce_to(&source);
+
+    let bf_value = BruteForce::new().solve(&source);
+    let ilp_solution = ILPSolver::new()
+        .solve(reduction.target_problem())
+        .expect("ILP should be solvable");
+    let extracted = reduction.extract_solution(&ilp_solution);
+    let extracted_value = source.evaluate(&extracted);
+
+    // The extracted center must be syntactically valid and match the BF optimum.
+    assert!(extracted_value.is_valid());
+    assert_eq!(extracted_value, bf_value);
+    // Sanity: the canonical instance has optimum radius 2.
+    assert_eq!(extracted_value, Min(Some(2)));
+}
+
+#[test]
+fn test_closeststring_to_ilp_bf_vs_ilp() {
+    let source = issue_instance();
+    let reduction = ReduceTo::<ILP<i32>>::reduce_to(&source);
+    assert_bf_vs_ilp(&source, &reduction);
+}
+
+#[test]
+fn test_closeststring_to_ilp_extract_known_center() {
+    // Build the binary encoding of the center 000 by hand:
+    // x_{0,0}=x_{1,0}=x_{2,0}=1, others 0, R = 2.
+    let source = issue_instance();
+    let reduction = ReduceTo::<ILP<i32>>::reduce_to(&source);
+
+    let mut target_solution = vec![0usize; reduction.target_problem().num_vars];
+    target_solution[0] = 1; // x_{0,0}
+    target_solution[2] = 1; // x_{1,0}
+    target_solution[4] = 1; // x_{2,0}
+    target_solution[6] = 2; // R = 2
+
+    let extracted = reduction.extract_solution(&target_solution);
+    assert_eq!(extracted, vec![0, 0, 0]);
+    assert_eq!(source.evaluate(&extracted), Min(Some(2)));
+}
+
+#[test]
+fn test_closeststring_to_ilp_ternary_alphabet() {
+    // q = 3, m = 2, three strings forcing a nonzero radius. The optimum
+    // radius is 1 (any center matches at least one position of every string).
+    let source = ClosestString::new(3, vec![vec![0, 1], vec![1, 2], vec![2, 0]]);
+    let reduction = ReduceTo::<ILP<i32>>::reduce_to(&source);
+    let ilp = reduction.target_problem();
+
+    // q * m + 1 = 3 * 2 + 1 = 7 variables; m + n = 2 + 3 = 5 constraints.
+    assert_eq!(ilp.num_vars, 7);
+    assert_eq!(ilp.constraints.len(), 5);
+
+    assert_bf_vs_ilp(&source, &reduction);
+}
+
+#[test]
+fn test_closeststring_to_ilp_single_string_zero_radius() {
+    // A single input string: the center equals the input and the optimum
+    // radius is 0. This guards against off-by-one errors in the radius
+    // constraints.
+    let source = ClosestString::new(2, vec![vec![1, 0, 1, 1]]);
+    let reduction = ReduceTo::<ILP<i32>>::reduce_to(&source);
+
+    let ilp_solution = ILPSolver::new()
+        .solve(reduction.target_problem())
+        .expect("ILP should be solvable");
+    let extracted = reduction.extract_solution(&ilp_solution);
+    assert_eq!(extracted, vec![1, 0, 1, 1]);
+    assert_eq!(source.evaluate(&extracted), Min(Some(0)));
+}