Add ClosestSubstring -> ILP reduction (#1035)

Position-character ILP encoding for ClosestSubstring: binary variables
x_{r,a} for the center substring's character at each position, window
choice indicators y_{i,p} (exactly one window per source string), and
an integer radius R bounded by per-window Hamming-distance constraints
plus a tight R <= ell upper-bound constraint critical for ILP solver
performance.

Closes #1035.

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

Author: isPANN

docs/paper/reductions.typ

@@ -13893,6 +13893,31 @@ The following reductions to Integer Linear Programming are straightforward formu
   _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("ClosestSubstring", "ILP")[
+  Integer variables select one alphabet symbol at each center position and one window start per input string. A conditional radius constraint is activated by the window-choice indicator and upper-bounds the Hamming distance between the center and the selected window of each string.
+][
+  _Construction._ Given alphabet $Sigma$ of size $q$, $n$ input strings $s_1, dots, s_n$ over $Sigma$, and window length $ell$ with $W_i = |s_i| - ell + 1$:
+
+  _Variables:_ (1) $x_(r, a) in {0, 1}$ for $r in {0, dots, ell - 1}$ and $a in {0, dots, q - 1}$: $x_(r, a) = 1$ iff the center has symbol $a$ at position $r$. (2) $y_(i, p) in {0, 1}$ for input string $s_i$ and window start $p in {0, dots, W_i - 1}$: $y_(i, p) = 1$ iff window $p$ is selected from $s_i$. (3) Nonnegative integer $R$: an upper bound on the worst-case Hamming distance.
+
+  _Constraints:_ (1) Center assignment: $sum_(a = 0)^(q - 1) x_(r, a) = 1$ for every position $r$. (2) Window choice: $sum_(p = 0)^(W_i - 1) y_(i, p) = 1$ for every input string $s_i$. (3) Conditional radius: $R + sum_(r = 0)^(ell - 1) x_(r, s_i [p + r]) - ell dot y_(i, p) >= 0$ for every $(i, p)$. When $y_(i, p) = 1$, this is equivalent to $R >= ell - sum_r x_(r, s_i [p + r]) = d_H (c, s_i [p .. p + ell))$; when $y_(i, p) = 0$, the constraint reduces to $R + (text("nonneg match count")) >= 0$, which holds automatically.
+
+  _Objective:_ Minimize $R$.
+
+  The ILP is:
+  $
+    "minimize" quad & R \
+    "subject to" quad & sum_(a = 0)^(q - 1) x_(r, a) = 1 quad forall r in {0, dots, ell - 1} \
+    & sum_(p = 0)^(W_i - 1) y_(i, p) = 1 quad forall i in {1, dots, n} \
+    & R + sum_(r = 0)^(ell - 1) x_(r, s_i [p + r]) - ell dot y_(i, p) >= 0 quad forall i, p \
+    & x_(r, a), y_(i, p) in {0, 1}, quad R in ZZ_(>= 0).
+  $
+
+  _Correctness._ ($arrow.r.double$) Given an optimal center $c^* in Sigma^ell$ and optimal window starts $p_1^*, dots, p_n^*$, set $x_(r, c^*[r]) = 1$, $y_(i, p_i^*) = 1$, and $R = max_i d_H (c^*, s_i [p_i^* .. p_i^* + ell))$. The assignment and window-choice constraints hold by construction. For each pair $(i, p_i^*)$ the radius constraint becomes $R >= d_H (c^*, s_i [p_i^* .. p_i^* + ell))$, which holds with equality at the worst case; for every other $(i, p)$ with $y_(i, p) = 0$ the constraint is redundant. ($arrow.l.double$) The assignment and window-choice constraints force each block of $x$ and each block of $y$ to be one-hot, encoding a center $c$ and one window per input string. The conditional radius constraint is active exactly on the selected windows and forces $R >= d_H (c, s_i [p_i .. p_i + ell))$ for every $i$, so $R$ is at least the worst-case selected Hamming distance. Minimizing $R$ therefore minimizes the maximum Hamming distance over chosen windows.
+
+  _Solution extraction._ For each position $r$, read the unique symbol $a$ with $x_(r, a) = 1$ as the center symbol; for each input string $s_i$, read the unique $p$ with $y_(i, p) = 1$ as the selected window start.
+]
+
 #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/closestsubstring_ilp.rs

@@ -0,0 +1,225 @@
+//! Reduction from ClosestSubstring to ILP (Integer Linear Programming).
+//!
+//! Given an alphabet of size `q`, `n` input strings `s_1, ..., s_n` (not
+//! necessarily of equal length), and a window length `ell`, the goal is to
+//! pick a center `c in Sigma^ell` and one length-`ell` window from each input
+//! string that together minimize the worst-case Hamming distance between the
+//! center and any chosen window. The ILP encoding combines the
+//! center-selection variables of ClosestString with one-hot window-choice
+//! indicators, plus a radius variable that is active only on each selected
+//! window.
+//!
+//! - Integer `x_{r, a}` for `r in {0, ..., ell - 1}` and
+//!   `a in {0, ..., q - 1}`: `x_{r, a} = 1` iff the center has symbol `a` at
+//!   position `r`. The non-negativity of ILP variables together with the
+//!   assignment constraint forces every `x_{r, a} in {0, 1}`.
+//! - Integer `y_{i, p}` for input string `s_i` and window start
+//!   `p in {0, ..., W_i - 1}` where `W_i = |s_i| - ell + 1`: `y_{i, p} = 1` iff
+//!   window `p` is selected from string `s_i`.
+//! - Nonnegative integer radius variable `R`.
+//! - Assignment constraint: `sum_a x_{r, a} = 1` for every position `r`.
+//! - Window-choice constraint: `sum_p y_{i, p} = 1` for every input string.
+//! - Conditional radius constraint per `(i, p)`:
+//!   `R + sum_{r} x_{r, s_i[p + r]} - ell * y_{i, p} >= 0`.
+//!   When `y_{i, p} = 1`, this becomes `R >= d_H(c, s_i[p..p + ell))`; when
+//!   `y_{i, p} = 0`, the constraint is automatically satisfied.
+//! - 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::ClosestSubstring;
+use crate::reduction;
+use crate::rules::traits::{ReduceTo, ReductionResult};
+
+/// Result of reducing ClosestSubstring to ILP.
+///
+/// Variable layout (`ILP<i32>`, all non-negative):
+/// - `x_{r, a}` at index `r * alphabet_size + a` for `r in [0, ell)` and
+///   `a in [0, q)`, forced into `{0, 1}` by the assignment constraints.
+/// - `y_{i, p}` at index `q * ell + window_offsets[i] + p` for input string
+///   `s_i` and window start `p in [0, W_i)`, forced into `{0, 1}` by the
+///   window-choice constraints.
+/// - `R` (radius) at index `q * ell + total_num_windows`, a non-negative
+///   integer in `[0, ell]`.
+#[derive(Debug, Clone)]
+pub struct ReductionClosestSubstringToILP {
+    target: ILP<i32>,
+    alphabet_size: usize,
+    substring_length: usize,
+    /// Prefix sums of per-string window counts: `window_offsets[i]` is the
+    /// number of `y_{j, p}` variables for `j < i`. Has length `num_strings`.
+    window_offsets: Vec<usize>,
+    /// `window_counts[i] = W_i = |s_i| - ell + 1`.
+    window_counts: Vec<usize>,
+}
+
+impl ReductionResult for ReductionClosestSubstringToILP {
+    type Source = ClosestSubstring;
+    type Target = ILP<i32>;
+
+    fn target_problem(&self) -> &ILP<i32> {
+        &self.target
+    }
+
+    /// Decode the integer ILP assignment into the source config layout.
+    ///
+    /// `ClosestSubstring::evaluate` expects `config` of length `ell + n`: the
+    /// first `ell` entries are the center symbols, the remaining `n` entries
+    /// are per-string window starts. For each center position `r`, we pick the
+    /// unique alphabet symbol `a` with `x_{r, a} = 1`; for each input string
+    /// `s_i`, we pick the unique window start `p` with `y_{i, p} = 1`. When no
+    /// indicator is set to 1 in some block (which only happens on partial /
+    /// infeasible ILP solutions), we fall back to 0 so the returned vector
+    /// still has the expected shape.
+    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
+        let q = self.alphabet_size;
+        let ell = self.substring_length;
+        let y_base = q * ell;
+
+        let mut out = Vec::with_capacity(ell + self.window_counts.len());
+
+        // Center symbols.
+        for r in 0..ell {
+            let symbol = (0..q)
+                .find(|&a| target_solution.get(r * q + a).copied().unwrap_or(0) == 1)
+                .unwrap_or(0);
+            out.push(symbol);
+        }
+
+        // Window starts.
+        for (i, &w_i) in self.window_counts.iter().enumerate() {
+            let start = (0..w_i)
+                .find(|&p| {
+                    target_solution
+                        .get(y_base + self.window_offsets[i] + p)
+                        .copied()
+                        .unwrap_or(0)
+                        == 1
+                })
+                .unwrap_or(0);
+            out.push(start);
+        }
+
+        out
+    }
+}
+
+#[reduction(
+    overhead = {
+        num_vars = "alphabet_size * substring_length + total_num_windows + 1",
+        num_constraints = "substring_length + num_strings + total_num_windows + 1",
+    }
+)]
+impl ReduceTo<ILP<i32>> for ClosestSubstring {
+    type Result = ReductionClosestSubstringToILP;
+
+    fn reduce_to(&self) -> Self::Result {
+        let q = self.alphabet_size();
+        let ell = self.substring_length();
+        let strings = self.strings();
+        let n = strings.len();
+
+        let window_counts: Vec<usize> = strings.iter().map(|s| s.len() - ell + 1).collect();
+        let mut window_offsets: Vec<usize> = Vec::with_capacity(n);
+        {
+            let mut acc = 0usize;
+            for &w in &window_counts {
+                window_offsets.push(acc);
+                acc += w;
+            }
+        }
+        let total_windows: usize = window_counts.iter().sum();
+
+        let x_idx = |r: usize, a: usize| -> usize { r * q + a };
+        let y_base = q * ell;
+        let y_idx = |i: usize, p: usize| -> usize { y_base + window_offsets[i] + p };
+        let r_idx = y_base + total_windows;
+        let num_vars = r_idx + 1;
+
+        let mut constraints: Vec<LinearConstraint> =
+            Vec::with_capacity(ell + n + total_windows + 1);
+
+        // Assignment constraints: exactly one symbol per center position.
+        // Together with the non-negativity built into ILP<i32>, this also
+        // forces every x_{r, a} to lie in {0, 1}.
+        for r in 0..ell {
+            let terms: Vec<(usize, f64)> = (0..q).map(|a| (x_idx(r, a), 1.0)).collect();
+            constraints.push(LinearConstraint::eq(terms, 1.0));
+        }
+
+        // Tight upper bound on R: the worst-case Hamming distance over a
+        // length-ell window is at most ell. Added as a single-term `<=`
+        // constraint so the solver's bound-tightening pass (which scans for
+        // exactly this pattern) picks it up. Without this, R defaults to the
+        // full i32 domain, which severely degrades HiGHS performance even on
+        // tiny instances.
+        constraints.push(LinearConstraint::le(vec![(r_idx, 1.0)], ell as f64));
+
+        // Window-choice constraints: exactly one window per input string.
+        // Combined with non-negativity, this forces every y_{i, p} in {0, 1}.
+        for (i, &w_i) in window_counts.iter().enumerate() {
+            let terms: Vec<(usize, f64)> = (0..w_i).map(|p| (y_idx(i, p), 1.0)).collect();
+            constraints.push(LinearConstraint::eq(terms, 1.0));
+        }
+
+        // Conditional radius constraints: for every (input string, window
+        // start) pair, R + sum_r x_{r, s_i[p + r]} - ell * y_{i, p} >= 0.
+        // - If y_{i, p} = 1: R >= ell - sum_r x_{r, s_i[p + r]} = d_H(c, window).
+        // - If y_{i, p} = 0: the LHS is R + (nonneg match count) >= 0,
+        //   automatically satisfied because R >= 0.
+        for (i, s) in strings.iter().enumerate() {
+            for p in 0..window_counts[i] {
+                let mut terms: Vec<(usize, f64)> = Vec::with_capacity(ell + 2);
+                terms.push((r_idx, 1.0));
+                for r in 0..ell {
+                    terms.push((x_idx(r, s[p + r]), 1.0));
+                }
+                terms.push((y_idx(i, p), -(ell as f64)));
+                constraints.push(LinearConstraint::ge(terms, 0.0));
+            }
+        }
+
+        // Objective: minimize R.
+        let objective = vec![(r_idx, 1.0)];
+
+        let target = ILP::new(num_vars, constraints, objective, ObjectiveSense::Minimize);
+
+        ReductionClosestSubstringToILP {
+            target,
+            alphabet_size: q,
+            substring_length: ell,
+            window_offsets,
+            window_counts,
+        }
+    }
+}
+
+#[cfg(feature = "example-db")]
+pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
+    vec![crate::example_db::specs::RuleExampleSpec {
+        id: "closestsubstring_to_ilp",
+        build: || {
+            // Canonical issue #1033 instance: binary alphabet, length-3
+            // windows on three length-5 strings. Optimum radius is 1; one
+            // optimal center is 010 with windows (0, 1, 0) selecting 000,
+            // 010, 110 from s_1, s_2, s_3 respectively.
+            let source = ClosestSubstring::new(
+                2,
+                vec![
+                    vec![0, 0, 0, 1, 1],
+                    vec![1, 0, 1, 0, 0],
+                    vec![1, 1, 0, 0, 1],
+                ],
+                3,
+            );
+            crate::example_db::specs::rule_example_via_ilp::<_, i32>(source)
+        },
+    }]
+}
+
+#[cfg(test)]
+#[path = "../unit_tests/rules/closestsubstring_ilp.rs"]
+mod tests;

src/rules/mod.rs

@@ -166,6 +166,8 @@ pub(crate) mod circuit_ilp;
 #[cfg(feature = "ilp-solver")]
 pub(crate) mod closeststring_ilp;
 #[cfg(feature = "ilp-solver")]
+pub(crate) mod closestsubstring_ilp;
+#[cfg(feature = "ilp-solver")]
 pub(crate) mod clustering_ilp;
 #[cfg(feature = "ilp-solver")]
 pub(crate) mod coloring_ilp;
@@ -551,6 +553,7 @@ pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::Ru
         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(closestsubstring_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());