Restore modern NumPy/SciPy compatibility, fix the long-range LD filter, compile the sampler, modernise packaging/CI

[bvilhjal]

Follow-up to a review of organisation and computational efficiency. Every item here was reproduced before being changed; claims that did not survive measurement were dropped rather than implemented (see §8).

8 commits, 26 files.

1. LDpred did not run at all on a current scientific Python stack

from scipy import isfinite raised ImportError, and 230 further call sites used NumPy aliases in the scipy namespace (sp.dot, sp.array, sp.clip, …) that SciPy has since removed. Verified on scipy 1.17.1 / numpy 2.4.6, where all 32 distinct attributes used are gone.

• Import numpy directly and rename every sp. to np. — all 33 map 1:1, so it is a pure rename. from scipy import stats / linalg are genuine SciPy imports and stay.
• dtype='bool8' (removed in NumPy 2.0) → dtype=bool.

Two latent bugs on the same path that NumPy 2 turns into hard errors:

• calc_risk_scores reshaped pval_derived_effects_prs in place to feed the covariate regressions, mutating the array shared with its caller; the score writer then formatted length-1 arrays with %0.6e. The regressions now use a separate column view and the vector stays 1-d.
• lstsq returns the residual sum of squares as a length-1 array for a column right-hand side, so every derived pred_r2 was an array rather than a float. Squeezed through a small _scalar helper.

2. Long-range LD handling — two bugs, one behaviour change, one new option

The exclusion never fired. load_lrld_dict initialised each chromosome as {'reg_dict': {}} but stored regions as siblings of reg_dict, so is_in_lrld read an always-empty dict and returned False for everything. Sweeping a 1 Mb grid genome-wide flagged 0 positions; chr6:30 Mb (the MHC) came back False. The loop body underneath was dead code that would have raised TypeError on a string key had it ever run.

Fixing that alone would have traded a silent no-op for a KeyError: callers passed a counter over util.chromosomes_list, which runs to 24 entries, while the region table is keyed 1–23. Lookups now key off the group name, and get_lrld_chrom_key resolves 'chrom_X' to 23. get_snp_lrld_status is vectorised over the region list (171× on 500k positions, identical mask).

Skipped variants were reported as zero. Once the exclusion actually ran, a second bug became reachable. A held-out variant stays frozen at whatever it starts with, so the starting value decides both what it contributes to its neighbours' residuals and what weight it is reported with. curr_post_means is only written inside the loop body the skip bypasses, so it kept its initial 0 and the output weight was zero — contradicting the help text and silently dropping ~5% of the genome from the score.

New --long-range-ld {include,inf,exclude}, default inf, replacing the boolean:

value	behaviour
inf (default)	hold out of the sampler, keep in the model at the LDpred-inf effect, report that weight
exclude	hold out with a zero effect, so they leave the model entirely, report a zero weight
include	sample them like any other SNP

--incl-long-range-ld still works as a deprecated spelling of --long-range-ld include and warns; combining it with an explicit --long-range-ld is an error rather than a silent precedence rule. The run summary reports which policy was used.

Behaviour change on real data: with ldpred gibbs at the default, variants in the 34 Price et al. (2008) regions are now held out of the sampler and reported with their LDpred-inf estimates. Scores over those regions will differ from previous releases. Recorded in CHANGES.txt; the policy choice is evidenced in §5.

3. Compiled Gibbs sampler (Numba)

Profiling put ldpred_gibbs at 66.7 s cumulative against 18.0 s for everything in ldpred_inf — the sampler is where a genome-wide run spends its time, and it ran as a per-SNP Python loop.

The sweep is now a standalone kernel compiled by Numba, an optional dependency (pip install ldpred[fast]). Without it the sampler uses an equivalent NumPy implementation. Median of 5 runs, m=20000, ld_radius=100, 60 iterations, with the provenance of each build verified:

	median	vs pre-numba
pre-numba	5.883 s	—
Numba off (fallback)	4.822 s	1.22× — no regression
Numba on	0.500 s	11.8×

Both backends are kept deliberately: running the kernel's scalar loop in the interpreter would be ~12× slower than the code it replaces, so the fallback expresses the inner product as one np.dot instead. LDPRED_DISABLE_NUMBA=1 forces it, which is how the two were compared. The active backend is reported in the run summary.

Random draws are still generated outside the sweep, so the random stream is byte-for-byte what the previous implementation consumed. The backends differ only in summation order: largest relative difference 2.0e-16 across the fixed-N, varying-N and long-range-LD paths.

4. Other measured efficiency fixes

Change	Speed-up	Agreement
calc_auc: O(n²) threshold sweep → rank/Mann-Whitney	504× at n=20000	1e-14 (tie-free)
ldpred_inf: pinv (SVD) → Cholesky solve	3–9.5×	1e-14
get_snp_lrld_status: Python loop → vectorised mask	171×	identical
Gibbs inner loop: dropped two tautological isreal() calls per SNP per sweep	—	bit-identical

