Fix #122: Add SteinerTree model

docs/paper/reductions.typ

@@ -100,6 +100,7 @@
   "MinimumFeedbackVertexSet": [Minimum Feedback Vertex Set],
   "ShortestCommonSupersequence": [Shortest Common Supersequence],
   "MinimumSumMulticenter": [Minimum Sum Multicenter],
+  "SteinerTree": [Steiner Tree],
   "SubgraphIsomorphism": [Subgraph Isomorphism],
   "PartitionIntoTriangles": [Partition Into Triangles],
   "FlowShopScheduling": [Flow Shop Scheduling],
@@ -786,6 +787,72 @@ Graph Partitioning is a core NP-hard problem arising in VLSI design, parallel co
     ]
   ]
 }
+#{
+  let x = load-model-example("SteinerTree")
+  let nv = graph-num-vertices(x.instance)
+  let ne = graph-num-edges(x.instance)
+  let edges = x.instance.graph.inner.edges.map(e => (e.at(0), e.at(1)))
+  let weights = x.instance.edge_weights
+  let terminals = x.instance.terminals
+  let sol = x.optimal.at(0)
+  let tree-edge-indices = sol.config.enumerate().filter(((i, v)) => v == 1).map(((i, _)) => i)
+  let tree-edges = tree-edge-indices.map(i => edges.at(i))
+  let cost = sol.metric.Valid
+  // Steiner vertices: in tree but not terminals
+  let tree-verts = tree-edges.map(e => (e.at(0), e.at(1))).fold((), (acc, pair) => {
+    let (u, v) = pair
+    let acc2 = if acc.contains(u) { acc } else { acc + (u,) }
+    if acc2.contains(v) { acc2 } else { acc2 + (v,) }
+  })
+  let steiner-verts = tree-verts.filter(v => not terminals.contains(v))
+  [
+    #problem-def("SteinerTree")[
+      Given an undirected graph $G = (V, E)$ with edge weights $w: E -> RR_(>= 0)$ and a set of terminal vertices $T subset.eq V$ with $|T| >= 2$, find a tree $S = (V_S, E_S)$ in $G$ such that $T subset.eq V_S$, minimizing $sum_(e in E_S) w(e)$. Vertices in $V_S backslash T$ are called _Steiner vertices_.
+    ][
+    One of Karp's 21 NP-complete problems @karp1972, foundational in network design with applications in telecommunications backbone routing, VLSI chip interconnect, pipeline planning, and phylogenetic tree construction. When $T = V$, the problem reduces to the minimum spanning tree (polynomial). The NP-hardness arises from choosing which Steiner vertices to include.
+
+    The best known exact algorithm runs in $O^*(3^(|T|) dot n + 2^(|T|) dot n^2)$ time via Dreyfus--Wagner dynamic programming over terminal subsets @dreyfuswagner1971. Byrka _et al._ achieved a $ln(4) + epsilon approx 1.39$-approximation @byrka2013; the classic 2-approximation uses the minimum spanning tree of the terminal distance graph.
+
+    // Find the unique direct terminal-terminal edge (both endpoints in T, not in the optimal tree)
+    #let terminal-set = terminals
+    #let direct-tt-edges = edges.enumerate().filter(((i, e)) => {
+      terminal-set.contains(e.at(0)) and terminal-set.contains(e.at(1)) and not tree-edge-indices.contains(i)
+    })
+    #let tt-edge = direct-tt-edges.at(0)
+    #let tt-idx = tt-edge.at(0)
+    #let tt-u = tt-edge.at(1).at(0)
+    #let tt-v = tt-edge.at(1).at(1)
+
+    *Example.* Consider $G$ with $n = #nv$ vertices, $m = #ne$ edges, and terminals $T = {#terminals.map(t => $v_#t$).join(", ")}$. The optimal Steiner tree uses edges ${#tree-edges.map(e => $(v_#(e.at(0)), v_#(e.at(1)))$).join(", ")}$ with Steiner vertices ${#steiner-verts.map(v => $v_#v$).join(", ")}$ acting as relay points. The total cost is #tree-edge-indices.map(i => $#(weights.at(i))$).join($+$) $= #cost$. Note the only direct terminal--terminal edge $(v_#tt-u, v_#tt-v)$ has weight #weights.at(tt-idx), equaling the entire Steiner tree cost.
+
+    #figure({
+      // Layout: v0 top-left, v1 top-center, v2 top-right, v3 bottom-center, v4 bottom-right
+      let verts = ((0, 1.2), (1.2, 1.2), (2.4, 1.2), (1.2, 0), (2.4, 0))
+      canvas(length: 1cm, {
+        for (idx, (u, v)) in edges.enumerate() {
+          let on-tree = tree-edge-indices.contains(idx)
+          g-edge(verts.at(u), verts.at(v),
+            stroke: if on-tree { 2pt + graph-colors.at(0) } else { 1pt + luma(200) })
+          let mx = (verts.at(u).at(0) + verts.at(v).at(0)) / 2
+          let my = (verts.at(u).at(1) + verts.at(v).at(1)) / 2
+          let dx = if u == 0 and v == 3 { -0.3 } else if u == 2 and v == 3 { 0.3 } else { 0 }
+          let dy = if u == 0 and v == 1 { 0.2 } else if u == 1 and v == 2 { 0.2 } else if u == 2 and v == 4 { 0.3 } else { 0 }
+          draw.content((mx + dx, my + dy), text(7pt, fill: luma(80))[#weights.at(idx)])
+        }
+        for (k, pos) in verts.enumerate() {
+          let is-terminal = terminals.contains(k)
+          g-node(pos, name: "v" + str(k),
+            fill: if is-terminal { graph-colors.at(0) } else { white },
+            stroke: if is-terminal { none } else { 1pt + graph-colors.at(0) },
+            label: text(fill: if is-terminal { white } else { black })[$v_#k$])
+        }
+      })
+    },
+    caption: [Steiner tree on #nv vertices with terminals $T = {#terminals.map(t => $v_#t$).join(", ")}$ (filled blue). Steiner vertices #steiner-verts.map(v => $v_#v$).join(", ") (outlined) relay connections. Blue edges form the optimal tree with cost #cost.],
+    ) <fig:steiner-tree>
+    ]
+  ]
+}
 #problem-def("OptimalLinearArrangement")[
   Given an undirected graph $G=(V,E)$ and a non-negative integer $K$, is there a bijection $f: V -> {0, 1, dots, |V|-1}$ such that $sum_({u,v} in E) |f(u) - f(v)| <= K$?
 ][

docs/paper/references.bib

@@ -459,6 +459,28 @@ @article{alber2004
   doi     = {10.1016/j.jalgor.2003.10.001}
 }
 
+@article{dreyfuswagner1971,
+  author  = {S. E. Dreyfus and R. A. Wagner},
+  title   = {The Steiner Problem in Graphs},
+  journal = {Networks},
+  volume  = {1},
+  number  = {3},
+  pages   = {195--207},
+  year    = {1971},
+  doi     = {10.1002/net.3230010302}
+}
+
+@article{byrka2013,
+  author  = {Jarosław Byrka and Fabrizio Grandoni and Thomas Rothvoß and Laura Sanità},
+  title   = {Steiner Tree Approximation via Iterative Randomized Rounding},
+  journal = {Journal of the ACM},
+  volume  = {60},
+  number  = {1},
+  pages   = {1--33},
+  year    = {2013},
+  doi     = {10.1145/2432622.2432628}
+}
+
 @article{horowitz1974,
   author  = {Ellis Horowitz and Sartaj Sahni},
   title   = {Computing Partitions with Applications to the Knapsack Problem},

docs/src/cli.md

@@ -45,6 +45,9 @@ pred create MIS --graph 0-1,1-2,2-3 -o problem.json
 # Create a weighted instance (variant auto-upgrades to i32)
 pred create MIS --graph 0-1,1-2,2-3 --weights 3,1,2,1 -o weighted.json
 
+# Create a Steiner Tree instance
+pred create SteinerTree --graph 0-1,0-3,1-2,1-3,2-3,2-4,3-4 --edge-weights 2,5,2,1,5,6,1 --terminals 0,2,4 -o steiner.json
+
 # Or start from a canonical model example
 pred create --example MIS/SimpleGraph/i32 -o example.json
 
@@ -272,6 +275,7 @@ pred create QUBO --matrix "1,0.5;0.5,2" -o qubo.json
 pred create KColoring --k 3 --graph 0-1,1-2,2-0 -o kcol.json
 pred create SpinGlass --graph 0-1,1-2 -o sg.json
 pred create MaxCut --graph 0-1,1-2,2-0 -o maxcut.json
+pred create SteinerTree --graph 0-1,0-3,1-2,1-3,2-3,2-4,3-4 --edge-weights 2,5,2,1,5,6,1 --terminals 0,2,4 -o steiner.json
 pred create Factoring --target 15 --bits-m 4 --bits-n 4 -o factoring.json
 pred create Factoring --target 21 --bits-m 3 --bits-n 3 -o factoring2.json
 pred create X3C --universe 9 --sets "0,1,2;0,2,4;3,4,5;3,5,7;6,7,8;1,4,6;2,5,8" -o x3c.json

docs/src/reductions/problem_schemas.json

@@ -636,6 +636,27 @@
       }
     ]
   },
