Metadata-Version: 2.4
Name: acdc-search
Version: 0.2.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: numba>=0.59
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 needs NumPy,
SciPy and numba.

For a development install from a checkout:

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

**numba.** The three inner loops that dominate the runtime are compiled
(`acdc/_kernels.py`); see [Notes on fidelity & performance](#notes-on-fidelity--performance).
They are cached to disk, so the first call in a fresh install pays about a
second of compilation and later runs start in ~0.1 s. If the package is
installed somewhere unwritable, numba silently recompiles once per process
instead — set `NUMBA_CACHE_DIR` to a writable path to avoid that.

## 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`, `gradient_entries`),
`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. `acdc._kernels` holds the compiled inner
loops and is private.

## 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` |
| — (compiled inner loops, no MATLAB counterpart) | `_kernels.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.
- **Three inner loops are compiled with numba** (`acdc/_kernels.py`): the
  gradient accumulation, the greedy pairwise-swap loop, and the
  support-restricted part of the swap-gain quadratic term. Each walks sparse
  adjacency structure one entry at a time over only a few dozen elements, so in
  NumPy they were dominated by per-call dispatch overhead rather than
  arithmetic. The kernels visit elements in the same order as the reference
  implementation, so results are **bit-identical** to the pure-NumPy version
  they replace, not merely close.
- 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 are
  applied cell by cell over 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). 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 scores each trial swap incrementally.** Exchanging the
  images of `i` and `j` can only change cells in rows and columns `i` and `j`,
  so the exact score change follows from `O(deg)` work where `make_swaps.m`
  recomputes the whole `O(nnz)` score after every trial. The accept/reject
  decision is identical — only the cost of reaching it differs — and the score
  is re-derived in full once per greedy iteration, so nothing accumulates.
  This was the single largest saving: on one greedy pass at `n = 1200` it took
  148.6 s down to 0.65 s.
- **The Frank–Wolfe line search never materializes the vertex gradient.** It
  reads `∇J(Q)` only on the support of `Q − P`, which has `O(n)` nonzeros, so
  those entries are evaluated pointwise instead of building a dense `n × n`
  array and discarding all but `O(n)` of it.
- `A` and `B` are needed in both CSR and CSC form by every gradient and every
  swap evaluation, but are fixed for a whole run; both orientations are built
  once and cached on the matrix rather than re-derived per call.
- 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.000    0.000     0.004    0.003    0.001    38.9   100.0   100.0
    200      1935   312.5KB    0.000    0.000     0.003    0.001    0.001    36.0   100.0   100.0
    400      7771     1.2MB    0.000    0.002     0.008    0.003    0.004    39.3   100.0   100.0
    800     31147     4.9MB    0.002    0.008     0.034    0.014    0.019    38.1   100.0   100.0
   2000    194913    30.5MB    0.108    0.157     0.695    0.371    0.323    38.2   100.0   100.0
   5000   1219064   190.7MB    2.722    3.094    14.461    8.295    6.165    38.2   100.0   100.0
- - - - - - - - - - - - - - - - - - - - - - - -
  extrapolated from the n=800..5000 rows (acdc ~ n^3.30, R2=1.000 in log-log) -- NOT measured
 ~10000  4.88e+06   762.9MB     45.5       29       142     94.5     53.4       -       -       -
 ~50000  1.22e+08    18.6GB  2.4e+04 5.31e+03  2.89e+04 2.53e+04 8.42e+03       -       -       -
~100000  4.88e+08    74.5GB 3.58e+05 5.01e+04  2.85e+05 2.81e+05 7.45e+04       -       -       -
```

| extrapolated `n` | full `acdc_match` | one dense `n × n` array |
|---|---|---|
| 10 000 | ~2.4 min | 762.9 MB |
| 50 000 | ~8 h | 18.6 GB |
| 100 000 | ~3.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 3.30 with
  `R² = 1.000`. (It reads higher than the `2.97` measured before the inner loops
  were compiled, because removing a large *constant* per-element overhead makes
  the small rows — which are dominated by fixed costs — relatively cheaper, and
  steepens the fitted line. The asymptotic behaviour is unchanged.) 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. Two tests pin the
compiled kernels to the definitions they stand in for: `gradient_entries`
against the dense gradient, and `make_swaps` against a literal
recompute-the-whole-score-per-trial reference.

## 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.