calc_auc also fixes a correctness bug: on tied scores the old sweep returned 0.889 where brute force over all pairs gives 0.722.

Scope note: on the integration suite these are worth ~30% end-to-end, but that workload over-weights ldpred_inf. On a genome-wide run with a --f grid, §4 alone is worth ~13%; §3 is what moves the total.

5. Evidence for the long-range LD policy

The exclusion had never actually run before this PR, so there was no evidence for how it should behave — and the goldens cannot supply any, since every simulated test SNP sits below 5000 bp while the earliest real region starts at 25 Mb.

simulations/lrld_policy.py compares five arms; three are reachable from the command line. Mean test R², 10 replicates per cell (simulations/README.md, per-replicate rows in lrld_policy_results.csv):

arm	loading 0.0 (control)	0.5	1.0	2.0
include	0.435	0.014	0.029	0.080
zero (previous behaviour, no longer reachable)	0.324	0.328	0.329	0.330
inf (default)	0.444	0.374	0.300	0.288
exclude	0.326	0.333	0.334	0.333
drop (coordination-step only, not exposed)	0.347	0.332	0.336	0.332

• include does not merely do worse — the sampler diverges, with Σβ² reaching 250–860× the true h², concentrated in the region. That is exactly what LDpred's existing convergence warning reports.
• Excluding is not free. In the control column — a region on the exclusion list that behaves normally — zero costs −0.111 R², exclude −0.109 and drop −0.087.
• inf is the only arm that is never bad, which is why it is the default.
• exclude tracks the full panel rebuild (drop) to within noise at every loading ≥ 0.5 and beats inf there by ~0.03–0.04, which is why it is the one exposed rather than drop. Rebuilding LD on a reduced panel is a coordination-step operation the sampler cannot perform.

Every arm sets its region weights explicitly rather than inheriting the sampler's default, so the arms stay independent of whatever the default happens to be.

Caveat: idealised Gaussian genotypes, one synthetic region, one h², one polygenicity. The relative ordering of inf vs exclude is exactly what could shift on real data, and the default depends on it. The README says plainly this should be repeated on real genotypes with the actual Price regions.

6. Organisation

• setup.cgf was a typo for setup.cfg, so its only content had never been read by anything. Removed.
• setup.py → pyproject.toml. Same distribution: version still single-sourced from ldpred.__version__, same three console scripts, same package data. Classifiers no longer advertise Python 2.7.
• requirements.txt carries lower bounds instead of bare names; new [fast] extra for Numba.
• CircleCI (pinned to a python:3.7.4 image, running only SimpleTests plus one ComplexTests case) → GitHub Actions matrix over 3.9/3.11/3.13 running the unit tests and both golden suites in full, each a second time with the NumPy fallback forced, plus a wheel/sdist build-and-import smoke job.
• Added a .gitignore; the repo had none, so 16 __pycache__/*.pyc files were tracked and every test run dirtied the tree.
• get_LDpred_sample_size existed byte-identically in LDpred_gibbs and LDpred_fast, with the LDpred_fast copy never called. Now lives once in util.

7. Tests

All 16 golden-file tests pass byte-identically on numpy 2.4.6 / scipy 1.17.1, with Numba enabled and with it disabled.

38 unit tests added — the repo had none.

file	count	covers
tests/test_util.py	12	long-range LD lookup, mask, calc_auc incl. ties
tests/test_sampler.py	19	packing helpers, backend equivalence, both sample-size paths, all three long-range LD policies, deprecated-flag resolution
tests/test_lrld_simulation.py	7	simulation stays runnable and reproduces the divergence result

The goldens are unaffected by the §2 behaviour change because none of their 31000 simulated SNPs falls in a real region — which is precisely why the unit tests exist. All three policies were additionally verified end to end through the CLI.

8. Claims that did not survive verification, and were therefore not implemented

• Tiled BLAS-3 LD table. Measured 2.0–3.3×, not the order of magnitude expected. It is sharply tile-sensitive (2.26× at tile 128, 0.49× at tile 2048) because a dense GEMM over a banded problem computes (tile+2r)/(2r+1) times more correlations than needed, and it perturbs D at ~6e-07. Decisively, the LD table is built once and cached: on a 7-value --f grid it is 3.0% of runtime, so even a 3.3× is worth ~2% overall.
• stats.norm.rvs → standard_normal. Only ~1.5× on modern SciPy.

9. Noted, not changed

• The tests package ships ~43 MB of golden fixtures into site-packages because the ldpred-unittest / ldpred-inttest entry points import it. Dropping those two scripts would let the wheel shed all of it.
• coord_genotypes.py:532 reads and parses the entire interpolated genetic-map file inside a loop whose body is entirely commented out, so genetic_map always comes back empty. A missing feature rather than a regression.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Fc3Ljp4v5ANQA7bqVXdGqR