Parallel tempering rewinds the iteration counter after a replica exchange, so iteration exchange_every+1 is processed twice and can write two samples.txt rows for one iteration

[wshlavacek]

What happens

replica_exchange calls try_to_choose_new_pset(j) to restart each chain, then does self.iteration[j] -= 1 with the comment that the counter "got off by 1 because try_to_choose_new_pset() was called twice". But try_to_choose_new_pset is not a pure counter bump: after self.iteration[index] += 1 it runs the sampling/histogram/output/diagnostics block (lines 221-243) — if self.iteration[index] > self.burn_in and self.iteration[index] % self.sample_every == 0 and self.should_sample(index): self.sample_pset(...), self.update_histograms(...), self.output_results(), and report_convergence_diagnostics. So the call inside replica_exchange executes all of that at iteration E+1, the decrement rewinds to E, and the next got_result increments back to E+1 and executes the whole block a second time. The rewind was written as if the second call were a no-op; it is not. The two sample rows are not even the same state (the accept/reject between them may have moved the chain), so one iteration of the chain contributes two draws to samples.txt, the histograms and the credible intervals.

Reproduction

Executed (read-only, in scratchpad): /private/tmp/claude-503/-Users-l119605-Code-PyBNF/eb55c5a3-1006-4c5d-9d1e-9d14fe8045f0/scratchpad/repro2.py

Config: fit_type=pt, population_size=2, beta=[0.5, 1.0], reps_per_beta=1, exchange_every=10, sample_every=1, burn_in=0, over tests/bngl_files/parabola.bngl. Drive the algorithm to the barrier by hand:

algo = BasicBayesMCMCAlgorithm(cfg); algo.start_run()
algo.current_pset = [a, b]; algo.ln_current_P = [0.0, 1.0]
algo.iteration = [10, 10] # both chains parked at the exchange barrier, E=10
algo.replica_exchange()
algo.try_to_choose_new_pset(1) # the got_result path for chain 1 (the max-beta, sampled replica)

Run with: uv run --extra tests --extra petab python

Observed:
rows before exchange: 1 (header only) iteration [10, 10]
rows after exchange: 2 iteration [10, 10] <- a draw was written at iteration 11, then the counter rewound to 10
rows after got_result for chain1: 3 iteration [10, 11] <- iteration 11 wrote a SECOND draw
samples.txt:
# Name Ln_probability v1__FREE v2__FREE v3__FREE
iter10run1 1.0 4.0 5.0 6.0
iter10run1 1.0 4.0 5.0 6.0

Expected: exactly one samples.txt row per sampled iteration of a chain — iteration 11 of chain 1 should contribute one draw, not two.

Affected configs are those where (k*exchange_every + 1) % sample_every == 0 for some k: any run with sample_every=1 (one extra draw per exchange cycle per max-beta replica), or e.g. exchange_every=19 with sample_every=20, or exchange_every=7 with sample_every=5. The shipped defaults (exchange_every=20, sample_every=100) never hit it.

Verification notes

I read /Users/l119605/Code/PyBNF/pybnf/algorithms/samplers/basic_mcmc.py (committed code — this file is NOT in the dirty working tree) and traced the barrier flow, then ran it.

The counter path is exactly as claimed. try_to_choose_new_pset (lines 194-260) does self.iteration[index] += 1 and then, unconditionally, runs the whole per-iteration block: sample_pset (line 223-225, gated on > burn_in, % sample_every == 0, should_sample), update_histograms, output_results, the progress prints and report_convergence_diagnostics. It is not a pure counter bump.

Sequence for a chain j sitting at the barrier iteration E (E % exchange_every == 0, set at line 256-259 with wait_for_sync=True):

1. got_result sees all wait_for_sync and calls replica_exchange (line ~172).
2. replica_exchange calls try_to_choose_new_pset(j) → counter E→E+1 → full block runs at E+1 with the POST-exchange state.
3. Line 349 self.iteration[j] -= 1 rewinds to E and names the pset iter{E}run{j}.
4. The pset comes back; got_result does accept/reject and calls try_to_choose_new_pset(index) again → counter E→E+1 → the same block runs a SECOND time at E+1, now with the post-accept/reject state.

I looked for the outs and none exist: sample_pset (base.py:382-398) is a bare append with no iteration dedupe; the MCMC-family postprocess (samplers/base.py:82-130) reconciles exchange_every/reps_per_beta/population_size but never relates exchange_every to sample_every; should_sample only restricts WHICH replicas sample (max-beta ones), not how often; the existing tests (tests/test_basic_mcmc_class.py TestReplicaExchange) assert only the swap of current_pset/ln_current_P and deliberately use exchange_every=10000 with sample_every=1000, so nothing pins the double-write as intended. The default pairing (exchange_every=20, sample_every=100) is genuinely safe — 20k+1 is never ≡ 0 mod 100 — which is why it has gone unnoticed.

The claim's one inaccuracy: the duplicated rows are NOT labeled iter{E+1}; sample_pset writes current_pset[index].name, i.e. the last accepted pset's name, so uniq -d on the label is not a reliable detector (and in a real run the two rows can hold different states, since accept/reject happens between them). That is a detail of the repro sketch, not of the defect: one chain iteration still contributes two rows to samples.txt, and therefore double weight in Results/Histograms/*.txt and the credible intervals, plus a doubled row in log_likelihood.txt (which base.py:400-419 documents as row-aligned with samples.txt — so the alignment claim still holds, but both files carry the extra draw).

I ran it and observed 2 rows written for the single iteration 11 of the sampled replica.

Where

pybnf/algorithms/samplers/basic_mcmc.py:349 — severity medium, bug class logic-error.