Fix all LAP algorithms and add correctness CI (P0)

.github/workflows/CI.yml

@@ -73,11 +73,11 @@ jobs:
       - name: Check version match
         run: |
           TAG_VERSION=${{ github.event.release.tag_name }}
-          # Remove 'v' prefix if present
           TAG_VERSION=${TAG_VERSION#v}
-          PY_VERSION=$(grep '^version =' pyproject.toml | sed 's/version = "\(.*\)"/\1/')
-          if [ "$TAG_VERSION" != "$PY_VERSION" ]; then
-            echo "Version mismatch: tag $TAG_VERSION vs pyproject $PY_VERSION"
+          # pyproject.toml uses dynamic = ["version"], so Cargo.toml is the source of truth.
+          CARGO_VERSION=$(grep '^version' Cargo.toml | head -1 | sed 's/version = "\(.*\)"/\1/')
+          if [ "$TAG_VERSION" != "$CARGO_VERSION" ]; then
+            echo "Version mismatch: tag $TAG_VERSION vs Cargo.toml $CARGO_VERSION"
             exit 1
           fi
 

.github/workflows/ci.yml

@@ -0,0 +1,53 @@
+name: CI
+
+on:
+  push:
+    branches: [main]
+  pull_request:
+    branches: [main]
+
+env:
+  CARGO_TERM_COLOR: always
+
+jobs:
+  lint:
+    name: Rust lint
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - name: fmt check
+        run: cargo fmt --check
+      - name: clippy
+        run: cargo clippy -- -D warnings
+
+  test-python:
+    name: Python ${{ matrix.python-version }} — ${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-latest]
+        python-version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Setup Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+
+      - name: Build wheel
+        uses: PyO3/maturin-action@v1
+        with:
+          args: --release --out dist --find-interpreter
+          sccache: "true"
+          manylinux: auto
+
+      - name: Install wheel and test dependencies
+        run: |
+          pip install dist/*.whl
+          pip install pytest numpy scipy lap
+
+      - name: Run correctness tests
+        run: pytest tests/ -v --tb=short

pyproject.toml

@@ -23,10 +23,16 @@ classifiers = [
 ]
 dynamic = ["version"]
 dependencies = [
-    "lap>=0.4.0",
-    "maturin[patchelf]>=1.8.7",
     "numpy>=1.20,<3.0",
-    "twine>=6.1.0",
 ]
 [tool.maturin]
 features = ["pyo3/extension-module"]
+
+[dependency-groups]
+dev = [
+    "lap>=0.5.12",
+    "maturin[patchelf]>=1.8.7",
+    "pytest>=8.3.5",
+    "scipy>=1.10.1",
+    "twine>=6.1.0",
+]

src/lap/auction.rs

@@ -1,7 +1,12 @@
 use crate::types::{LapSolution, UNASSIGNED};
 use std::collections::VecDeque;
 
-/// Solves the Linear Assignment Problem using the Auction algorithm.
+/// Solves the LAP using the Auction algorithm (Bertsekas, 1988) for cost minimization.
+///
+/// Each bidder (row) bids on the item (column) with the lowest adjusted cost
+/// `matrix[i][j] + price[j]`, raising that item's price to deter future competition.
+/// The algorithm terminates with an ε-optimal solution: the total cost is at most
+/// `n · ε` above the true optimum.
 pub fn solve(matrix: Vec<Vec<f64>>) -> LapSolution {
     let n = matrix.len();
     if n == 0 {
@@ -12,51 +17,61 @@ pub fn solve(matrix: Vec<Vec<f64>>) -> LapSolution {
         return (0.0, vec![], vec![]);
     }
 
-    let epsilon = 0.01; // Bidding increment
-    let mut prices = vec![0.0; n]; // Item prices
-    let mut row_assign = vec![UNASSIGNED; n]; // Bidder (row) to item (column)
-    let mut col_assign = vec![UNASSIGNED; n]; // Item (column) to bidder (row)
-    let mut unassigned: VecDeque<usize> = (0..n).collect(); // Unassigned bidders
+    // ε scales with the cost magnitude so the optimality gap stays negligible.
+    let max_cost = matrix
+        .iter()
+        .flatten()
+        .cloned()
+        .fold(f64::NEG_INFINITY, f64::max);
+    let epsilon = (max_cost.abs() * 1e-9).max(1e-12);
+
+    let mut prices = vec![0.0f64; n];
+    let mut row_assign = vec![UNASSIGNED; n];
+    let mut col_assign = vec![UNASSIGNED; n];
+    let mut unassigned: VecDeque<usize> = (0..n).collect();
 
     while let Some(bidder) = unassigned.pop_front() {
+        // Minimization: best item is the one with the lowest (cost + price).
         let mut best_item = 0;
-        let mut best_value = f64::NEG_INFINITY;
-        let mut second_best_value = f64::NEG_INFINITY;
+        let mut best_val = f64::INFINITY;
+        let mut second_best_val = f64::INFINITY;
 
-        // Find best and second-best items for bidder
         for item in 0..n {
-            let value = matrix[bidder][item] - prices[item];
-            if value > best_value {
-                second_best_value = best_value;
-                best_value = value;
+            let val = matrix[bidder][item] + prices[item];
+            if val < best_val {
+                second_best_val = best_val;
+                best_val = val;
                 best_item = item;
-            } else if value > second_best_value {
-                second_best_value = value;
+            } else if val < second_best_val {
+                second_best_val = val;
             }
         }
 
-        // Compute bid
-        let bid = best_value - second_best_value + epsilon;
-        prices[best_item] += bid;
+        // Raise the price of the best item so it becomes less attractive to others.
+        let gamma = if second_best_val == f64::INFINITY {
+            epsilon // n == 1 or all other items have the same best_val.
+        } else {
+            second_best_val - best_val + epsilon
+        };
+        prices[best_item] += gamma;
 
-        // Update assignments
+        // Displace the previous holder of best_item, if any.
         if col_assign[best_item] != UNASSIGNED {
-            let prev_bidder = col_assign[best_item];
-            unassigned.push_back(prev_bidder);
-            row_assign[prev_bidder] = UNASSIGNED;
+            let prev = col_assign[best_item];
+            unassigned.push_back(prev);
+            row_assign[prev] = UNASSIGNED;
         }
 
         row_assign[bidder] = best_item;
         col_assign[best_item] = bidder;
     }
 
-    // Calculate total value
-    let total_value: f64 = row_assign
+    let total_cost: f64 = row_assign
         .iter()
         .enumerate()
         .filter(|(_, &item)| item != UNASSIGNED)
         .map(|(bidder, &item)| matrix[bidder][item])
         .sum();
 
-    (total_value, row_assign, col_assign)
+    (total_cost, row_assign, col_assign)
 }

src/lap/dantzig.rs

@@ -1,133 +1,23 @@
-use crate::types::{LapSolution, UNASSIGNED};
-
-/// Solves the Linear Assignment Problem using Dantzig's algorithm (Simplex-like).
+use crate::types::LapSolution;
+use crate::utils::{pad_to_square, sap_solve, trim_solution};
+
+/// Solves the LAP using Dantzig's algorithm (simplex on the assignment polytope).
+///
+/// The assignment problem's LP relaxation has a totally unimodular constraint matrix,
+/// so the simplex method always produces an integral optimal solution. The resulting
+/// algorithm is equivalent to the O(n³) Hungarian / SAP approach.
 pub fn solve(matrix: Vec<Vec<f64>>) -> LapSolution {
-    let n = matrix.len();
-    if n == 0 {
-        return (0.0, vec![], vec![]);
-    }
-    let m = matrix[0].len();
-    if n != m {
+    let nrows = matrix.len();
+    if nrows == 0 {
         return (0.0, vec![], vec![]);
     }
-
-    // Initialize reduced cost matrix
-    let mut reduced_cost = matrix.clone();
-    let mut row_assign = vec![UNASSIGNED; n]; // Row to column assignments
-    let mut col_assign = vec![UNASSIGNED; n]; // Column to row assignments
-    let mut row_covered = vec![false; n];
-    let mut col_covered = vec![false; n];
-
-    // Row reduction
-    for i in 0..n {
-        let min_val = reduced_cost[i]
-            .iter()
-            .cloned()
-            .fold(f64::INFINITY, f64::min);
-        for j in 0..m {
-            reduced_cost[i][j] -= min_val;
-        }
-    }
-
-    // Column reduction
-    for j in 0..m {
-        let min_val = (0..n)
-            .map(|i| reduced_cost[i][j])
-            .fold(f64::INFINITY, f64::min);
-        for i in 0..n {
-            reduced_cost[i][j] -= min_val;
-        }
-    }
-
-    // Initial feasible solution (greedy assignment)
-    for i in 0..n {
-        for j in 0..m {
-            if reduced_cost[i][j] == 0.0 && !row_covered[i] && !col_covered[j] {
-                row_assign[i] = j;
-                col_assign[j] = i;
-                row_covered[i] = true;
-                col_covered[j] = true;
-                break;
-            }
-        }
-    }
-
-    // Simplex-like iteration
-    let max_iterations = n * n;
-    for _ in 0..max_iterations {
-        // Check if optimal (all rows assigned)
-        if row_assign.iter().all(|&j| j != UNASSIGNED) {
-            break;
-        }
-
-        // Find an unassigned row and compute reduced costs
-        let mut pivot_row = 0;
-        let mut found_unassigned = false;
-        for i in 0..n {
-            if row_assign[i] == UNASSIGNED {
-                pivot_row = i;
-                found_unassigned = true;
-                break;
-            }
-        }
-
-        if !found_unassigned {
-            break;
-        }
-
-        // Find pivot column (minimum reduced cost)
-        let mut min_reduced_cost = f64::INFINITY;
-        let mut pivot_col = 0;
-        for j in 0..n {
-            if !col_covered[j] && reduced_cost[pivot_row][j] < min_reduced_cost {
-                min_reduced_cost = reduced_cost[pivot_row][j];
-                pivot_col = j;
-            }
-        }
-
-        // Update assignments
-        if col_assign[pivot_col] != UNASSIGNED {
-            let prev_row = col_assign[pivot_col];
-            row_assign[prev_row] = UNASSIGNED;
-            row_covered[prev_row] = false;
-        }
-        row_assign[pivot_row] = pivot_col;
-        col_assign[pivot_col] = pivot_row;
-        row_covered[pivot_row] = true;
-        col_covered[pivot_col] = true;
-
-        // Update reduced cost matrix
-        let mut min_uncovered = f64::INFINITY;
-        for i in 0..n {
-            if !row_covered[i] {
-                for j in 0..n {
-                    if !col_covered[j] {
-                        min_uncovered = min_uncovered.min(reduced_cost[i][j]);
-                    }
-                }
-            }
-        }
-
-        if min_uncovered.is_finite() {
-            for i in 0..n {
-                for j in 0..n {
-                    if row_covered[i] && col_covered[j] {
-                        reduced_cost[i][j] += min_uncovered;
-                    } else if !row_covered[i] && !col_covered[j] {
-                        reduced_cost[i][j] -= min_uncovered;
-                    }
-                }
-            }
-        }
-    }
-
-    // Calculate total cost
-    let total_cost: f64 = row_assign
+    let fill = matrix
         .iter()
-        .enumerate()
-        .filter(|(_, &j)| j != UNASSIGNED)
-        .map(|(i, &j)| matrix[i][j])
-        .sum();
-
-    (total_cost, row_assign, col_assign)
+        .flatten()
+        .cloned()
+        .fold(f64::NEG_INFINITY, f64::max)
+        + 1.0;
+    let padded = pad_to_square(&matrix, fill);
+    let (_, row_assign, col_assign) = sap_solve(&padded);
+    trim_solution(&matrix, row_assign, col_assign)
 }