feat: Population Annealing baseline + rebrand to PQQA

- Add qqa.population_annealing (Hukushima-Iba/Machta), with systematic and
  multinomial resampling, log-stable Boltzmann reweighting, ESS tracking,
  and the same SAResult-shaped contract as qqa.simulated_annealing.
- Refactor qqa.sa to share single-sweep MCMC primitives (QUBO Glauber and
  generic single-spin MH) with the new PA backend, eliminating the risk of
  ΔE-formula drift between SA and PA.
- Compare page: extend the shootout into a 3-way race PQQA vs SA vs PA;
  speedup vs the stronger of the two baselines.
- README: rename top-line title to "Parallel Quasi-Quantum Annealer (PQQA)"
  and reposition the project as an OSS toolkit unifying continuous-relaxation
  CO methods (PQQA / CRA-PI-GNN / CPRA) plus matched MCMC baselines.
- Tests: 8 new pytest cases (tests/test_pa.py), updated GUI AppTest for
  the renamed shootout mode, public-API regression covers PA exports.
- Verified: ruff clean, 206 passed / 2 skipped, verify_pa.py 7/7,
  verify_sa_compare.py 7/7.

Author: Yuma-Ichikawa

README.md

@@ -1,11 +1,17 @@
-# QQA4CO — Parallel Quasi-Quantum Annealing for Combinatorial Optimisation
+# Parallel Quasi-Quantum Annealer (PQQA)
 
-A research-grade PyTorch toolkit organised around **Parallel Quasi-Quantum
-Annealing (PQQA)** — a gradient-based, GPU-friendly framework for
-combinatorial optimisation (CO) that unifies several recent
-unsupervised-learning solvers under a single API. One installation gives you
-the PQQA solver, an optional GNN backend that plugs in CRA-PI-GNN-style
-methods, a 20+ class problem catalogue, a Streamlit dashboard and a CLI.
+A research-grade PyTorch toolkit for **continuous-relaxation combinatorial
+optimisation**. PQQA unifies the recent line of work that *lifts* a
+discrete CO problem into a continuous space and then solves it with
+gradient descent — **PQQA** (ICLR 2025), the **CRA-PI-GNN** baseline
+(NeurIPS 2024) and the **CPRA** diverse-solution framework (TMLR 2025) —
+behind a single, GPU-friendly API. One installation gives you the PQQA
+solver, an optional GNN backend that plugs in CRA-PI-GNN-style methods,
+GPU-parallel **SA** and **Population Annealing** baselines for honest
+comparisons, a 20+ class problem catalogue, a Streamlit dashboard and a CLI.
+
+The package is published on PyPI as **`qqa`** for backwards compatibility
+with earlier QQA4CO releases (``import qqa``).
 
 <p align="center">
   <a href="https://pypi.org/project/qqa/"><img src="https://img.shields.io/pypi/v/qqa.svg?logo=pypi&logoColor=white&label=PyPI" alt="PyPI version"></a>
@@ -51,7 +57,7 @@ parallel population in 3D solution-space PCA / loss-spectrogram views.
 | **Home** | Pick from 19 problems (MIS, Max-Cut, MaxClique, Vertex Cover, GraphBisection, MinDominatingSet, Coloring, BalancedGraphPartition, TSP, QAP, NQueens, Knapsack, NumberPartitioning, MaxSAT3, 1D Ising, Edwards–Anderson, SK, p-spin, RFIM, BinaryPerceptron, HopfieldMemory). Auxiliary sliders are problem-aware. |
 | **Solve** | Run PQQA / CRA-PI-GNN / CPRA with live progress, polish (1-flip local search) and warm-start toggles, and a per-problem solution viewer (TSP tour, NQueens board, highlighted IS, colouring, …). |
 | **Visualize** | 10 tabs of post-hoc plots — **solution-space PCA** of the final parallel population (3D, replicas coloured by loss, global best highlighted), **loss spectrogram** over time, diversity, replica fate, schedule, ridgeline. |
