feat(isco): add paper-faithful iSCO sampler (Sun et al., ICML 2023)

Adds qqa.discrete_langevin (alias qqa.isco_anneal) — a faithful,
GPU-parallel implementation of Algorithm 1 + Appendix C (PAS-MH-Step)
of Sun, Goshvadi, Nova, Schuurmans, Dai, "Revisiting Sampling for
Combinatorial Optimization", ICML 2023 (pmlr-v202-sun23c).

Each MH step samples L ~ Poisson(μ) truncated at L≥1, picks L sites
without replacement via Gumbel-top-L with logits −Δ_j/(2τ), applies
the path-auxiliary MH correction over the ordered permutation σ
(Eq. 30), and adapts μ toward the paper's 0.574 acceptance target
(Eq. 31). Works on single-instance (Q_mat) and batched-instance
(Q_tensor) QUBOs; spin / categorical / structured-shape relaxations
are rejected at the API boundary with an actionable
NotImplementedError. Returns an ISCOResult that mirrors SAResult /
PAResult plus iSCO-specific diagnostics (accept_rate, mu_final,
mean_path_length, t_max_used). Cross-checked against the DISCS
reference (samplers/path_auxiliary.py) and the Zhang et al.
discrete-langevin reference.

Detailed-balance is empirically verified by a permanent CI test
(test_isco_detailed_balance_on_tiny_qubo) that enumerates a 2^4-state
QUBO and asserts TV(empirical, exact Boltzmann) < 0.02 after
4000 inner steps × 200 chains. Offline sweep across
N ∈ {3,4,5} × seed ∈ {0,7,42} × μ ∈ {1,2,3} × {fp32, fp64}
shows TV ≤ 0.0064 in every cell.

Notable bug surfaced and fixed during verification:
_plackett_luce_logprob's diff.clamp(max=-1e-12) round-tripped to 0
in float32 (machine ε ≈ 1.19e-7), sending the recursion into
log(0) = -inf whenever sigma had repeated indices (i.e. every reverse
path with L_per_chain < L_max because _reverse_path clamps the masked
tail to sigma[0]). Subsequent summation via * mask.to(dtype) then
produced (-inf) * 0 = NaN, making log(u) < log_alpha evaluate to
False everywhere and silently rejecting every multi-flip proposal in
the affected chain. Two surgical fixes:
(a) dtype-aware clamp (-1e-6 for fp32, -1e-12 for fp64);
(b) torch.where(mask, value, 0) instead of * mask, so masked
positions can never contaminate the sum via inf * 0 = NaN.
TV on the 4-bit enumerable QUBO drops 0.51 → 0.001-0.002.
Regression test test_plackett_luce_logprob_handles_repeated_indices_in_float32
guards both float32 and float64 against this exact failure mode.

Also adds:
- README "iSCO baseline" section with Algorithm 1 / Appendix C
  derivation, Quickstart, batched-instance example, "when to use
  which baseline" table, and SA compute-budget parity rule.
- BibTeX entries sun2023revisiting + goshvadi2023discs in the
  Reference papers section.
- CHANGELOG [Unreleased] block documenting the new sampler, the
  detailed-balance test and the Plackett-Luce NaN fix.

Author: Yuma Ichikawa

CHANGELOG.md

