Fix #1200: [Rule] Unit ComparativeContainment to DecisionQUBO

docs/paper/reductions.typ

@@ -12343,6 +12343,156 @@ _Solution extraction and bounds._ Copy the subset bits; aggregate recovery is id
  The canonical instance and target coefficients are shown below. Its source witness is #repr(ex.solutions.at(0).source_config); the displayed target witness decodes to it and meets the objective or feasibility certificate described above.
 ]
 
+// ComparativeContainment → DecisionQUBO
+#let cc-qubo = load-example("ComparativeContainment", "DecisionQUBO",
+  source-variant: (weight: "One"), target-variant: (weight: "i64"))
+#let cc-source = cc-qubo.source.instance
+#let cc-target = cc-qubo.target.instance
+#let cc-sol = cc-qubo.solutions.at(0)
+#let cc-n = cc-source.universe_size
+#let cc-r = cc-source.r_sets.len()
+#let cc-s = cc-source.s_sets.len()
+#let cc-m = cc-r + 1
+#let cc-outside(a) = range(cc-n).filter(v => not a.contains(v))
+#let cc-g = (cc-source.r_sets + cc-source.s_sets).map(a => calc.max(0, cc-outside(a).len()-1)).sum(default: 0)
+#let cc-selected(bits) = range(cc-n).filter(v => bits.at(v))
+#let cc-count(family, bits) = family.filter(a => cc-selected(bits).all(v => a.contains(v))).len()
+#let cc-energy(bits) = cc-target.inner.entries.map(e => if bits.at(e.at(0)) and bits.at(e.at(1)) { e.at(2) } else { 0 }).sum(default: 0)
+#let cc-extend(bits) = {
+  let aux = ()
+  for a in cc-source.r_sets + cc-source.s_sets {
+    let outside = cc-outside(a)
+    if outside.len() >= 2 {
+      for k in range(2, outside.len()+1) {
+        aux.push(outside.slice(0,k).all(v => not bits.at(v)))
+      }
+    }
+  }
+  bits + aux
+}
+#assert(cc-target.inner.num_vars == cc-n + cc-g)
+#assert(cc-extend(cc-sol.source_config) == cc-sol.target_config)
+#assert(cc-source.s_sets.all(a => range(cc-n).all(v => a.contains(v))))
+#assert(cc-source.r_sets.all(a => cc-outside(a).len() == 2))
+#assert(cc-target.bound == -(cc-m * cc-g + cc-s))
+#assert(cc-count(cc-source.r_sets,cc-sol.source_config) >= cc-count(cc-source.s_sets,cc-sol.source_config))
+#assert(cc-energy(cc-sol.target_config) <= cc-target.bound)
+#let cc-faulty = cc-sol.target_config.enumerate().map(pair => if pair.at(0) == cc-n { not pair.at(1) } else { pair.at(1) })
+#let cc-candidate = cc-sol.source_config.enumerate().map(pair => pair.at(1) or pair.at(0) == cc-source.r_sets.at(0).at(0))
+#assert(cc-energy(cc-faulty) > cc-target.bound)
+#assert(cc-count(cc-source.r_sets,cc-candidate) < cc-count(cc-source.s_sets,cc-candidate))
+#assert(cc-energy(cc-extend(cc-candidate)) > cc-target.bound)
+#let cc-no-s = cc-source.s_sets + cc-source.s_sets
+#let cc-no-bound = cc-target.bound - cc-source.s_sets.len()
+#assert(cc-no-s.len() > cc-r)
+#assert(cc-energy(cc-sol.target_config) > cc-no-bound)
+#reduction-rule("ComparativeContainment", "DecisionQUBO",
+  example: true,
+  example-source-variant: (weight: "One"),
+  example-target-variant: (weight: "i64"),
+  extra: [
+    #pred-commands(
+      "pred create --example " + rule-spec(cc-qubo) + " -o containment.json",
+      "pred path ComparativeContainment/One DecisionQUBO/i64 -o paths.json",
+      "jq '.paths[] | select(.steps == 1)' paths.json > route.json",
+      "pred reduce containment.json --via route.json -o bundle.json",
+      "pred solve bundle.json --solver brute-force",
+      "pred evaluate containment.json --config " + cli-config(cc-sol.source_config),
+    )
+  ],
+)[
+  Unit Comparative Containment reduces to integer decision QUBO by chaining
+  containment indicators. The indicator viewpoint follows the
+  #link("https://github.com/CodingThrust/problem-reductions/issues/401")[antecedent containment formulation];
+  the quadratic product penalty is Rosenberg's, as presented in @borosgruber2014quadratization.
+][
+  _Construction._ Let $U={0,dots,n-1}$ be the universe, and $R,S$ the indexed
+  families of sizes $r,s$. Put the universe selection bits $x_v$ first. For
+  each entry $A$, its containment indicator is $t_A=product_(v in U without A)(1-x_v)$.
+  Enumerate its complement increasingly. An empty complement uses constant 1;
+  one element uses its complemented bit; otherwise use a fresh chain of
+  $|U without A|-1$ product bits. Process R before S in input order and never
+  share chain bits across entries. For each binary gate $z=a b$, add
+  $P(a,b,z)=a b-2a z-2b z+3z$. With $M=r+1$, form
+  $H=M sum_"gates" P + sum_(A in S)t_A-sum_(A in R)t_A=C+Q$.
+  Expand into unique upper-triangular entries (diagonals are linear, each
+  coefficient used once), discard zero coefficients, and use bound $-C$.
+
+  _Correctness._ ($arrow.r.double$) A feasible source subset extends by its
+  true prefix products. Every penalty vanishes and $H$ is its S count minus
+  its R count, at most zero; hence $Q<=-C$. ($arrow.l.double$) On any target
+  assignment, each indicator remains binary, so the signed indicator sum is
+  at least $-r$. The gate penalty is nonnegative, vanishes exactly at $z=a b$,
+  and is at least 1 otherwise. One faulty gate would give $H>=M-r=1$.
+  A target witness therefore has every gate correct. Induction along each
+  chain recovers actual containment, and $H<=0$ gives the source inequality.
+
+  _Solution extraction._ Validate the target predicate once and return its
+  first $n$ bits. Reject wrong lengths and assignments above the bound.
+  Completed `Or` values map identically, including NO without a witness.
+  An empty family entry contains only the empty selected subset; a full-universe
+  entry always contributes 1 and needs no gates. With empty universe the
+  single empty witness works exactly when $r>=s$. Empty families and repeated
+  entries retain their count semantics; repeated elements inside a set do
+  not change membership. In particular, empty R with one full-universe S entry
+  is NO: $Q=0$ and bound $-1$.
+
+  _Size and arithmetic._ Write $g=sum_A max(0,n-|A|-1)$ over all indexed
+  entries. There are exactly $n+g$ variables. A first gate has six nonconstant
+  entries and three off-diagonal terms; a later gate has five and three.
+  Each indicator adds at most one diagonal. Thus there are at most $6g+r+s$
+  stored entries and $3g$ off-diagonal terms; coefficient collection can
+  cancel entries. Since $g<=n(r+s)$, the registered bounds are $n+n(r+s)$
+  variables and $6n(r+s)+r+s$ off-diagonal terms. The latter is conservative;
+  neither registered formula asserts equality. Complement scans cost
+  $O(n(r+s))$ plus reading the raw listings; sparse collection adds a
+  logarithmic factor in the number of gate and indicator records.
+
+  Checked size arithmetic and conversions precede expansion. Require
+  $M=r+1<=2^63-1$ and $B=8M g+r+s<=2^63-1$, computing $M g$ before multiplying
+  by 8 so zero gates remain valid. Expanded gates have nonconstant absolute
+  coefficient sums at most $8M$ and constants at most $M$. Consequently the
+  nonconstant norm is at most $B$, and separately $|C|<=M g+r+s<=B$.
+  Every coefficient, accumulation prefix and evaluation on every assignment
+  fits `i64`, even for infeasible assignments. All operations remain checked;
+  an unrepresentable target returns a typed error.
+
+  *YES example.* The loaded source has universe size #cc-n, R entries
+  #repr(cc-source.r_sets) and S entries #repr(cc-source.s_sets). Here
+  $r=#cc-r$, $s=#cc-s$, $M=#cc-m$, $g=#cc-g$: the #cc-n universe bits and
+  #cc-g chain bits give #cc-target.inner.num_vars variables. The target has
+  #cc-target.inner.entries.len() stored entries, of which
+  #(cc-target.inner.entries.filter(e => e.at(0) != e.at(1)).len()) are
+  off-diagonal. Its constant is #(-cc-target.bound), bound #cc-target.bound,
+  coefficient norm #(cc-target.inner.entries.map(e => calc.abs(e.at(2))).sum(default: 0)),
+  and arithmetic certificate #(8*cc-m*cc-g+cc-r+cc-s).
+  Its R complement pairs are #repr(cc-source.r_sets.map(cc-outside)).
+  The full polynomial is specified by triples #repr(cc-target.inner.entries).
+  The selected subset #repr(cc-selected(cc-sol.source_config)) has R count
+  #cc-count(cc-source.r_sets,cc-sol.source_config) and S count
+  #cc-count(cc-source.s_sets,cc-sol.source_config). Its extension
+  #repr(cc-sol.target_config) has $Q=#cc-energy(cc-sol.target_config)$ and
+  $H=#(cc-energy(cc-sol.target_config)-cc-target.bound)$, and extracts back
+  to #repr(cc-sol.source_config). Flipping its first chain bit gives
+  $Q=#cc-energy(cc-faulty)$ and $H=#(cc-energy(cc-faulty)-cc-target.bound)$,
+  which fails the bound. The subset #repr(cc-selected(cc-candidate)) has
+  counts #cc-count(cc-source.r_sets,cc-candidate) and
+  #cc-count(cc-source.s_sets,cc-candidate); its true extension has
+  $Q=#cc-energy(cc-extend(cc-candidate))$, also rejected.
+
+  *NO example.* Double the loaded S family, keeping the same universe and R.
+  Its #cc-no-s.len() full-universe entries always contribute that count,
+  while at most #cc-r R entries can contain any selected subset. The count
+  difference is therefore at least #(cc-no-s.len()-cc-r), strictly positive.
+  There are still #cc-g gates and #cc-target.inner.num_vars variables, and
+  the same nonconstant Q. The constant becomes #(-cc-no-bound), bound
+  #cc-no-bound, and arithmetic certificate #(8*cc-m*cc-g+cc-r+cc-no-s.len()).
+  Every consistent target misses the bound by at least
+  #(cc-no-s.len()-cc-r); every inconsistent target has $H>=M-r=1$.
+  The previous vector still gives $Q=#cc-energy(cc-sol.target_config)$,
+  missing the new bound by #(cc-energy(cc-sol.target_config)-cc-no-bound).
+]
+
 #reduction-rule("ThreeDimensionalMatching", "DecisionQUBO", example: true,
  extra: [#direct-qubo-commands("ThreeDimensionalMatching", "DecisionQUBO")],
 )[

docs/paper/references.bib

@@ -2230,6 +2230,7 @@ @misc{borosgruber2014quadratization
   eprint = {1404.6538},
   archivePrefix = {arXiv},
   primaryClass = {math.OC},
+  url = {https://arxiv.org/abs/1404.6538},
   note = {Paper presented at ISAIM 2012; arXiv version 2014},
 }
 

src/rules/comparativecontainment_qubo.rs

@@ -0,0 +1,169 @@
+//! Unit containment indicators encoded by penalized binary product chains.
+use crate::models::{ComparativeContainment, Decision, QUBO};
+use crate::rules::qubo_helpers::QuboBuilder;
+use crate::rules::{ReduceTo, ReductionError, ReductionResult};
+use crate::types::One;
+use std::collections::HashSet;
+
+type Builder = QuboBuilder<ComparativeContainment<One>, Decision<QUBO<i64>>>;
+
+#[derive(Debug, Clone)]
+pub struct ReductionComparativeContainmentToQUBO {
+    target: Decision<QUBO<i64>>,
+    universe_size: usize,
+}
+
+impl ReductionResult for ReductionComparativeContainmentToQUBO {
+    type Source = ComparativeContainment<One>;
+    type Target = Decision<QUBO<i64>>;
+
+    fn target_problem(&self) -> &Self::Target {
+        &self.target
+    }
+
+    fn extract_solution(&self, bits: &Vec<bool>) -> crate::rules::ExtractionResult<Vec<bool>> {
+        crate::rules::traits::validate_target_witness(
+            &self.target,
+            bits,
+            |value| value.0,
+            "QUBO does not meet the containment threshold",
+        )?;
+        Ok(bits[..self.universe_size].to_vec())
+    }
+}
+
+#[crate::aggregate_reduction(identity)]
+impl crate::rules::AggregateReductionResult for ReductionComparativeContainmentToQUBO {}
+
+#[crate::reduction(transform = upper_bound {
+    num_vars = "universe_size + universe_size * (num_r_sets + num_s_sets)",
+    num_quadratic_terms = "6 * universe_size * (num_r_sets + num_s_sets) + num_r_sets + num_s_sets",
+})]
+impl ReduceTo<Decision<QUBO<i64>>> for ComparativeContainment<One> {
+    type Result = ReductionComparativeContainmentToQUBO;
+
+    fn reduce_to(&self) -> Result<Self::Result, ReductionError> {
+        construct(self, i64::MAX)
+    }
+}
+
+fn checked_sizes(
+    n: usize,
+    r: usize,
+    s: usize,
+    g: usize,
+    limit: i64,
+) -> Result<(usize, i64), ReductionError> {
+    let count = r.checked_add(s).ok_or_else(Builder::overflow)?;
+    let num_vars = n.checked_add(g).ok_or_else(Builder::overflow)?;
+    let r = i64::try_from(r).map_err(|_| Builder::overflow())?;
+    let count = i64::try_from(count).map_err(|_| Builder::overflow())?;
+    let g = i64::try_from(g).map_err(|_| Builder::overflow())?;
+    let penalty = r.checked_add(1).ok_or_else(Builder::overflow)?;
+    // Multiply by the gate count first: a gate-free target needs no 8*M.
+    let certificate = penalty
+        .checked_mul(g)
+        .and_then(|v| v.checked_mul(8))
+        .and_then(|v| v.checked_add(count))
+        .ok_or_else(Builder::overflow)?;
+    if penalty > limit || certificate > limit {
+        return Err(Builder::overflow());
+    }
+    Ok((num_vars, penalty))
+}
+
+fn construct(
+    source: &ComparativeContainment<One>,
+    limit: i64,
+) -> Result<ReductionComparativeContainmentToQUBO, ReductionError> {
+    let n = source.universe_size();
+    let mut g = 0usize;
+    for entry in source.r_sets().iter().chain(source.s_sets()) {
+        let distinct: HashSet<_> = entry.iter().collect();
+        let k = n - distinct.len();
+        g = g
+            .checked_add(k.saturating_sub(1))
+            .ok_or_else(Builder::overflow)?;
+    }
+    let (num_vars, penalty) = checked_sizes(n, source.num_r_sets(), source.num_s_sets(), g, limit)?;
+    let p = i128::from(penalty);
+    let mut q = Builder::new(num_vars);
+    let mut next = n;
+    for (family, sign) in [(source.r_sets(), -1i128), (source.s_sets(), 1)] {
+        for entry in family {
+            let present: HashSet<_> = entry.iter().copied().collect();
+            let mut complement = (0..n).filter(|v| !present.contains(v));
+            let Some(u) = complement.next() else {
+                q.add_constant(sign)?;
+                continue;
+            };
+            let Some(v) = complement.next() else {
+                q.add_constant(sign)?;
+                q.add(u, u, -sign)?;
+                continue;
+            };
+            let mut z = next;
+            next = next.checked_add(1).ok_or_else(Builder::overflow)?;
+            // P(1-x_u,1-x_v,z) = 1-x_u-x_v+x_u*x_v-z+2*x_u*z+2*x_v*z.
+            q.add_constant(p)?;
+            for (a, b, c) in [
+                (u, u, -p),
+                (v, v, -p),
+                (u, v, p),
+                (z, z, -p),
+                (u, z, 2 * p),
+                (v, z, 2 * p),
+            ] {
+                q.add(a, b, c)?;
+            }
+            for v in complement {
+                let previous = z;
+                z = next;
+                next = next.checked_add(1).ok_or_else(Builder::overflow)?;
+                // P(w,1-x_v,z) = w-w*x_v-2*w*z+z+2*x_v*z.
+                for (a, b, c) in [
+                    (previous, previous, p),
+                    (previous, v, -p),
+                    (previous, z, -2 * p),
+                    (z, z, p),
+                    (v, z, 2 * p),
+                ] {
+                    q.add(a, b, c)?;
+                }
+            }
+            q.add(z, z, sign)?;
+        }
+    }
+    let (target, constant) = q.finish()?;
+    let bound = constant.checked_neg().ok_or_else(Builder::overflow)?;
+    Ok(ReductionComparativeContainmentToQUBO {
+        target: Decision::new(target, bound),
+        universe_size: n,
+    })
+}
+
+#[cfg(feature = "example-db")]
+pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
+    vec![crate::example_db::specs::RuleExampleSpec {
+        id: "comparativecontainment_qubo",
+        build: || {
+            let source = ComparativeContainment::<One>::new(
+                4,
+                vec![vec![0, 1], vec![1, 2], vec![1, 3]],
+                vec![vec![0, 1, 2, 3]; 2],
+            )
+            .expect("canonical containment instance must be valid");
+            crate::example_db::specs::rule_example_with_witness::<_, Decision<QUBO<i64>>>(
+                source,
+                crate::export::SolutionPair {
+                    source_config: serde_json::json!([false, true, false, false]),
+                    target_config: serde_json::json!([false, true, false, false, true, true, true]),
+                },
+            )
+        },
+    }]
+}
+
+#[cfg(test)]
+#[path = "../unit_tests/rules/comparativecontainment_qubo.rs"]
+mod tests;

src/rules/mod.rs

@@ -45,6 +45,7 @@ pub(crate) mod twodimensionalconsecutivesets_ilp;
 pub(crate) use graph::{RecoveryStep, recover_completed_result};
 mod binpacking_casts;
 mod comparativecontainment_casts;
+pub(crate) mod comparativecontainment_qubo;
 pub(crate) mod graph_helpers;
 pub(crate) mod graphpartitioning_maxcut;
 pub(crate) mod graphpartitioning_qubo;
@@ -498,6 +499,7 @@ pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::Ru
         optimallineararrangement_sequencingtominimizeweightedcompletiontime::canonical_rule_example_specs(),
     );
     specs.extend(setsplitting_betweenness::canonical_rule_example_specs());
+    specs.extend(comparativecontainment_qubo::canonical_rule_example_specs());
     specs.extend(setsplitting_naesatisfiability::canonical_rule_example_specs());
     specs.extend(satisfiability_integralflowhomologousarcs::canonical_rule_example_specs());
     specs.extend(satisfiability_naesatisfiability::canonical_rule_example_specs());