Optimise QIS sampling; compute degeneracy in log space

[virgesmith]

Four independent changes to the sampling hot path and the summary statistics. Every one was benchmarked and checked for output equivalence against main before being kept.

QIS::sample — drop the per-call std::map

sample is called once per person per marginal, and each call constructed a std::map<int64_t, int64_t> to track how dimensions are remapped as the array is sliced, then passed it by value into the recursion. The map is a dense lookup over dims.size() entries, so a vector indexed by the original dimension does the same job with no allocation. Also reserves the two dims_to_* vectors and skips the array copy when there are no fixed dimensions to slice.

pick — don't re-sum the distribution

pick began with an accumulate over the distribution. But the sum of a slice is exactly the reduced value the parent just picked from, so it can be threaded down the recursion instead.

This is done for QIS only. The same change in QISI is not made: there the marginals are double, the reused sum differs from a re-accumulation in the last bits, and that occasionally flips a pick at a boundary — at pop=51200 over 8^4 it moved chiSq from 400.11 to 394.57 and produced a different population for the same seed. It also measured no faster there (938ms vs 942ms), so there is nothing to trade for it. QIS's marginals are int64_t and the sums are exact in double, so no such drift is possible.

QISI does keep the signature change — pick takes a pointer and length rather than a const std::vector&, which removes a copy of the sliced array in the 1-D base case. Behaviour there is unchanged.

Results

solve_m, min of 5 runs, 4-D problem with three 2-D marginals:

	main	this branch
pop=100k, 8^4	138.8 ms	116.6 ms	-16%
pop=200k, 16^4	571.6 ms	529.6 ms	-7%
pop=50k, 4^4	48.5 ms	38.4 ms	-21%

Output — population array, chiSq, pValue, degeneracy — is identical to main in every case. QISI output is identical to main for every skips value tested (0, 7, 100, 1000).

degeneracy — fix a silent zero

degeneracy computed n!/prod((a_k+1)!) as a running product of tgamma calls. tgamma overflows to inf for any argument above ~170, i.e. any cell count of 170 or more, and one inf in the denominator collapses the entire product to zero. A 1000-cell array holding a single count of 500 returned 0 for a value that is in fact astronomically large.

Computing the same expression via lgamma and exponentiating once removes the intermediate overflow. Results that were representable before are bit-identical (20 cells of 5 each: 1.7357e-39 both ways).

The stale // not convinced that this is correct comment is replaced with a statement of what the function actually computes — the formula itself is unchanged, only how it is evaluated.

chiSq

Walks the underlying storage rather than incrementing an Index. Equivalent for the same-shaped arrays it is called with, and it is only called once per solve, so this is tidying rather than a speedup.

Testing

• pytest: 21 passed
• C++ _unittest: 1094 tests, 0 failures
• Equivalence and timings above measured by building main and this branch side by side as standalone harnesses driving QIS::solve and QISI::solve directly.

🤖 Generated with Claude Code

[codecov]

Codecov Report

❌ Patch coverage is 93.75000% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 94.68%. Comparing base (48186d9) to head (3aa7458).

Files with missing lines	Patch %	Lines
src/QIS.cpp	94.44%	1 Missing ⚠️
src/QISI.cpp	85.71%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #66      +/-   ##
==========================================
- Coverage   95.03%   94.68%   -0.35%     
==========================================
  Files          26       26              
  Lines        1027     1035       +8     
==========================================
+ Hits          976      980       +4     
- Misses         51       55       +4     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.