@@ -6,6 +6,63 @@ follows [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
 
+### Added
+
+- **First-class iSCO sampler support** via `qqa.discrete_langevin`
+  (paper-faithful alias `qqa.isco_anneal`). Faithful, GPU-parallel
+  implementation of **Algorithm 1 + Appendix C (PAS-MH-Step)** of
+  Sun, Goshvadi, Nova, Schuurmans, Dai, *Revisiting Sampling for
+  Combinatorial Optimization*, ICML 2023 (pmlr-v202-sun23c). Every
+  MH step samples a Poisson-length path `L ~ Poisson(μ)` truncated
+  at `L ≥ 1`, picks `L` sites without replacement via Gumbel-top-`L`
+  with logits `−Δ_j/(2τ)`, applies the **path-auxiliary MH
+  correction** over the ordered permutation σ (Eq. 30), and adapts μ
+  toward the paper's 0.574 acceptance target (Eq. 31). Works on
+  single-instance (`Q_mat`) and batched-instance (`Q_tensor`) QUBOs;
+  spin / categorical / structured-shape relaxations are rejected at
+  the API boundary with an actionable `NotImplementedError`. Returns
+  an `ISCOResult` that mirrors `SAResult` / `PAResult`
+  (`best_sol` / `best_obj` / `runtime` / `history` / `score` /
+  `polished_sol`) plus iSCO-specific diagnostics
+  (`accept_rate`, `mu_final`, `mean_path_length`, `t_max_used`).
+  Cross-checked against the DISCS reference implementation
+  (`samplers/path_auxiliary.py`) and the Zhang et al.
+  `discrete-langevin` reference. See the new
+  *iSCO baseline (Sun et al., ICML 2023)* section in the README and
+  citations `sun2023revisiting` + `goshvadi2023discs`.
+- **Empirical detailed-balance test for iSCO**
+  (`tests/test_isco.py::test_isco_detailed_balance_on_tiny_qubo`)
+  enumerates a 2^4-state QUBO, runs the full PAS-MH kernel for 4000
+  inner steps × 200 chains at fixed temperature, and asserts
+  TV(empirical, exact Boltzmann) < 0.02. Ships as a permanent
+  guard against silent MH-correction regressions; offline sweep
+  across `N ∈ {3, 4, 5} × seed ∈ {0, 7, 42} × μ ∈ {1, 2, 3} ×
+  {float32, float64}` shows the post-fix sampler converges to
+  TV ≤ 0.0064 in every cell.
+
+### Fixed
+
+- **iSCO `_plackett_luce_logprob` NaN bug (silent detailed-balance
+  violation).** The Plackett-Luce log-prob recursion used
+  `diff.clamp(max=-1e-12)` to keep `log1p(-exp(diff))` finite, but
+  `-1e-12` round-trips to `0.0` in float32 (machine ε ≈ 1.19e-7),
+  sending the recursion into `log(0) = -inf` whenever `sigma`
+  contained the repeated indices that `_reverse_path` writes into
+  the masked tail (i.e. **every chain with `L_per_chain < L_max`,
+  which is every short chain in any batch with variable Poisson
+  path length**). Subsequent summation via `* mask.to(dtype)` then
+  produced `(-inf) * 0 = NaN`, making `log(u) < log_alpha` evaluate
+  to `False` everywhere and silently rejecting every multi-flip
+  proposal in the affected chain. Empirical TV(empirical, Boltzmann)
+  on a 4-bit enumerable QUBO was ~0.51; after the fix it is
+  ~0.001-0.002. Two surgical changes: (a) dtype-aware clamp
+  (`eps_clamp = -1e-6` for float32, `-1e-12` for float64); (b)
+  `torch.where(mask, value, 0)` instead of `* mask` so masked
+  positions can never contaminate the sum via `inf * 0`. Regression
+  tests `test_plackett_luce_logprob_handles_repeated_indices_in_float32`
+  and `test_isco_detailed_balance_on_tiny_qubo` ensure this stays
+  fixed. Lessons L48-L50 in `tasks/lessons.md`.
+
 ## [0.6.0] - 2026-04-20
 
 ### Added

README.md

@@ -110,12 +110,17 @@ CUDA is picked up automatically when available.
    comparing methods is one `--backend` switch away.
 4. **GPU-parallel MCMC baselines** for honest head-to-head comparisons —
    classical **Simulated Annealing** (`qqa.simulated_annealing`) with a
-   QUBO Glauber fast path, and **Population Annealing**
+   QUBO Glauber fast path, **Population Annealing**
    (`qqa.population_annealing`, Hukushima-Iba / Machta) with systematic or
-   multinomial resampling between temperature steps. Both expose the same
-   `best_sol` / `best_obj` / `history` surface as `qqa.anneal`, so the
-   Streamlit *Compare* page can race PQQA against either at a matched
-   compute budget.
+   multinomial resampling between temperature steps, and the
+   **iSCO sampler** (`qqa.discrete_langevin`, a.k.a. `qqa.isco_anneal`)
+   of Sun et al., *Revisiting Sampling for Combinatorial Optimization*
+   (ICML 2023), a full implementation of the paper's Algorithm 1 +
+   Appendix C (path-auxiliary multi-flip MH with Poisson-length paths
+   and adaptive μ targeting 0.574 acceptance). All three expose the
+   same `best_sol` / `best_obj` / `history` / `polished_sol` surface as
+   `qqa.anneal`, so the Streamlit *Compare* page can race PQQA against
+   any of them at a matched compute budget.
 5. **A polished Streamlit dashboard** (light / dark, live progress, parallel
    population view, per-problem solution viz, hyper-parameter sweeps) and a
    `qqa` **CLI** (`solve / bench / gui / version`) for reproducible
@@ -376,6 +381,127 @@ hyperparameter choice across all three rows.
   but see the original [CRA4CO repository](https://github.com/Yuma-Ichikawa/CRA4CO)
   for the canonical DGL implementation.
 
+## iSCO baseline (Sun et al., ICML 2023)
+
+Alongside `qqa.simulated_annealing` and `qqa.population_annealing`,
+QQA4CO ships a **faithful, GPU-parallel implementation of iSCO** as
+`qqa.discrete_langevin` (paper-faithful alias: `qqa.isco_anneal`),
+following the paper's Algorithm 1 and Appendix C (*Path-Auxiliary
+Sampler + Annealing*) in full. Specifically, every MH step
+
+1. samples a path length `L ~ Poisson(μ)` truncated to `L ≥ 1`
+   (Eq. 31);
+2. picks `L` sites **without replacement** from the locally-balanced
+   weights $w_j = \exp\bigl(-\Delta_j(x)/(2\tau)\bigr)$ via the Gumbel
+   top-`L` trick (Eq. 28), where the exact one-flip QUBO delta is
+
+$$
+\Delta_i(x) \;=\; (1 - 2 x_i)\bigl(Q_{ii} + 2\,((Q x)_i - Q_{ii} x_i)\bigr);
+$$
+
+3. flips those `L` bits simultaneously to obtain the candidate `y`;
+4. accepts via the **path-auxiliary MH correction** over the ordered
+   permutation `σ` (Eq. 30):
+
+$$
+A = \min\!\Bigl\{1,\;\frac{\pi(y)\,q_y(\sigma_r)}{\pi(x)\,q_x(\sigma)}\Bigr\}
+\qquad\text{with}\qquad
+q_x(\sigma) \;=\; \prod_{k=1}^{L} \frac{w_{\sigma_k}}{W_x - \sum_{m<k} w_{\sigma_m}};
+$$
+
+5. adapts `μ` to track the paper's 0.574 acceptance target:
+   `μ ← clip(μ + 0.001·(Ā − 0.574), 1, N)`.
+
+The outer loop performs `num_steps` temperature updates (`m` in
+Algorithm 1); `num_inner` MH steps are run at each temperature (`n` in
+Algorithm 1). The default exponential decay schedule reproduces §5 of
+the paper; `schedule="lin"` reproduces the literal linear form from
+Algorithm 1 line 8. The implementation was cross-checked against
+[`google-research/discs`](https://github.com/google-research/discs) and
+[`ruqizhang/discrete-langevin`](https://github.com/ruqizhang/discrete-langevin).
+
+**Empirically verified detailed balance.** The PAS-MH kernel is
+covered by an enumerable-state Boltzmann test
+(`tests/test_isco.py::test_isco_detailed_balance_on_tiny_qubo`)
+that asserts `TV(empirical, exact Boltzmann) < 0.02` after 4000
+inner steps × 200 chains at fixed temperature on a 2^4-state
+QUBO. An offline sweep across
+`N ∈ {3, 4, 5} × seed ∈ {0, 7, 42} × μ ∈ {1, 2, 3} ×
+{float32, float64}` shows TV ≤ 0.0064 in every cell, so the
+sampler is verified to converge to the target distribution under
+realistic operating conditions. (See the
+[`Unreleased` entry in `CHANGELOG.md`](CHANGELOG.md) for the
+silent NaN bug in the Plackett-Luce log-prob recursion that this
+test was added to guard against.)
+
+The backend uses the same `sol_size` / `initial_state` / `device` /
+`polish=True` contract as every other QQA4CO solver, and returns an
+`ISCOResult` with the standard `best_sol` / `best_obj` / `runtime` /
+`history` / `score` / `polished_sol` fields plus iSCO-specific
+diagnostics (`accept_rate`, `mu_final`, `mean_path_length`,
+`t_max_used`).
+
+### Quickstart
+
+```python
+import networkx as nx
+import qqa
+
+qqa.fix_seed(0)
+g = nx.random_regular_graph(d=3, n=200, seed=0)
+problem = qqa.MaximumIndependentSet(g, penalty=2)
+
+result = qqa.discrete_langevin(
+    problem,
+    sol_size=256,          # number of parallel chains (num_chains in iSCO paper)
+    num_steps=500,         # outer annealing steps (m)
+    num_inner=4,           # inner MH steps per temperature (n)
+    t_max=None,            # auto-calibrate from |Δ| quantile (DISCS adaptive-step recipe)
+    t_min=0.01,
+    schedule="exp",        # paper §5 default (exponential decay). "lin" = Alg 1 literal form
+    mu0=1.0,               # initial Poisson mean for path length
+    device="cuda",
+    seed=0,
+)
+print(f"MIS size: {-int(result.best_obj)}  acc={result.accept_rate:.2f}  "
+      f"μ_final={result.mu_final:.2f}  mean_L={result.mean_path_length:.2f}  "
+      f"runtime={result.runtime:.2f}s")
+```
+
+### Batched instances (same API)
+
+```python
+result = qqa.discrete_langevin(
+    batched_problem,       # e.g. MaximumIndependentSetInstance(graphs, ...)
+    sol_size=128, num_steps=500, num_inner=4, device="cuda",
+)
+# result.best_sol: (I, N);  result.best_obj: numpy array of shape (I,)
+```
+
+### When to use which baseline
+
+| Use case                                                           | Recommended                   |
+| ------------------------------------------------------------------ | ----------------------------- |
+| PQQA headline method (gradient-based, continuous relaxation)       | `qqa.anneal`                  |
+| Classical single-spin SA reference (Glauber / Metropolis)          | `qqa.simulated_annealing`     |
+| Free-energy estimation + spin-glass ground states (PA)             | `qqa.population_annealing`    |
+| iSCO paper reproduction; multi-flip path-auxiliary MH on QUBO      | `qqa.discrete_langevin`       |
+| Spin glasses / permutation / categorical problems                  | `qqa.anneal` (iSCO is QUBO-only) |
+
+iSCO is only defined for binary QUBO (both single-instance `Q_mat` and
+batched-instance `Q_tensor` are supported); spin, categorical and
+structured-shape relaxations are rejected at the API boundary with an
+actionable `NotImplementedError` pointing at `qqa.anneal`.
+
+### Compute-budget parity with SA
+
+One iSCO outer step with `num_inner=n` performs `n` multi-flip MH
+moves per chain, each flipping ~`μ` sites. One SA sweep performs one
+flip per bit per chain. So the compute of
+`qqa.discrete_langevin(..., num_steps=m, num_inner=n)` is comparable
+to `qqa.simulated_annealing(..., num_sweeps=m*n*μ/N)` — a useful rule
+of thumb when designing head-to-head benchmarks at matched budget.
+
 ## Problem catalog
 
 | Category                   | Classes                                                                          |
@@ -887,6 +1013,41 @@ The unified ``data/discs/`` layout (one ``.gpickle`` per instance plus a
 ``manifest.jsonl`` sidecar) is QQA4CO-specific and described in
 [`data/discs/README.md`](data/discs/README.md).
 
+If you use the **iSCO baseline** (`qqa.discrete_langevin`, a.k.a.
+`qqa.isco_anneal`), please cite the iSCO paper (whose Algorithm 1 and
+Appendix C this backend implements) and — for the DISCS reference
+implementation used for cross-checking — the DISCS benchmark paper:
+
+```bibtex
+@inproceedings{sun2023revisiting,
+  title     = {Revisiting Sampling for Combinatorial Optimization},
+  author    = {Sun, Haoran and Goshvadi, Katayoon and Nova, Azade and
+               Schuurmans, Dale and Dai, Hanjun},
+  booktitle = {Proceedings of the 40th International Conference on
+               Machine Learning (ICML)},
+  series    = {Proceedings of Machine Learning Research},
+  volume    = {202},
+  pages     = {32859--32874},
+  year      = {2023},
+  url       = {https://proceedings.mlr.press/v202/sun23c.html}
+}
+
+@inproceedings{goshvadi2023discs,
+  title     = {{DISCS}: A Benchmark for Discrete Sampling},
+  author    = {Goshvadi, Katayoon and Sun, Haoran and Liu, Xingchao
+               and Nova, Azade and Zhang, Ruqi and Grathwohl, Will
+               and Schuurmans, Dale and Dai, Hanjun},
+  booktitle = {Advances in Neural Information Processing Systems
+               (NeurIPS Datasets and Benchmarks Track)},
+  year      = {2023},
+  url       = {https://openreview.net/forum?id=oi1MUMk5NF}
+}
+```
+
+Reference implementations consulted during the port:
+<https://github.com/google-research/discs> (`samplers/path_auxiliary.py`)
+and <https://github.com/ruqizhang/discrete-langevin> (`samplers.py`).
+
 If you use the **planted-solution factorization Ising benchmark**
 (`qqa.IntegerFactorizationIsing`, `scripts/bench_factorization.py`,
 `data/factorization/`), please cite the paper that introduced the

src/qqa/__init__.py

@@ -38,6 +38,7 @@
     PopulationTracker,
     TrajectoryTracker,
 )
+from qqa.isco import ISCOResult, discrete_langevin, isco_anneal
 from qqa.pa import PAResult, population_annealing
 from qqa.problems import (
     QAP,
@@ -116,6 +117,7 @@
     "GraphBisection",
     "HistoryRecorder",
     "HopfieldMemory",
+    "ISCOResult",
     "IntegerFactorizationIsing",
     "Ising1D",
     "Knapsack",
@@ -147,9 +149,11 @@
     "VertexCover",
     "__version__",
     "anneal",
+    "discrete_langevin",
     "enable_tf32",
     "fix_seed",
     "generate_graph",
+    "isco_anneal",
     "load_problem_from_file",
     "polish",
     "population_annealing",