+  {
+    "name": "SteinerTree",
+    "description": "Find minimum weight tree connecting terminal vertices",
+    "fields": [
+      {
+        "name": "graph",
+        "type_name": "G",
+        "description": "The underlying graph G=(V,E)"
+      },
+      {
+        "name": "edge_weights",
+        "type_name": "Vec<W>",
+        "description": "Edge weights w: E -> R"
+      },
+      {
+        "name": "terminals",
+        "type_name": "Vec<usize>",
+        "description": "Terminal vertices T that must be connected"
+      }
+    ]
+  },
   {
     "name": "SubgraphIsomorphism",
     "description": "Determine if host graph G contains a subgraph isomorphic to pattern graph H",

docs/src/reductions/reduction_graph.json

@@ -491,6 +491,26 @@
       "doc_path": "models/graph/struct.SpinGlass.html",
       "complexity": "2^num_spins"
     },
+    {
+      "name": "SteinerTree",
+      "variant": {
+        "graph": "SimpleGraph",
+        "weight": "One"
+      },
+      "category": "graph",
+      "doc_path": "models/graph/struct.SteinerTree.html",
+      "complexity": "3^num_terminals * num_vertices + 2^num_terminals * num_vertices^2"
+    },
+    {
+      "name": "SteinerTree",
+      "variant": {
+        "graph": "SimpleGraph",
+        "weight": "i32"
+      },
+      "category": "graph",
+      "doc_path": "models/graph/struct.SteinerTree.html",
+      "complexity": "3^num_terminals * num_vertices + 2^num_terminals * num_vertices^2"
+    },
     {
       "name": "SubgraphIsomorphism",
       "variant": {},
@@ -713,7 +733,7 @@
     },
     {
       "source": 21,
-      "target": 56,
+      "target": 58,
       "overhead": [
         {
           "field": "num_elements",
@@ -1341,7 +1361,7 @@
       "doc_path": "rules/spinglass_casts/index.html"
     },
     {
-      "source": 57,
+      "source": 59,
       "target": 12,
       "overhead": [
         {
@@ -1356,7 +1376,7 @@
       "doc_path": "rules/travelingsalesman_ilp/index.html"
     },
     {
-      "source": 57,
+      "source": 59,
       "target": 49,
       "overhead": [
         {

problemreductions-cli/src/cli.rs

@@ -233,6 +233,7 @@ Flags by problem type:
   X3C (ExactCoverBy3Sets)         --universe, --sets (3 elements each)
   BicliqueCover                   --left, --right, --biedges, --k
   BMF                             --matrix (0/1), --rank
+  SteinerTree                     --graph, --edge-weights, --terminals
   CVP                             --basis, --target-vec [--bounds]
   OptimalLinearArrangement        --graph, --bound
   RuralPostman (RPP)              --graph, --edge-weights, --required-edges, --bound
@@ -367,6 +368,9 @@ pub struct CreateArgs {
     /// Variable bounds for CVP as "lower,upper" (e.g., "-10,10") [default: -10,10]
     #[arg(long, allow_hyphen_values = true)]
     pub bounds: Option<String>,
+    /// Terminal vertices for SteinerTree (comma-separated indices, e.g., "0,2,4")
+    #[arg(long)]
+    pub terminals: Option<String>,
     /// Tree edge list for IsomorphicSpanningTree (e.g., 0-1,1-2,2-3)
     #[arg(long)]
     pub tree: Option<String>,

problemreductions-cli/src/commands/create.rs

@@ -18,7 +18,7 @@ use problemreductions::topology::{
     UnitDiskGraph,
 };
 use serde::Serialize;
-use std::collections::BTreeMap;
+use std::collections::{BTreeMap, BTreeSet};
 
 /// Check if all data flags are None (no problem-specific input provided).
 fn all_data_flags_empty(args: &CreateArgs) -> bool {
@@ -51,6 +51,7 @@ fn all_data_flags_empty(args: &CreateArgs) -> bool {
         && args.basis.is_none()
         && args.target_vec.is_none()
         && args.bounds.is_none()
+        && args.terminals.is_none()
         && args.tree.is_none()
         && args.required_edges.is_none()
         && args.bound.is_none()
@@ -239,6 +240,7 @@ fn example_for(canonical: &str, graph_type: Option<&str>) -> &'static str {
         }
         "PartitionIntoTriangles" => "--graph 0-1,1-2,0-2",
         "Factoring" => "--target 15 --m 4 --n 4",
+        "SteinerTree" => "--graph 0-1,1-2,1-3,3-4 --edge-weights 2,2,1,1 --terminals 0,2,4",
         "OptimalLinearArrangement" => "--graph 0-1,1-2,2-3 --bound 5",
         "MinimumFeedbackArcSet" => "--arcs \"0>1,1>2,2>0\"",
         "RuralPostman" => {
@@ -360,6 +362,19 @@ pub fn create(args: &CreateArgs, out: &OutputConfig) -> Result<()> {
             create_vertex_weight_problem(args, canonical, graph_type, &resolved_variant)?
         }
 
+        // SteinerTree (graph + edge weights + terminals)
+        "SteinerTree" => {
+            let (graph, _) = parse_graph(args).map_err(|e| {
+                anyhow::anyhow!(
+                    "{e}\n\nUsage: pred create SteinerTree --graph 0-1,1-2,1-3,3-4 --edge-weights 2,2,1,1 --terminals 0,2,4"
+                )
+            })?;
+            let edge_weights = parse_edge_weights(args, graph.num_edges())?;
+            let terminals = parse_terminals(args, graph.num_vertices())?;
+            let data = ser(SteinerTree::new(graph, edge_weights, terminals))?;
+            (data, resolved_variant.clone())
+        }
+
         // Graph partitioning (graph only, no weights)
         "GraphPartitioning" => {
             let (graph, _) = parse_graph(args).map_err(|e| {
@@ -1298,6 +1313,32 @@ fn parse_vertex_weights(args: &CreateArgs, num_vertices: usize) -> Result<Vec<i3
     }
 }
 
+/// Parse `--terminals` as comma-separated vertex indices.
+fn parse_terminals(args: &CreateArgs, num_vertices: usize) -> Result<Vec<usize>> {
+    let s = args
+        .terminals
+        .as_deref()
+        .ok_or_else(|| anyhow::anyhow!("SteinerTree requires --terminals (e.g., \"0,2,4\")"))?;
+    let terminals: Vec<usize> = s
+        .split(',')
+        .map(|t| t.trim().parse::<usize>())
+        .collect::<std::result::Result<Vec<_>, _>>()
+        .context("invalid terminal index")?;
+    for &t in &terminals {
+        anyhow::ensure!(
+            t < num_vertices,
+            "terminal {t} >= num_vertices ({num_vertices})"
+        );
+    }
+    let distinct_terminals: BTreeSet<_> = terminals.iter().copied().collect();
+    anyhow::ensure!(
+        distinct_terminals.len() == terminals.len(),
+        "terminals must be distinct"
+    );
+    anyhow::ensure!(terminals.len() >= 2, "at least 2 terminals required");
+    Ok(terminals)
+}
+
 /// Parse `--edge-weights` as edge weights (i32), defaulting to all 1s.
 fn parse_edge_weights(args: &CreateArgs, num_edges: usize) -> Result<Vec<i32>> {
     match &args.edge_weights {
@@ -1680,6 +1721,34 @@ fn create_random(
             (data, variant)
         }
 
+        // SteinerTree
+        "SteinerTree" => {
+            anyhow::ensure!(
+                num_vertices >= 2,
+                "SteinerTree random generation requires --num-vertices >= 2"
+            );
+            let edge_prob = args.edge_prob.unwrap_or(0.5);
+            if !(0.0..=1.0).contains(&edge_prob) {
+                bail!("--edge-prob must be between 0.0 and 1.0");
+            }
+            let mut state = util::lcg_init(args.seed);
+            let graph = util::create_random_graph(num_vertices, edge_prob, Some(state));
+            // Advance state past the graph generation
+            for _ in 0..num_vertices * num_vertices {
+                util::lcg_step(&mut state);
+            }
+            let edge_weights: Vec<i32> = (0..graph.num_edges())
+                .map(|_| (util::lcg_step(&mut state) * 9.0) as i32 + 1)
+                .collect();
+            let num_terminals = std::cmp::max(2, num_vertices * 2 / 5);
+            let terminals = util::lcg_choose(&mut state, num_vertices, num_terminals);
+            let variant = variant_map(&[("graph", "SimpleGraph"), ("weight", "i32")]);
+            (
+                ser(SteinerTree::new(graph, edge_weights, terminals))?,
+                variant,
+            )
+        }
+
         // SpinGlass
         "SpinGlass" => {
             let edge_prob = args.edge_prob.unwrap_or(0.5);
@@ -1730,7 +1799,7 @@ fn create_random(
             "Random generation is not supported for {canonical}. \
              Supported: graph-based problems (MIS, MVC, MaxCut, MaxClique, \
              MaximumMatching, MinimumDominatingSet, SpinGlass, KColoring, TravelingSalesman, \
-             OptimalLinearArrangement, HamiltonianPath)"
+             SteinerTree, OptimalLinearArrangement, HamiltonianPath)"
         ),
     };
 

problemreductions-cli/src/util.rs

@@ -214,6 +214,20 @@ pub fn create_random_float_positions(num_vertices: usize, seed: Option<u64>) ->
         .collect()
 }
 
+/// Choose `k` distinct elements from `0..n` using Fisher-Yates partial shuffle.
+/// Returns a sorted vector of chosen indices.
+pub fn lcg_choose(state: &mut u64, n: usize, k: usize) -> Vec<usize> {
+    assert!(k <= n, "k={k} exceeds n={n}");
+    let mut indices: Vec<usize> = (0..n).collect();
+    for i in 0..k {
+        let j = i + (lcg_step(state) * (n - i) as f64) as usize % (n - i);
+        indices.swap(i, j);
+    }
+    let mut chosen: Vec<usize> = indices[..k].to_vec();
+    chosen.sort_unstable();
+    chosen
+}
+
 // ---------------------------------------------------------------------------
 // Small shared helpers
 // ---------------------------------------------------------------------------