-| **Compare** | Hyper-parameter grid sweep with parallel-coordinates view, **and** a head-to-head **PQQA vs. SA shootout** that reports the wall-clock speed-up at matched compute budget. |
+| **Compare** | Hyper-parameter grid sweep with parallel-coordinates view, **and** a head-to-head **PQQA vs. SA vs. Population-Annealing shootout** that reports the wall-clock speed-up at matched compute budget. |
 
 ### Or run it locally — same UI, your hardware
 
@@ -92,12 +98,20 @@ CUDA is picked up automatically when available.
    diverse-solution framework (TMLR 2025) all share `AnnealResult`,
    `score_summary`, the same problem builders and the same CLI flags — A/B
    comparing methods is one `--backend` switch away.
-4. **A polished Streamlit dashboard** (light / dark, live progress, parallel
+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**
+   (`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.
+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
    experiments. A hosted instance lives at
    <https://parallelquasiquantum4co.streamlit.app/>.
-5. **MkDocs + Material reference docs** with auto-generated API pages, a
+6. **MkDocs + Material reference docs** with auto-generated API pages, a
    runnable `examples/` notebook gallery (Open-in-Colab badges) and a
    `scripts/verify_all_problems.py` correctness sweep that benchmarks every
    problem against ground truth or a strong baseline (29 / 29 instances pass

app/pages/3_Compare.py

@@ -69,12 +69,12 @@ def _scalar_best(obj) -> float:
 with st.sidebar:
     mode = st.radio(
         "Compare mode",
-        ["QQA hyper-parameter sweep", "PQQA vs SA shootout"],
+        ["QQA hyper-parameter sweep", "PQQA vs SA vs PA shootout"],
         index=0,
         help=(
             "Sweep: grid-search QQA hyper-parameters on the current problem. "
-            "Shootout: race PQQA against the SA baseline at a matched compute "
-            "budget — useful for showing how much the relaxation buys you."
+            "Shootout: race PQQA against the SA and Population-Annealing "
+            "baselines at a matched compute budget."
         ),
     )
 
@@ -151,25 +151,27 @@ def _scalar_best(obj) -> float:
 
 
 # =============================================================================
-# Mode 2 — PQQA vs SA shootout (new)
+# Mode 2 — PQQA vs SA vs PA shootout
 # =============================================================================
-elif mode == "PQQA vs SA shootout":
-    # SA only supports binary / spin relaxations; refuse upfront on Categorical.
+elif mode == "PQQA vs SA vs PA shootout":
+    # SA / PA only support binary / spin relaxations; refuse on Categorical.
     probe_problem = build_problem(cfg)
-    sa_supported = not isinstance(getattr(probe_problem, "relaxation", None), CategoricalRelaxation)
+    chain_supported = not isinstance(
+        getattr(probe_problem, "relaxation", None), CategoricalRelaxation
+    )
 
     st.markdown(
-        "Race **Parallel QQA** against a textbook **Simulated Annealing** "
-        "baseline on the same problem. Both run on the device you picked on "
-        "Home, with a comparable compute budget — the convergence plot below "
-        "exposes the wall-clock gap directly."
+        "Race **Parallel QQA** against two textbook MCMC baselines — "
+        "**Simulated Annealing (SA)** and **Population Annealing (PA, with "
+        "resampling)** — on the same problem. Total compute is matched: "
+        "PA's ``num_temps × sweeps_per_temp`` ≈ SA's ``num_sweeps``."
     )
 
-    if not sa_supported:
+    if not chain_supported:
         st.warning(
-            "SA is not supported for this problem (categorical relaxation). "
-            "Switch the problem on Home to a QUBO / Ising / spin family "
-            "(e.g. MIS, Max-Cut, SK) to use the shootout."
+            "SA / PA do not support categorical-relaxation problems. "
+            "Pick a QUBO / Ising / spin family (e.g. MIS, Max-Cut, SK) on "
+            "**Home** to use the shootout."
         )
 
     with st.sidebar:
@@ -187,9 +189,17 @@ def _scalar_best(obj) -> float:
         sa_beta_end = st.slider("SA β_end", 1.0, 100.0, 10.0, step=1.0, key="sa_b1")
         sa_schedule = st.selectbox("SA schedule", ["geometric", "linear"], key="sa_sched")
 
-        run_shootout = st.button("▶ Run shootout", type="primary", disabled=not sa_supported)
+        st.caption("Population Annealing baseline")
+        pa_temps = st.slider("PA num_temps", 10, 500, 100, step=10, key="pa_temps")
+        pa_sweeps = st.slider("PA sweeps_per_temp", 1, 50, 10, key="pa_sweeps")
+        pa_chains = st.slider("PA population (sol_size)", 4, 256, 64, key="pa_chains")
+        pa_beta_start = st.slider("PA β_start", 0.01, 1.0, 0.1, step=0.01, key="pa_b0")
+        pa_beta_end = st.slider("PA β_end", 1.0, 100.0, 10.0, step=1.0, key="pa_b1")
+        pa_resample = st.selectbox("PA resample", ["systematic", "multinomial"], key="pa_resample")
+
+        run_shootout = st.button("▶ Run shootout", type="primary", disabled=not chain_supported)
 
-    if sa_supported and run_shootout:
+    if chain_supported and run_shootout:
         with st.spinner("Running PQQA…"):
             problem_pqqa = build_problem(cfg)
             res_pqqa = qqa.anneal(
@@ -215,11 +225,25 @@ def _scalar_best(obj) -> float:
                 verbose=False,
             )
 
+        with st.spinner("Running Population Annealing baseline…"):
+            problem_pa = build_problem(cfg)
+            res_pa = qqa.population_annealing(
+                problem_pa,
+                sol_size=int(pa_chains),
+                num_temps=int(pa_temps),
+                sweeps_per_temp=int(pa_sweeps),
+                beta_start=float(pa_beta_start),
+                beta_end=float(pa_beta_end),
+                resample=pa_resample,
+                device=cfg["device"],
+                verbose=False,
+            )
+
         pqqa_best = _scalar_best(res_pqqa.best_obj)
         sa_best = _scalar_best(res_sa.best_obj)
+        pa_best = _scalar_best(res_pa.best_obj)
 
-        # Side-by-side summary cards.
-        col_p, col_s = st.columns(2)
+        col_p, col_s, col_pa = st.columns(3)
         col_p.metric(
             "PQQA best_obj",
             f"{pqqa_best:.4f}",
@@ -232,60 +256,66 @@ def _scalar_best(obj) -> float:
             delta=f"runtime {res_sa.runtime:.2f} s",
             delta_color="off",
         )
+        col_pa.metric(
+            "PA best_obj",
+            f"{pa_best:.4f}",
+            delta=f"runtime {res_pa.runtime:.2f} s",
+            delta_color="off",
+        )
 
-        # "How much faster did PQQA reach SA's best?" by walking the SA history.
-        sa_history = res_sa.history.get("best_obj", []) or [sa_best]
-        # Time per SA sweep (rough), so we can map sweep index -> wall-clock.
-        sa_dt = float(res_sa.runtime) / max(1, len(sa_history))
-        # Find the first sweep at which SA reached the PQQA best.
-        sa_reach_idx = next(
-            (i for i, v in enumerate(sa_history) if float(v) <= pqqa_best),
+        # "How much faster did PQQA reach the best baseline's best?"
+        baseline_best = min(sa_best, pa_best)
+        baseline_label = "SA" if sa_best <= pa_best else "PA"
+        baseline_res = res_sa if sa_best <= pa_best else res_pa
+        baseline_history = baseline_res.history.get("best_obj", []) or [baseline_best]
+        baseline_dt = float(baseline_res.runtime) / max(1, len(baseline_history))
+        baseline_reach = next(
+            (i for i, v in enumerate(baseline_history) if float(v) <= pqqa_best),
             None,
         )
 
-        if sa_reach_idx is None:
+        if baseline_reach is None:
             st.warning(
-                f"SA never matched PQQA's best ({pqqa_best:.4f}). "
-                f"SA stalled at {sa_best:.4f} after {res_sa.runtime:.2f}s. "
-                "Increase `SA num_sweeps` or `β_end` to give SA more room."
+                f"Neither SA nor PA matched PQQA's best ({pqqa_best:.4f}). "
+                f"Best baseline = {baseline_label} at {baseline_best:.4f} "
+                "after the chosen budget. Increase the baseline's compute or "
+                "tune β to give it more room."
             )
         else:
-            sa_time_to_pqqa = sa_dt * (sa_reach_idx + 1)
-            speedup = sa_time_to_pqqa / max(res_pqqa.runtime, 1e-6)
+            baseline_t = baseline_dt * (baseline_reach + 1)
+            speedup = baseline_t / max(res_pqqa.runtime, 1e-6)
             st.success(
                 f"PQQA reached `best_obj = {pqqa_best:.4f}` in "
-                f"**{res_pqqa.runtime:.2f}s**, while SA needed "
-                f"~**{sa_time_to_pqqa:.2f}s** to match it "
-                f"(speedup ≈ **{speedup:.1f}×**)."
+                f"**{res_pqqa.runtime:.2f}s**, while {baseline_label} "
+                f"(the stronger baseline) needed ~**{baseline_t:.2f}s** to "
+                f"match it (speedup ≈ **{speedup:.1f}×**)."
             )
 
-        # Convergence overlay: x = epoch / sweep index, y = best_obj.
         try:
             import plotly.graph_objects as go
 
             fig = go.Figure()
-            pqqa_hist = res_pqqa.history.get("best_obj", []) or [pqqa_best]
-            fig.add_trace(
-                go.Scatter(
-                    x=list(range(len(pqqa_hist))),
-                    y=[float(v) for v in pqqa_hist],
-                    name="PQQA",
-                    mode="lines",
-                    line={"width": 3},
-                )
-            )
-            fig.add_trace(
-                go.Scatter(
-                    x=list(range(len(sa_history))),
-                    y=[float(v) for v in sa_history],
-                    name="SA",
-                    mode="lines",
-                    line={"width": 3, "dash": "dash"},
+            for name, hist, fallback, dash in (
+                ("PQQA", res_pqqa.history.get("best_obj", []), pqqa_best, None),
+                ("SA", res_sa.history.get("best_obj", []), sa_best, "dash"),
+                ("PA", res_pa.history.get("best_obj", []), pa_best, "dot"),
+            ):
+                hist = hist or [fallback]
+                line = {"width": 3}
+                if dash:
+                    line["dash"] = dash
+                fig.add_trace(
+                    go.Scatter(
+                        x=list(range(len(hist))),
+                        y=[float(v) for v in hist],
+                        name=name,
+                        mode="lines",
+                        line=line,
+                    )
                 )
-            )
             fig.update_layout(
                 title="Best-objective trajectory",
-                xaxis_title="iteration (PQQA epoch / SA sweep)",
+                xaxis_title="iteration (PQQA epoch / SA sweep / PA temperature step)",
                 yaxis_title="best_obj (lower is better)",
                 **plotly_layout(),
             )
@@ -294,7 +324,6 @@ def _scalar_best(obj) -> float:
         except Exception as e:
             st.info(f"Plot unavailable: {e}")
 
-        # Final results table for export.
         st.subheader("Summary table")
         st.dataframe(
             [
@@ -312,6 +341,13 @@ def _scalar_best(obj) -> float:
                     "iterations": int(sa_sweeps),
                     "sol_size": int(sa_chains),
                 },
+                {
+                    "backend": "PA",
+                    "best_obj": pa_best,
+                    "runtime_s": float(res_pa.runtime),
+                    "iterations": int(pa_temps) * int(pa_sweeps),
+                    "sol_size": int(pa_chains),
+                },
             ],
             width="stretch",
         )

src/qqa/__init__.py

@@ -38,6 +38,7 @@
     PopulationTracker,
     TrajectoryTracker,
 )
+from qqa.pa import PAResult, population_annealing
 from qqa.problems import (
     QAP,
     TSP,
@@ -122,6 +123,7 @@
     "MinimumDominatingSet",
     "NQueens",
     "NumberPartitioning",
+    "PAResult",
     "PSpinGlass",
     "PopulationTracker",
     "QUBOProblem",
@@ -141,6 +143,7 @@
     "generate_graph",
     "load_problem_from_file",
     "polish",
+    "population_annealing",
     "simulated_annealing",
     "user_problem_from_source",
     "warmstart",