Fix #300: Add QuadraticAssignment model

docs/paper/reductions.typ

@@ -110,6 +110,7 @@
   "SubsetSum": [Subset Sum],
   "MinimumFeedbackArcSet": [Minimum Feedback Arc Set],
   "MinimumFeedbackVertexSet": [Minimum Feedback Vertex Set],
+  "QuadraticAssignment": [Quadratic Assignment],
   "SequencingWithReleaseTimesAndDeadlines": [Sequencing with Release Times and Deadlines],
   "MultipleChoiceBranching": [Multiple Choice Branching],
   "ShortestCommonSupersequence": [Shortest Common Supersequence],
@@ -1850,6 +1851,93 @@ NP-completeness was established by Garey, Johnson, and Stockmeyer @gareyJohnsonS
   ]
 }
 
+#{
+  let x = load-model-example("QuadraticAssignment")
+  let C = x.instance.cost_matrix
+  let D = x.instance.distance_matrix
+  let n = C.len()
+  let m = D.len()
+  let sol = (config: x.optimal_config, metric: x.optimal_value)
+  let fstar = sol.config
+  let cost-star = sol.metric.Valid
+  // Convert integer matrix to math.mat content
+  let to-mat(m) = math.mat(..m.map(row => row.map(v => $#v$)))
+  // Compute identity assignment cost
+  let id-cost = range(n).fold(0, (acc, i) =>
+    range(n).fold(acc, (acc2, j) =>
+      if i != j { acc2 + C.at(i).at(j) * D.at(i).at(j) } else { acc2 }
+    )
+  )
+  // Format optimal assignment as 1-indexed
+  let fstar-display = fstar.map(v => str(v + 1)).join(", ")
+  // Find the highest-flow off-diagonal pair
+  let max-flow = 0
+  let max-fi = 0
+  let max-fj = 0
+  for i in range(n) {
+    for j in range(i + 1, n) {
+      if C.at(i).at(j) > max-flow {
+        max-flow = C.at(i).at(j)
+        max-fi = i
+        max-fj = j
+      }
+    }
+  }
+  let assigned-li = fstar.at(max-fi)
+  let assigned-lj = fstar.at(max-fj)
+  let dist-between = D.at(assigned-li).at(assigned-lj)
+  [
+    #problem-def("QuadraticAssignment")[
+      Given $n$ facilities and $m$ locations ($n <= m$), a flow matrix $C in ZZ^(n times n)$ representing flows between facilities, and a distance matrix $D in ZZ^(m times m)$ representing distances between locations, find an injective assignment $f: {1, dots, n} -> {1, dots, m}$ that minimizes
+      $ sum_(i != j) C_(i j) dot D_(f(i), f(j)). $
+    ][
+    The Quadratic Assignment Problem was introduced by Koopmans and Beckmann (1957) to model the optimal placement of economic activities (facilities) across geographic locations, minimizing total transportation cost weighted by inter-facility flows. It is NP-hard, as shown by Sahni and Gonzalez (1976) via reduction from the Hamiltonian Circuit problem. QAP is widely regarded as one of the hardest combinatorial optimization problems: even moderate instances ($n > 20$) challenge state-of-the-art exact solvers. Best exact approaches use branch-and-bound with Gilmore--Lawler bounds or cutting-plane methods; the best known general algorithm runs in $O^*(n!)$ by exhaustive enumeration of all permutations#footnote[No algorithm significantly improving on brute-force permutation enumeration is known for general QAP.].
+
+    Applications include facility layout planning, keyboard and control panel design, scheduling, VLSI placement, and hospital floor planning. As a special case, when $D$ is a distance matrix on a line (i.e., $D_(k l) = |k - l|$), QAP reduces to the Optimal Linear Arrangement problem.
+
+    *Example.* Consider $n = m = #n$ with flow matrix $C$ and distance matrix $D$:
+    $ C = #to-mat(C), quad D = #to-mat(D). $
+    The identity assignment $f(i) = i$ gives cost #id-cost. The optimal assignment is $f^* = (#fstar-display)$ with cost #cost-star: it places the heavily interacting facilities $F_#(max-fi + 1)$ and $F_#(max-fj + 1)$ (highest flow $= #max-flow$) at locations $L_#(assigned-li + 1)$ and $L_#(assigned-lj + 1)$ (distance $= #dist-between$).
+
+    #figure(
+      canvas(length: 1cm, {
+        import draw: *
+        let fac-x = 0
+        let loc-x = 5
+        let ys = range(n).rev()
+        // Draw facility nodes
+        for i in range(n) {
+          circle((fac-x, ys.at(i)), radius: 0.3, fill: graph-colors.at(0), stroke: 0.8pt + graph-colors.at(0), name: "f" + str(i))
+          content("f" + str(i), text(fill: white, 8pt)[$F_#(i+1)$])
+        }
+        // Draw location nodes
+        for j in range(m) {
+          circle((loc-x, ys.at(j)), radius: 0.3, fill: graph-colors.at(1), stroke: 0.8pt + graph-colors.at(1), name: "l" + str(j))
+          content("l" + str(j), text(fill: white, 8pt)[$L_#(j+1)$])
+        }
+        content((fac-x, n - 0.3), text(9pt, weight: "bold")[Facilities])
+        content((loc-x, m - 0.3), text(9pt, weight: "bold")[Locations])
+        // Draw optimal assignment arrows
+        for (fi, li) in fstar.enumerate() {
+          line("f" + str(fi) + ".east", "l" + str(li) + ".west",
+            mark: (end: "straight"), stroke: 1.2pt + luma(80))
+        }
+        // Highlight highest-flow pair
+        on-layer(-1, {
+          let y0 = calc.min(ys.at(max-fi), ys.at(max-fj)) - 0.55
+          let y1 = calc.max(ys.at(max-fi), ys.at(max-fj)) + 0.55
+          rect((-0.55, y0), (0.55, y1),
+            fill: graph-colors.at(0).transparentize(92%),
+            stroke: (dash: "dashed", paint: graph-colors.at(0).transparentize(50%), thickness: 0.6pt))
+        })
+        content((fac-x, -0.9), text(6pt, fill: luma(100))[flow$(F_#(max-fi + 1), F_#(max-fj + 1)) = #max-flow$])
+      }),
+      caption: [Optimal assignment $f^* = (#fstar-display)$ for the $#n times #m$ QAP instance. Facilities (blue, left) are assigned to locations (red, right) by arrows. Facilities $F_#(max-fi + 1)$ and $F_#(max-fj + 1)$ (highest flow $= #max-flow$) are assigned to locations $L_#(assigned-li + 1)$ and $L_#(assigned-lj + 1)$ (distance $= #dist-between$). Total cost $= #cost-star$.],
+    ) <fig:qap-example>
+    ]
+  ]
+}
+
 #{
   let x = load-model-example("ClosestVectorProblem")
   let basis = x.instance.basis

problemreductions-cli/src/cli.rs

@@ -257,6 +257,7 @@ Flags by problem type:
   LCS                             --strings, --bound [--alphabet-size]
   FAS                             --arcs [--weights] [--num-vertices]
   FVS                             --arcs [--weights] [--num-vertices]
+  QAP                             --matrix (cost), --distance-matrix
   StrongConnectivityAugmentation  --arcs, --candidate-arcs, --bound [--num-vertices]
   FlowShopScheduling              --task-lengths, --deadline [--num-processors]
   StaffScheduling                 --schedules, --requirements, --num-workers, --k
@@ -467,6 +468,9 @@ pub struct CreateArgs {
     /// Directed arcs for directed graph problems (e.g., 0>1,1>2,2>0)
     #[arg(long)]
     pub arcs: Option<String>,
+    /// Distance matrix for QuadraticAssignment (semicolon-separated rows, e.g., "0,1,2;1,0,1;2,1,0")
+    #[arg(long)]
+    pub distance_matrix: Option<String>,
     /// Weighted potential augmentation edges (e.g., 0-2:3,1-3:5)
     #[arg(long)]
     pub potential_edges: Option<String>,

problemreductions-cli/src/commands/create.rs

@@ -84,6 +84,7 @@ fn all_data_flags_empty(args: &CreateArgs) -> bool {
         && args.pattern.is_none()
         && args.strings.is_none()
         && args.arcs.is_none()
+        && args.distance_matrix.is_none()
         && args.candidate_arcs.is_none()
         && args.potential_edges.is_none()
         && args.budget.is_none()
@@ -298,6 +299,7 @@ fn example_for(canonical: &str, graph_type: Option<&str>) -> &'static str {
         "Satisfiability" => "--num-vars 3 --clauses \"1,2;-1,3\"",
         "KSatisfiability" => "--num-vars 3 --clauses \"1,2,3;-1,2,-3\" --k 3",
         "QUBO" => "--matrix \"1,0.5;0.5,2\"",
+        "QuadraticAssignment" => "--matrix \"0,5;5,0\" --distance-matrix \"0,1;1,0\"",
         "SpinGlass" => "--graph 0-1,1-2 --couplings 1,1",
         "KColoring" => "--graph 0-1,1-2,2-0 --k 3",
         "HamiltonianCircuit" => "--graph 0-1,1-2,2-3,3-0",
@@ -1010,6 +1012,54 @@ pub fn create(args: &CreateArgs, out: &OutputConfig) -> Result<()> {
             util::ser_ksat(num_vars, clauses, k)?
         }
 
+        // QuadraticAssignment
+        "QuadraticAssignment" => {
+            let cost_str = args.matrix.as_deref().ok_or_else(|| {
+                anyhow::anyhow!(
+                    "QuadraticAssignment requires --matrix (cost) and --distance-matrix\n\n\
+                     Usage: pred create QAP --matrix \"0,5;5,0\" --distance-matrix \"0,1;1,0\""
+                )
+            })?;
+            let dist_str = args.distance_matrix.as_deref().ok_or_else(|| {
+                anyhow::anyhow!(
+                    "QuadraticAssignment requires --distance-matrix\n\n\
+                     Usage: pred create QAP --matrix \"0,5;5,0\" --distance-matrix \"0,1;1,0\""
+                )
+            })?;
+            let cost_matrix = parse_i64_matrix(cost_str).context("Invalid cost matrix")?;
+            let distance_matrix = parse_i64_matrix(dist_str).context("Invalid distance matrix")?;
+            let n = cost_matrix.len();
+            for (i, row) in cost_matrix.iter().enumerate() {
+                if row.len() != n {
+                    bail!(
+                        "cost matrix must be square: row {i} has {} columns, expected {n}",
+                        row.len()
+                    );
+                }
+            }
+            let m = distance_matrix.len();
+            for (i, row) in distance_matrix.iter().enumerate() {
+                if row.len() != m {
+                    bail!(
+                        "distance matrix must be square: row {i} has {} columns, expected {m}",
+                        row.len()
+                    );
+                }
+            }
+            if n > m {
+                bail!("num_facilities ({n}) must be <= num_locations ({m})");
+            }
+            (
+                ser(
+                    problemreductions::models::algebraic::QuadraticAssignment::new(
+                        cost_matrix,
+                        distance_matrix,
+                    ),
+                )?,
+                resolved_variant.clone(),
+            )
+        }
+
         // QUBO
         "QUBO" => {
             let matrix = parse_matrix(args)?;
@@ -2878,6 +2928,37 @@ fn parse_matrix(args: &CreateArgs) -> Result<Vec<Vec<f64>>> {
         .collect()
 }
 
+/// Parse a semicolon-separated matrix of i64 values.
+/// E.g., "0,5;5,0"
+fn parse_i64_matrix(s: &str) -> Result<Vec<Vec<i64>>> {
+    let matrix: Vec<Vec<i64>> = s
+        .split(';')
+        .enumerate()
+        .map(|(row_idx, row)| {
+            row.trim()
+                .split(',')
+                .enumerate()
+                .map(|(col_idx, v)| {
+                    v.trim().parse::<i64>().map_err(|e| {
+                        anyhow::anyhow!("Invalid value at row {row_idx}, col {col_idx}: {e}")
+                    })
+                })
+                .collect()
+        })
+        .collect::<Result<_>>()?;
+    if let Some(first_len) = matrix.first().map(|r| r.len()) {
+        for (i, row) in matrix.iter().enumerate() {
+            if row.len() != first_len {
+                bail!(
+                    "Ragged matrix: row {i} has {} columns, expected {first_len}",
+                    row.len()
+                );
+            }
+        }
+    }
+    Ok(matrix)
+}
+
 fn parse_potential_edges(args: &CreateArgs) -> Result<Vec<(usize, usize, i32)>> {
     let edges_str = args.potential_edges.as_deref().ok_or_else(|| {
         anyhow::anyhow!("BiconnectivityAugmentation requires --potential-edges (e.g., 0-2:3,1-3:5)")
@@ -3552,6 +3633,7 @@ mod tests {
             pattern: None,
             strings: None,
             arcs: None,
+            distance_matrix: None,
             potential_edges: None,
             budget: None,
             candidate_arcs: None,

src/lib.rs

@@ -42,7 +42,7 @@ pub mod variant;
 /// Prelude module for convenient imports.
 pub mod prelude {
     // Problem types
-    pub use crate::models::algebraic::{BMF, QUBO};
+    pub use crate::models::algebraic::{QuadraticAssignment, BMF, QUBO};
     pub use crate::models::formula::{CNFClause, CircuitSAT, KSatisfiability, Satisfiability};
     pub use crate::models::graph::{
         BalancedCompleteBipartiteSubgraph, BicliqueCover, BiconnectivityAugmentation,

src/models/algebraic/mod.rs

@@ -5,15 +5,18 @@
 //! - [`ILP`]: Integer Linear Programming
 //! - [`ClosestVectorProblem`]: Closest Vector Problem (minimize lattice distance)
 //! - [`BMF`]: Boolean Matrix Factorization
+//! - [`QuadraticAssignment`]: Quadratic Assignment Problem
 
 pub(crate) mod bmf;
 pub(crate) mod closest_vector_problem;
 pub(crate) mod ilp;
+pub(crate) mod quadratic_assignment;
 pub(crate) mod qubo;
 
 pub use bmf::BMF;
 pub use closest_vector_problem::{ClosestVectorProblem, VarBounds};
 pub use ilp::{Comparison, LinearConstraint, ObjectiveSense, VariableDomain, ILP};
+pub use quadratic_assignment::QuadraticAssignment;
 pub use qubo::QUBO;
 
 #[cfg(feature = "example-db")]
@@ -23,5 +26,6 @@ pub(crate) fn canonical_model_example_specs() -> Vec<crate::example_db::specs::M
     specs.extend(ilp::canonical_model_example_specs());
     specs.extend(closest_vector_problem::canonical_model_example_specs());
     specs.extend(bmf::canonical_model_example_specs());
+    specs.extend(quadratic_assignment::canonical_model_example_specs());
     specs
 }