Metadata-Version: 2.4
Name: acdc-search
Version: 0.1.0
Summary: AC/DC search - alternating continuous and discrete combinatorial graph matching
Project-URL: Homepage, https://github.com/flyconnectome/acdc_py
Project-URL: Source, https://github.com/flyconnectome/acdc_py
Project-URL: Issues, https://github.com/flyconnectome/acdc_py/issues
Author-email: Philipp Schlegel <observing@web.de>
License-Expression: MIT
License-File: LICENSE
License-File: NOTICE
Keywords: connectomics,frank-wolfe,graph-matching,quadratic-assignment
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Mathematics
Requires-Python: >=3.11
Requires-Dist: numpy>=1.23.2
Requires-Dist: scipy>=1.9.2
Provides-Extra: test
Requires-Dist: pytest>=7.0; extra == 'test'
Description-Content-Type: text/markdown

# AC/DC search (Python)

A Python translation of the MATLAB reference implementation of **AC/DC search** —
the **A**lternating **C**ontinuous and **D**iscrete **C**ombinatorial optimization
behind the winning solution to the FlyWire Ventral Nerve Cord Matching
Challenge.

> D. D. Lee, A. Matsliah & L. K. Saul, *"AC/DC search: behind the winning
> solution to the FlyWire graph-matching challenge"*, Transactions on Machine
> Learning Research (01/2026).
> [OpenReview](https://openreview.net/forum?id=8MjCOMyaDf)

## What it does

Given two weighted directed graphs `A` and `B` (same number of nodes), AC/DC
searches for the permutation `P` (a one-to-one node correspondence) that
maximizes the **overlapping edge weight**

```
J(P) = sum_ij  min( A_ij , B_{p(i), p(j)} )
```

where `p(i)` is the node of `B` matched to node `i` of `A`. It alternates:

- a **continuous** phase — Frank–Wolfe optimization over the doubly stochastic
  relaxation (Birkhoff polytope), with an exact closed-form line search and a
  linear-assignment subproblem at each step; and
- a **discrete** phase — greedy pairwise-swap search over permutations, applying
  the exact-gain swaps that improve the score.

Each phase warm-starts the other.

## Install

```bash
pip install acdc-search
```

The distribution is named `acdc-search` (the names `acdc` and `acdc-py` are
taken on PyPI by unrelated projects); the import name is `acdc`. It is pure
Python and needs only NumPy and SciPy.

For a development install from a checkout:

```bash
pip install -e ".[test]"        # numpy, scipy + pytest
```

## Usage

```python
import numpy as np
import scipy.sparse as sp
from acdc import acdc_match, score, matrix_to_perm

# --- build two graphs A, B as scipy.sparse n x n matrices ---
rng = np.random.default_rng(0)
A = sp.random(200, 200, density=0.05, random_state=0).tocsr()

# (toy example) B is A relabeled by a hidden permutation we will recover
perm = rng.permutation(200)
P = sp.csc_matrix((np.ones(200), (np.arange(200), perm)), shape=(200, 200))
B = (P.T @ A @ P).tocsr()

# --- run the full AC/DC algorithm ---
M = acdc_match(A, B, max_iter=5, num_frank_wolfe=10)

print("score:", score(M, A, B), " (optimum:", A.sum(), ")")
matching = matrix_to_perm(M)          # 0-based: node i of A -> matching[i] of B
```

### High-level functions

| Function | MATLAB equivalent | Description |
|---|---|---|
| `acdc_match(A, B, P0=None, max_iter=5, num_frank_wolfe=10, max_swap=inf, solver='dense')` | `main_acdc.m` | Full alternating algorithm. |
| `frank_wolfe_search(A, B, P0=None, num_updates=40, solver='dense')` | `main_continuous.m` | Continuous (Frank–Wolfe) phase only. |
| `greedy_match(A, B, P0=None, max_swap=inf)` | `main_discrete.m` | Discrete (greedy-swap) phase only. |

All three return a sparse permutation matrix.

**Progress output.** The phases emit their MATLAB-style tables through the
`acdc` logger at `INFO`. `verbose=True` (the default) attaches a stdout handler
for the duration of the call, so they just appear; `verbose=False` leaves your
logging configuration alone, which means an application that has configured
`logging` at `INFO` still receives them and one that has not stays silent:

```python
import logging
logging.basicConfig(level=logging.INFO)         # or logging.getLogger("acdc")
M = acdc_match(A, B, verbose=False)             # tables go to your handlers
```

**Inputs.** `A` and `B` may be a SciPy sparse matrix, a dense `numpy.ndarray`,
or an edge-list tuple `(rows, cols, weights)` / `(rows, cols, weights, n)` with
**0-based** node indices. `P0` may be a permutation vector, a sparse permutation
matrix, or `None` (identity start).

### Public API

`acdc` exports exactly six names: the three entry points above, plus `score`
to evaluate a matching and `perm_to_matrix` / `matrix_to_perm` to convert
between a permutation matrix and a 0-based permutation vector.

The individual algorithm phases live in the submodules that mirror the MATLAB
files — `acdc.objective` (`compute_gradient`), `acdc.frank_wolfe`
(`do_frank_wolfe`), `acdc.swaps` (`evaluate_swaps`, `make_swaps`,
`greedy_search`), `acdc.matching` (`permutation_match`, `as_graph`, …) — and
can be imported from there. They are **internal**: their signatures follow the
call graph rather than any user-facing contract and may change without a major
version bump.

## File map (MATLAB → Python)

| MATLAB (`../src/`) | Python (`acdc/`) |
|---|---|
| `compute_gradient.m` + score expression | `objective.py` |
| `permutation_match.m` | `matching.py` |
| `do_frank_wolfe.m` | `frank_wolfe.py` |
| `evaluate_swaps.m`, `make_swaps.m`, `greedy_search.m` | `swaps.py` |
| `main_acdc.m`, `main_continuous.m`, `main_discrete.m` | `core.py` |

## Notes on fidelity & performance

- **`solver='dense'` (default)** uses SciPy's exact Hungarian/LAPJV solver
  (`linear_sum_assignment`) — the faithful analogue of the MATLAB
  `perfectMatching`. **`solver='sparse'`** uses
  `min_weight_full_bipartite_matching`, which is much faster and lighter at
  scale but only considers stored edges (a perfect matching must exist on the
  sparsity pattern).
- The MATLAB warm-start preconditioner for the assignment solver is omitted on
  purpose: its transform only adds per-row/column constants and does not change
  the optimal matching — it existed solely to speed MATLAB's sparse solver.
- Like the MATLAB original, the gradient and swap-gain matrices are dense
  `n × n`, so memory is `O(n²)`. At challenge scale (`n = 18524`, see below)
  that needs ~16 GB+ of RAM and is best run with `solver='dense'` on a machine
  with enough memory; small and medium graphs run comfortably anywhere.
- The quadratic part of the swap-gain matrix has two implementations. Six of
  its eight terms are `min(x, y)` with one argument taken from `A` or `P B P'`,
  so for **nonnegative** edge weights they vanish off those supports and can be
  scattered onto the union of the supports instead of swept densely. On sparse
  graphs that holds 2 dense `n × n` arrays instead of 5 (≈5.5 GB rather than
  ≈13.7 GB at challenge scale). Scattered writes are slower per element, so dense
  or negatively-weighted graphs use the literal dense expression; the two paths
  agree bit for bit and the choice is automatic.
- The greedy phase recomputes the full score after each trial swap (faithful to
  the MATLAB original); incremental scoring is a possible future optimization.
- Both loops stop early at a fixed point: a Frank–Wolfe update that leaves the
  iterate unchanged, or an AC/DC alternation that leaves the matching
  unchanged, makes every remaining iteration a deterministic replay. `max_iter`
  and `num_frank_wolfe` are therefore upper bounds, not exact counts. The
  result is unchanged — only the wasted work is skipped.

## TODO

### Pinned (fixed) matches

Neither the MATLAB original nor this translation supports **pinning** a subset of
matches so they are never swapped out. Every pair is free to move at every stage:
the discrete phase picks its swap from a `max` over the whole gain matrix
(`make_swaps.m` → `swaps.py`), and each Frank–Wolfe step solves an
*unconstrained* linear assignment problem over the full gradient
(`permutation_match.m` → `matching.py`) — the `P0` argument there is only a
warm-start accelerator and does not restrict the solution. This would be useful
whenever part of the correspondence is known a priori, e.g. neurons matched with
high confidence by cell type.

Adding it means two hooks, and would surface as a `pinned=` argument on the three
entry points:

- **Discrete phase** — mask the swap-gain matrix `D`: zero the row *and* column
  of each pinned index where `D` is assembled in `evaluate_swaps`, so that `dMax`
  also respects the pins and the `while dMax > 0` loop still terminates.
- **Continuous phase** — constrain the assignment problem. Prefer *reducing* it
  (drop the pinned rows/columns, solve on the free submatrix, reinsert the pinned
  pairs) over penalising forbidden entries with a large negative cost: the
  reduced LAP is strictly smaller and avoids `-inf` handling in
  `linear_sum_assignment` entirely. Only the two `permutation_match` calls need
  constraining — the interpolation `P + step*(Q - P)` with `step` in `(0, 1]`
  preserves any entry that is `1` in both `P` and `Q`, so pinned entries carry
  through the continuous iterate on their own.

## Benchmark

[`benchmark.py`](benchmark.py) runs AC/DC on synthetic *planted* instances of
increasing size to give a feel for scaling. Each instance relabels a random
graph by a hidden permutation, corrupts that alignment to make a warm start
(40 % of nodes scrambled by default), then measures runtime and how much of the
planted optimum is recovered.

```bash
python benchmark.py                                  # default sweep, ~1 min
python benchmark.py --sizes 100 200 400 800          # quick
python benchmark.py --sizes 1000 2000 --density 0.02 --solver sparse
```

Example run (laptop, `solver='dense'`, `max_iter=3`, `num_frank_wolfe=10`):

```
      n     edges       mem     grad     eval      acdc     cont     disc  start%  final%  nodes%
-------------------------------------------------------------------------------------------------
    100       481    78.1KB    0.001    0.002     0.014    0.010    0.004    38.9   100.0   100.0
    200      1935   312.5KB    0.003    0.003     0.021    0.013    0.008    36.0   100.0   100.0
    400      7771     1.2MB    0.007    0.009     0.060    0.038    0.022    39.3   100.0   100.0
    800     31147     4.9MB    0.028    0.035     0.219    0.139    0.080    38.1   100.0   100.0
   2000    194913    30.5MB    0.522    0.572     3.483    2.255    1.227    38.2   100.0   100.0
   5000   1219064   190.7MB    8.149    8.837    50.284   32.975   17.308    38.2   100.0   100.0
- - - - - - - - - - - - - - - - - - - - - - - -
  extrapolated from the n=800..5000 rows (acdc ~ n^2.97, R2=1.000 in log-log) -- NOT measured
 ~10000  4.88e+06   762.9MB     71.6     72.2       399      266      134       -       -       -
 ~50000  1.22e+08    18.6GB 1.04e+04 9.29e+03  4.73e+04 3.25e+04 1.52e+04       -       -       -
~100000  4.88e+08    74.5GB 8.89e+04 7.52e+04   3.7e+05 2.57e+05 1.16e+05       -       -       -
```

| extrapolated `n` | full `acdc_match` | one dense `n × n` array |
|---|---|---|
| 10 000 | ~6.7 min | 762.9 MB |
| 50 000 | ~13 h | 18.6 GB |
| 100 000 | ~4.3 d | 74.5 GB |

- `grad`/`eval` are per-call times (in seconds) for `compute_gradient` / `evaluate_swaps`
  (best of as many repeats as fit in a short budget, so small sizes are not
  dominated by noise). `acdc` is the full run, split into the continuous
  (`cont`) and discrete (`disc`) phases — roughly a 2:1 split throughout.
  `mem` is the size of one dense `n × n` float64 matrix; the `n = 5000` row
  peaks at about 1.1 GB resident.
- **Rows prefixed `~` are extrapolations, not measurements.** Each timing column
  is fitted as a power law over the asymptotic tail of the measured rows
  (`n ≥ 800`) and evaluated at the larger size. Fitting the *small* rows instead
  gives an exponent near 1.3 and under-predicts `n = 5000` by ~25×, because
  those sizes are dominated by fixed overhead — the benchmark warns when the
  fitted exponent lands below 2.5. Held out against a real measurement, a fit
  over `n = 400..2000` predicts 33 s at `n = 5000` against 42 s measured.
- **Scaling is `O(n³)` at fixed density**, not `O(n²)`: the gradient loops over
  the `n` matched nodes and does work quadratic in the degree, which at fixed
  density is itself proportional to `n`. The measured exponent is 2.97 with
  `R² = 1.000`. At fixed *average degree* — the realistic regime for a
  connectome — the same argument gives `O(n²)`, dominated by the dense
  matrices. Which of the two applies to your graph is the single biggest factor
  in what it will cost.
- The extrapolated rows hold density fixed at 0.05, so they are an **upper
  bound** for a graph of that size rather than a prediction for a real one. The
  challenge connectomes are `n = 18524` at density 0.012 (male) and 0.006
  (female) — 4–8× sparser than this sweep. At these sizes memory binds well
  before time does: 74.5 GB for a *single* array at `n = 100 000`, and the
  greedy phase holds two to three of them.
- From a corrupted-but-informative warm start the planted alignment is
  recovered essentially perfectly (`final% ≈ 100`). Quality from a *cold*
  (identity) start is much lower — AC/DC is a local search and, like the paper,
  relies on a reasonable warm start.

## The challenge data

The challenge was to align the connectomes of the **ventral nerve cords (VNCs)
of a male and a female fruit fly** — not the FlyWire brain connectome, despite
the challenge carrying the FlyWire name. Per the paper, each VNC connectome is
a directed weighted graph with `n = 18524` nodes (neurons) and millions of
edges (synapse counts); the organizers also supplied a cell-type-derived
baseline match, scoring 5 154 247, which teams could use as a warm start.

The two graphs are that size in the data shipped here as well: the male edge
list spans node ids `m1..m18524` and the female `f1..f18524` (the female list
has 18 523 nodes with at least one edge — `f9574` is isolated, which is why
`_coerce_pair` pads to a common `n`). Their densities are 0.012 and 0.006
respectively, so both are far sparser than the benchmark's default 0.05.

## CSV / challenge I/O

This package translates the **core algorithm only**. The challenge-specific CSV
readers/writers (`read_connectome.m`, `read_solution.m`, `save_solution.m`) and
the figure scripts are intentionally out of scope; build `A`, `B` as sparse
matrices (or pass edge lists) as shown above.

## Tests

```bash
pytest -q
```

Covers: a finite-difference gradient check, the score identity, self-match
optimality, recovery of a known permutation, exactness of the swap-gain matrix,
dense/sparse backend agreement, and end-to-end smoke tests.

## Attribution & licence

This package is a **translation** of the MATLAB reference implementation that
Lee, Matsliah & Saul published as supplementary material to their TMLR paper:

> D. D. Lee, A. Matsliah & L. K. Saul, *"AC/DC search: behind the winning
> solution to the FlyWire graph-matching challenge"*, Transactions on Machine
> Learning Research (01/2026).
> <https://openreview.net/forum?id=8MjCOMyaDf>

That original MATLAB code is **MIT-licensed**, Copyright (c) 2025 Daniel Lee and
Lawrence Saul. This Python translation is therefore also offered under the MIT
licence; [LICENSE](LICENSE) carries both notices, with the upstream licence
reproduced in full as it requires.

[NOTICE](NOTICE) summarises how the translation departs from the original —
0-based indexing, SciPy in place of MATLAB's internal `perfectMatching`, the
omitted assignment preconditioner, the untranslated challenge I/O, and the added
tests and benchmark. The "Notes on fidelity & performance" section above covers
the same ground in more detail.

If you use this package, please cite the paper above — the algorithm is the
authors' work; this is only a translation of it.
