Repository navigation
v8.2.0
Added
-
GEMMA_DIVERGENCES.md§6 cites the jamma-lean proof that the accepted
Newton step contracts toward the root on a well-conditioned smooth peak
(accepts_of_wellCond,newtonLoop_contracts). -
GEMMA_DIVERGENCES.md§6 explains why lambda refinement is bounded by the
coarse grid bracket rather than the golden-section bracket, states the
worst-case bound from jamma-lean, and records measured lambda error by
sample size and peak curvature, with its effect on Wald p-values. -
tests/test_lean_proven_identities.pychecks the production code against
each identity proved in Lean in
jamma-lean, and
GEMMA_EQUIVALENCE.mdlinks the proofs. -
GEMMA parity on BGEN input.
tests/fixtures/bgen_parityholds a BGEN built
from GEMMA's mouse_hs1940 example (160 samples, 1,508 chromosome 1-2
variants, fractional dosages and missing samples on a share of them) and
GEMMA 0.98.5's-gk 1and-lmm 1-4outputs on a BIMBAM copy of the
decoded dosages.scripts/generate_bgen_parity_inputs.pybuilds the inputs,
andgenerate_gemma_fixtures.shregenerates the GEMMA outputs byte for
byte. -
BGEN input from the CLI and
gwas().-bgen FILEreads a BGEN v1.2 file,
with-sample(default: the.bgenpath with suffix.sample) and-bgi
(default:<bgen>.bgi); exactly one of-bfileand-bgenis required.
-p FILEreads a GEMMA phenotype file (whitespace-separated, no header, one
row per sample in genotype order,NAand-9missing) with either input,
and-nselects its column; BGEN input requires it, and a row count that
differs from the genotype sample count is an error.-info FLOATfilters
kinship and association SNPs on INFO over the analysed samples, and enters
the LOCO eigen cache key when it is on.-hweand--backend numpyare
rejected with-bgen, and-infowith-bfile;autoruns BGEN on the
streaming backend.gwas()takesbgen=,sample=,bgi=,
phenotype_file=andinfo=, and its PLINK prefix is now optional.
PipelineConfiggains the matching fields, andPipelineConfig.genotypes()
parses them into aPlinkInputorBgenInput. BGEN LOCO eigen caches key
on the BGEN file identity; PLINK keys are byte-identical. -
GenotypeDataset.open_bgen(bgen, sample, bgi)reads BGEN v1.2 layout-2
files (biallelic, unphased diploid, bit depth 1 to 16; zlib, zstd or
uncompressed) through the dataset API. The counted allele is the first
allele, so the dosage is 2·P(11) + P(12). Variant metadata and block offsets
come from the bgenix.bgiindex and samples from the.samplefile; the
header, index and sample counts, and any embedded sample IDs, are checked
against each other at open. The decoder is C in the_lmm_accelextension
(which now links the system zlib), parallel across variants. zstd files
need the newjamma[zstd]extra below Python 3.14. -
An imputation INFO filter for BGEN input.
SnpStats.infoholds each SNP's
INFO, GCTA's--info(the IMPUTE2 information measure) over the
non-missing analysed samples, computed from the stored integer
probabilities in exact integer sums, so it equals GCTA's value bit for bit.
It is not clamped and can be negative; a monomorphic or all-missing SNP
gets 1. Hard-call genotypes (PLINK) have INFO 1.SnpFilterSpec,
KinshipSnpFilter,LmmRunSpec,compute_kinship_streamingand
compute_loco_kinship_streamingtake aninfo_threshold(default 0.0,
off) that keeps SNPs with INFO >= the threshold, for kinship and
association SNPs alike; a threshold on genotypes without INFO raises
ValueError. The LOCO eigen cache key includes the threshold only when it
is on, so existing PLINK keys are unchanged.
Changed
-
plan_associationandExecutableAssociationPlantake a required
genotype_encoding, so a caller cannot price BGEN as hard calls by
omission. The plan'sstats_block_sizeis the only statistics-pass width;
LmmRunSpec.stats_block_sizeis gone. -
A zstd-compressed BGEN without the zstd module, or any BGEN without the
_lmm_accelC extension, now ends a CLI run with a one-lineError:and
exit code 1 instead of a traceback. Both raise the new
BgenDependencyError, anImportError; the C-extension case was a
RuntimeError. -
The BGEN decoder accumulates the exact INFO sums while it decodes, over the
rowsGenotypeDataset.statsis given or the newinfo_rowsargument of
GenotypeDataset.blocksnames, so INFO no longer costs a second pass over
the quantised probabilities. On 100k samples x 250 variants (8-bit zlib, 18
threads)dataset.statsfell from 75 to 33 ms per block over all rows and
from 155 to 73 ms over 90k rows; the kinship pass with an INFO threshold
fell the same. INFO values are unchanged bit for bit, andblock.infofor
rows other than the block was read with still recomputes from the stored
quantised values.decode_bgen_probabilities_ctakes optional
info_rowsandinfo_sumsarguments (_lmm_accelABI 25). -
PipelineRunner.validate_inputschecks only that the input files exist;
opening the dataset checks the PLINK.bedsize, which was previously
checked and logged twice. The covariate row-count error now says "genotype
data" instead of "PLINK data".EigenCacheComponentsis a plain dict,
since its file-identity keys depend on the genotype format. -
run_lmm_association_numpy_streamingandrun_lmm_locotake a
GenotypeDataset(GenotypeDataset.open_plink(prefix)) instead of a PLINK
prefix. The streaming runner loses itssnp_infoparameter, since the
dataset carries the variant metadata, and itsvalidate_genotypes
parameter, since hard-call validation now follows the dataset's encoding.
One function,prepare_genotypes(dataset, samples, filters, *, stats=None),
prepares genotypes for every run, and batch mode reads the dataset into
memory with the newGenotypeDataset.materialize().run_lmm_association_numpy
keeps its signature; its in-memory matrix is now also checked for values
outside {0, 1, 2, NaN}, and a count of them is logged, as the.bed
paths already did. Results are unchanged. -
compute_kinship_streamingandcompute_loco_kinship_streamingtake a
GenotypeDataset(GenotypeDataset.open_plink(prefix)) instead of a PLINK
prefix, andcompute_loco_kinship_streamingloses itsmetaparameter.
compute_eigen_cache_keytakes the dataset too and reads its file identity
fromdataset.fingerprint(); PLINK cache keys are byte-identical, so
existing LOCO eigen caches still hit.LocoRuncarries the dataset in place
ofbed_pathandmeta. The pipeline and-gkopen the genotypes once,
so-gknow also checks the.bedsize against the.famand.bim.
Results are unchanged. -
SnpMetaandSnpInfoRecordmoved fromjamma.lmm.schemato
jamma.genotype.variants, with no alias at the old path, so genotype
readers can build variant metadata without importing the LMM package.
Results are unchanged. -
jamma.corenow holds only leaf infrastructure and imports no other jamma
package. Domain modules moved to the packages that own them, with no aliases
left at the old paths:jamma.core.snp_statsandjamma.core.snp_filter
are nowjamma.genotype.snp_statsandjamma.genotype.snp_filter;
jamma.core.eigen_planisjamma.lmm.eigen_plan;jamma.core.recompileis
jamma._native;jamma.core.estimatesisjamma.estimates; and
estimate_kinship_memorymoved fromjamma.core.memoryto
jamma.kinship.memory.write_gemma_logmoved fromjamma.utilsto
jamma.gemma_log. Results are unchanged. -
The NumPy reference kernels the C accelerator is held to are now public:
jamma.lmm.compute_numpy.compute_wald_numpy,compute_lrt_numpyand
compute_score_numpy, andjamma.lmm.stats.batch_lrt_pvalues_numpy. They
were underscore-prefixed although tests import them directly. Results do not
change. -
Without the jlinalg C extension,
jamma.jlinalg.get_n_threads()and
set_n_threads()return 1, because the NumPy SNP-statistics fallback is
unthreaded. They used to store a thread count that nothing read.
set_n_threadsstill rejects counts below 1. -
jamma.jlinalgdrops the private_dgemm_numpy,_dgemm_numpy_impl,
_dsyrk_numpy,_dsyrk_numpy_impl, and_eigh_numpyaliases. The NumPy
implementations stay injamma.jlinalg._dgemm,_dsyrk, and_eigh. -
jamma._nativeis the only code that evicts a rebuilt extension from
sys.modules.compile_and_link.compile_extensionno longer does it too. -
read_matrix_parallelraisesValueErrorforn_workers < 1, as
write_matrix_parallelalready did, instead of clamping to one worker. It
scans the input once rather than twice to find chunk boundaries. -
The parallel matrix reader and writer remove their
.jamma_mread_*and
.jamma_mwrite_*temp directories and everything in them on every exit.
A failed early delete of the writer's memmap no longer leaves the directory,
with the memmap inside, beside the output. -
scripts/bench_memory.pyreplacesbench_impute_memory.py,
bench_snp_stats_layout.pyandbench_kinship_memory.pyas the
subcommandsimpute,snp-statsandkinship. The digest scripts share
one hashing and--out/--diffmodule,scripts/_digest_common.py, which
the digest workflows stage beside the head script. Digest keys and values
are unchanged. -
The native build resolves its compile and link flags once per extension,
and every build step reports through oneBuildReportwith a verbose
channel and an always-visible warning channel. Compiler command lines are
byte-identical to before. A runtime recompile after an ABI mismatch now
logs compiler failures and OpenMP retries through the logger instead of
printing them to stderr.python -m jamma.lmm._compile_acceland
python -m jamma.jlinalg._compile_jlinalgprint the same output as
before; acompile_extension()call withoutverbose=Trueno longer
prints the success line. -
Without the C extension, the REML and MLE likelihoods each have one
finisher shared by the grid search and the per-SNP refinement, and the
identity-weighted log-determinant is computed once per chunk. The REML
logl_H1that the Python API returns for-lmm 1can move by one ulp
(at most 2.1e-16 relative, 11 of 500 SNPs ongemma_synthetic); lambdas,
every p-value, and the.assoc.txtbytes onmouse_hs1940are unchanged. -
Without the C extension, the Wald statistics of an intercept-only run
(-lmm 1or-lmm 4, no-c) now come from the same NumPy REML optimiser
as every other mode and covariate count. The separate split-Uab optimiser
and thenumpy_walddispatch path are removed. Results onmouse_hs1940
are byte-identical. On that fixture's one 10,768-SNP chunk the fallback's
association pass is 12 percent slower for-lmm 1and 5 percent slower for
-lmm 4, because the full Uab doubles the per-SNP memory traffic. At
2,000-SNP chunks the general optimiser is 4 to 9 percent faster. For
-lmm 1the per-SNP chunk quote rises by threen_samplesrows plus 18 Iab
cells, so the planner picks smaller chunks. -
The
_lmm_accelextension has one workspace creator and one chunk compute
for every covariate count,create_workspace_candcompute_lmm_chunk_c,
replacing the separate n_cvt=1 and general entry points (ABI 23). The C
side picks the kernel family fromn_cvt, soDispatchPath.FUSED_GENERAL
is gone andDispatchPath.FUSEDcovers every covariate count. The n_cvt=1
workspace now allocates its per-thread scratch once at creation, as the
general one does, and reads the null-modelwcolumn straight fromUtW
instead of a copy.workspace_sizes_creturns persistent and per-thread
bytes only, and an n_cvt=1 quote falls by 8 bytes per sample. Results are
bit-identical. -
The likelihood-ratio test with more than one covariate (
-lmm 2and
-lmm 4) now reads the same per-SNP column layout as the REML search,
instead of copying each SNP's inputs into a second per-thread buffer of
n_index * n_samplesdoubles. That buffer took 22.4 MB per thread at
100,000 samples and 5 covariates, and the memory planner no longer prices
it. The MLE coarse grid also reuses the grid's cached invariant sums.
Results are bit-identical. -
LOCO kinship no longer runs a separate SNP-statistics pass over the BED.
Standard and LOCO kinship now share one chunk generator that filters each
chunk as it is read, so a LOCO run reads the genotypes once fewer. The
batch planner sizes passes against every chromosome the BIM and-ksnps
leave SNPs on, and allocates a per-chromosome accumulator only for a
chromosome that keeps SNPs after filtering.LocoKinshipStream.snp_stats
is readable once the first matrix has been yielded, and a filtering error
surfaces on the stream's first advance. Kinship matrices and statistics
are bit-identical. -
LOCO association reuses the SNP statistics kinship PASS 1 computed over the
analysed samples, so a run with missing phenotypes no longer reads every
chromosome's genotypes a second time. A run on cached eigenpairs computes
the statistics in one streamed pass before the chromosome loop. Results are
unchanged. -
Score-only runs (
-lmm 3) with two or more covariates no longer build the
lambda grid, which only Wald and LRT read: at 100,000 samples and the
default 50 grid points that is 40 MB less per workspace. The memory quote
now comes from the same per-family layout each native workspace allocates
from, so it drops that grid, the REML-only identity sums in-lmm 2and
-lmm 3, and one of four scratch vectors then_cvt=1LRT loop never
allocated. Results are unchanged. -
jamma.lmm.plan_associationtakes anLmmConfigasconfig=in place of
the separatelmm_mode,n_grid,n_refineandmem_budgetarguments,
and onebackend=literal ("auto","numpy","numpy-streaming"or
"loco") in place ofrequested=plusloco=True. The
log_dispatch_choicesargument is gone: the dispatch path is always logged
at debug level, once. -
The pipeline logs the backend request and its source as one line,
Backend request: numpy (config)orBackend request: numpy-streaming (JAMMA_BACKEND), in place of the threeBackend: numpy (...)variants. -
run_lmm_locono longer takesexecutionoreigen_plan. A caller with
its own plan builds aLocoRunand callsrun_loco(run, output_path)
instead;LocoRunrejects a plan wider thanloco.col_chunk_sizeor
without a resolved kinship shape at construction. A pipeline LOCO run now
resolves its samples, plan and kinship shape once, and reads the.bim
once instead of twice. -
PipelineConfigrejects an illegal kinship or eigen source at
construction:-kwith-loco, an unpaired-d/-u,-d/-uwith
-loco, and-widvwith-locoor with-d/-u. These used to wait for
validate_inputs(), after the PLINK and input-file checks, so they now win
over a missing file. On the CLI they exit with the usage-error code 2
instead of 1, with the same message. -
gwas(eigen_dir=...)andPipelineConfig(eigen_dir=...)without
loco=TrueraiseValueError. They used to accept the directory and
ignore it; only the CLI rejected--eigen-diroutside-loco. -
gwas(phenotype_columns=...)defaults to(1,)instead ofNone. -
PipelineResulttakesphenotype_results,n_samples,timing,
n_covariatesandanalyzed_sample_indices.associations,
n_snps_tested,assoc_path,assoc_paths,pve_estimateandpve_se
are now read-only properties derived fromphenotype_results, so reading
them is unchanged but passing them to the constructor is an error. -
AssocComparisonResultholds one comparison per column the mode carries,
plusaf, incolumns; read one withcomparison["l_remle"]. A column the
mode does not carry is absent instead of a vacuous pass, a SNP count
mismatch reports the single columnn_snps, andpassedis derived.
ComparisonResultis frozen. -
compare_assoc_resultsraisesValueErrorfor LRT rows whose beta or se is
not NaN, instead of returning a failed verdict with every column passing. -
load_gemma_assocrejects duplicate SNP IDs, and takesmode=and
require_logl=to check the header, so the math-validation wrapper no
longer parses headers or rows itself. -
The memory plan prices
-lmm 2at three result columns per SNP, the
logl_H1,l_mleandp_lrtarrays the kernels allocate, instead of two.
The chunk sizer and the NumPy and native workspace estimates all read the
column count from the mode's specification. Chunk width can shrink by a few
SNPs in LRT runs; results are unchanged. -
A failing native kernel names its mode as
-lmm Nin the error, for example
Fused general -lmm 4 dispatch failed ..., in place of per-mode labels such
asFused general mode-4 Uab dispatch. -
The C accelerator decodes
lmm_modeonce into the set of tests it runs,
and both covariate families share one null-model argument contract, thread
clamp, and result builder. Two behaviours now match across families: the
general (n_cvt >= 2) path accepts an empty genotype chunk and returns
zero-row columns instead of raisingValueError, and it warns about
betainc non-convergence for standalone Score as the n_cvt = 1 path already
did. The n_cvt = 1 LRT runs through the same loop as Wald and mode 4. No
result bit moves. -
Whole-genome and LOCO eigen generations share one member naming scheme, one
manifest validator and one commit protocol (EigenGenerationin
jamma.lmm.eigen_io). The whole-genome.eigen_manifest.jsonis now
fsynced before it is published, so like the LOCO cache manifest it survives
a power cut. A whole-genome manifest must name its members exactly; a member
name with extra parts is rejected. Manifests written by 8.1.0 still resolve.
LocoConfig.eigen_stemandLocoConfig.eigen_paths, which named a
pre-generation layout nothing wrote, are removed. -
scripts/bench_all_backends.pytimes a second kinship row,kinship_npy,
which runs-gk 1without--legacy-textand so writes the binary.npy
matrix JAMMA ships by default. The existing text row remains the
like-for-like GEMMA comparison, and the new row has no GEMMA counterpart. -
jamma.jlinalg.jlinalg_isareports the SIMD ISA the extension was compiled
for, read from the preprocessor, instead of probing the CPU with CPUID,
XGETBV, andgetauxvalat import. A baseline Linux x86_64 wheel built
without-mavx2said"AVX2"on any AVX2 CPU; it now says"generic".
The AVX2 release wheels still say"AVX2", and AArch64 builds, including
Apple Silicon, still say"NEON". The attribute's name, its possible
values, andABI_VERSIONare unchanged.
Fixed
-
The mutation runner's repository copy no longer fails on a test file
another worker is planting or removing.tests/math_validation/mutations.py
skipstest_planted_*and*.partwhen copyingtests/. The
tests/test_planted_*.pyfiles thattests/test_conftest_c_seam.pywrites
could be listed and then deleted beforecopytreereached them, failing
test_single_mutation_is_isolated_and_detectedwith[Errno 2] No such file or directory(master run 36030355261). -
GEMMA_EQUIVALENCE.md§8 no longer says the chi-squared CDF is linear near
LRT = 0; its slope is unbounded there, which amplifies small LRT
differences.GEMMA_NUMERICAL_EQUIVALENCE_BOUND.mdassumption 7 says why the
Lipschitz range must exclude a zero statistic, and §2 shows that eigenvector
error, which depends on the eigengap, does not reach Pab row 0 orlog|H|. -
GEMMA_EQUIVALENCE.md§4 states the REML formula JAMMA computes, including
thelog|Z'Z|term. §3 says which quantities are invariant to eigenvector
sign flips (Pab andlog|H|, notU'y). §5 and the optimizer docstrings
say up to three Newton steps, not one, and that the accept rule bounds the
result by the coarse grid bracket.GEMMA_NUMERICAL_EQUIVALENCE_BOUND.md§3
gives the lambda tolerance as relative, since the optimizer works in log
lambda. -
Documentation matches the code at
ecbbfd7d.GEMMA_DIVERGENCES.md§14
now covers BGEN dosage precision, missing samples, supported layouts and
INFO edge cases.--backend numpyis documented as rejected with-bgen
only without-loco. The mouse_hs1940 table inGEMMA_EQUIVALENCE.mdis
re-measured:logl_H1matches GEMMA at its printed precision, replacing the
stale ~1.35e-3. Docker-oexamples, compile entry points, the CI matrix,
error messages and API snippets are corrected. -
A stale BGEN index fails before decoding. When the
.bgichromosome,
rsid, position or allele order differs from the variant block, the error
names the field and both values, so the index's labels are never attached
to another variant's dosages. An identifier or allele that is not UTF-8 is
reported as such, not as a header overrun. -
Memory quotes price every genotype read with one model, for PLINK and BGEN
alike: the decoded float64 block plus the larger of the BGEN decode buffers
and the consumer's copy of the block over the analysed rows. LOCO kinship
reserves the kinship chunk's working set of three float64 blocks and two
masks, not one block: tracemalloc measured a PLINK LOCO pass at 34.80 MB
against a 21.76 MB quote. The streaming statistics pass reserves the NumPy
statistics kernel's working set, and both it and the association pass
reserve the read overmin(chunk, n_snps)columns. PLINK quotes rise as
well as BGEN ones;tests/test_memory_ledger_digest.pyrecords every moved
row. -
A genotype block stream frees each block before it reads the next, with
row selection or a progress bar, so a chunk stream holds one chunk at a
time. -
The Docker image loads one Intel OpenMP runtime. It moves to Python 3.12
because numpy-mkl publishes no rebuilt cp311 wheel, and the old cp311 wheel
bundles a secondlibiomp5that can crash at exit.mkl-servicemoves to
2.8.0, matching the fixed wheels, and a build step fails when a process maps
more than one OpenMP runtime. -
The build now uses clang with Intel OpenMP (libiomp5) on hosts that have
clang but no LLVMlibomp-dev, as in images that install clang with
--no-install-recommends. The clang probe passed-fopenmpat link time,
so clang added-lompand the link failed. The build then fell back to GCC,
whoseGOMP_*calls go through libiomp5's compatibility shim and can hit
OMP: Error #13. The probe now compiles with-fopenmpand links libiomp5
by path, as the real build does. It also adds-Ifor theomp.hthat the
intel-openmpwheel installs next to libiomp5, because clang has noomp.h
of its own withoutlibomp-dev. -
Eigen cache members are flushed to disk, and their directory with them,
before the manifest that names them is committed. The manifest was already
fsynced, so after a power cut it could survive while its members came back
empty. A LOCO eigen cache whose.npymember is empty or cut short now
counts as incomplete and is recomputed, instead of failing on first read. -
IncrementalAssocWriterno longer publishes an incomplete.assoc.txt, or
finishes silently without publishing, when a caller catches a write error and
carries on. Once a write, flush,tell()or rollback fails, later
write_arrays_batchcalls raise aRuntimeErrornaming that failure, and a
clean exit raises instead of publishing, so the previous output stays in
place. An interrupt or out-of-memory error during a write cuts the retained
.partialback to its last complete row, so the logged result count matches
the rows on disk. When an earlier failure already discarded the output, the
interrupt warning says no partial output was retained instead of naming a
deleted temp file. Runs that raise no error publish the same bytes as before. -
A kinship or eigen text file read while another run replaced it no longer
leaves a sidecar of the old contents that later reads prefer. The.npy
sidecar now carries the modification time of the text it was built from,
not the time it was written, and is not written at all when the text changed
during the parse. A binary.npywritten after the text still takes
precedence, as before. On filesystems with one-second timestamps, a rewrite
within the same second can still pass as current. -
The Linux build links the C extension against Intel OpenMP when it comes from
theintel-openmppip distribution, as it does with the current numpy-mkl
wheels, instead of falling back to GCC'slibgomp. The finder now also checks
that distribution's installed files and<sys.prefix>/lib, which pip's build
isolation leaves pointing at the installing environment. A new build-time
variable,JAMMA_LIBIOMP5, pins the path explicitly; see
docs/CONFIGURATION.md. -
The runtime rebuild of a missing or stale C extension works again. Since the
loader moved tojamma._native, it passed the directory above the installed
jamma/package, so every rebuild failed with "C source files missing" and
the run fell back to NumPy. After an ABI mismatch, the loader now rebuilds and
asks for a restart instead of retrying: CPython keeps the stale extension
loaded for the life of the process, so only a new process can use the rebuild. -
Reading a kinship or eigen file no longer deletes its
.npyon a read
error. A permission or I/O error was treated as corruption, and a binary-only
write leaves the.npyas the only copy, so one unreadable read destroyed
the matrix. Only a corrupt file (NumPy'sValueError) is removed now. -
LOCO eigendecomposition no longer keeps a chromosome's eigenvectors alive
after association releases them. The worker thread and the submitting loop
each kept a local reference to the last Future, which holds the
eigenvectors, so one extran x nmatrix per idle worker stayed live
outside whatplan_loco_workersprices. With the DSYEVR driver, that
pushed the peak above the priced figure. -
The MLE lambda optimizer (
-lmm 2and-lmm 4, native C and NumPy) now
applies the safeguarded analytic-score Newton refinement REML already had.
On a flat likelihood, golden section alone stopped where rounding decided
its comparisons, so the reportedl_mledepended on the CPU: the weighted
mode-4 contract'sboundary5gave4.729012e-03to4.729206e-03against
a stationary point of4.72909239e-03, and failed on some Linux CI runners.
It now lands on the stationary point regardless ofn_refine. On
mouse_hs1940mode 4,l_mlechanges for 99% of SNPs, by up to1.9e-5
relative: against dense-solve MLE roots for 200 sampled SNPs, the median
error falls from8.7e-6to6.9e-8, which is the output's print
precision.p_lrtchanges for about 8% of SNPs, by up to8.9e-6
relative. REML outputs, Wald and Score are bit-identical. -
jamma -gk 1 -locono longer reserves memory for an eigendecomposition it
never runs. The standalone LOCO kinship path only writes each matrix to
disk, but its pass planner held back the DSYEVR peak (160 GB at 100,000
samples), so a 100,000-sample run on a node with about 480 GB free planned
one chromosome per disk pass (22 passes) where three fit. The LOCO
association path still reserves its eigen workers' peak. The written
matrices are unchanged. -
The association pass no longer collapses to one SNP per chunk on large
sample counts. Since 8.0.0 the chunk sizer subtracted the whole eigenvector
matrix from the per-chunk budget, whose auto ceiling is 40 GB, so above
about 70,000 samples (80 GB of eigenvectors) nothing was left and every SNP
became its own chunk, each re-streaming the matrix through the rotation
GEMM: a 100,000 x 50,000-lmm 1run spent 6 hours in the association
loop. Per-run allocations are priced once by the memory preflight, and the
sizer floors the chunk at 100 SNPs again; a budget that cannot hold the
floored chunk is refused by the preflight rather than planned narrower.
Results are unchanged. -
jlinalgnow wires vendordsyrk,dsyevd, anddsyevrfrom the
scipy-openblas64 library that NumPy's Linux wheels bundle. Its exports carry
ascipy_prefix (scipy_dsyrk_64_) that only thedgemmname table
listed, so on those installsjlinalg.dsyrkandjlinalg.eighfell back to
NumPy whileblas_backendreportedOpenBLAS-ILP64. Linux results move
in the last bits: kinship by at most 7.4e-16 and association outputs by at
most 1.2e-11, relative to each output's largest magnitude. -
A SNP whose projected genotype variance is not positive (a constant
genotype, or one collinear with a covariate) now gets a NaN REML
log-likelihood andl_remleat the lower bound on every path. The
n_cvt=1C path and the NumPy path used to skip the non-positive Pab
diagonal and report a finite likelihood and an optimised lambda, while the
general C path returned NaN. One rule,logdet_diag_termin
_lmm_types.hand_logdet_diaginlikelihood_numpy.py, now holds for
all three, matching GEMMA'sLogRL_f, which takes an unguardedlogof
each diagonal entry. beta, se and the p-values were already NaN; results
for SNPs with positive projected variance do not change. Such SNPs never
pass the polymorphism filter, sojammaCLI output on mouse_hs1940 is
byte-identical. -
The
forced-lower-boundmathematical mutation plants its defect again. Its
patch targeted the grid-evaluation call that the likelihood finisher
refactor rewrote, soscripts/mathematical_mutations.py --allreported it
INCONCLUSIVE and the post-merge slow suite failed. The runner now checks
every patch against the source before running any test and reports a
non-matching one asSTALE PATCH, naming the mutation and file, and a
tier0 test pins each Python patch to exactly one match. -
jlinalg.dsyrkno longer segfaults under MKL ILP64 (libmkl_rt, pipmkl
on Linux x86_64). The dispatcher stored the Fortran symboldsyrk_64_in
the CBLAS slot and called it with CBLAS enums, so MKL dereferenced the
integer101as theuplostring. Fortran and CBLAS names now resolve into
separately typed slots, and a name whose prefix contradicts its slot is left
unwired with a warning instead of being called. The MKL and OpenBLAS ILP64
names take the Fortran path, whichdgemmalready used; Accelerate wiring
is unchanged. -
With
JAMMA_FORCE_NUMPY_FALLBACK, SNP means from float32 genotype chunks
are now accumulated in float64, as the C kernel does, instead of rounded to
float32. Allele frequencies and missing-genotype imputation under the
fallback now match the C path; results with the C kernel are unchanged. -
The build prints each OpenMP retry notice once.
execute_buildreported a
failed OpenMP compile or link through bothon_retryandverbose_print,
so the wheel build andpython -m jamma.lmm._compile_accel, which point
both at stderr, logged the same line twice. The notice now goes to
on_retry, or toerror_printwhen no callback is given, so a caller that
silencesverbose_printstill sees the single-threaded downgrade. -
scripts/demonstrate_equivalence.pyjudges mouse_hs1940 at
tests.fixture_paths.NUMPY_GEMMA_TOLERANCES, the tolerances the parity
suite enforces, and reads every fixture path from the same module. #321
tightened the suite tolambda_rtol=5e-5, pvalue_rtol=2e-4, se_rtol=1e-5, logl_rtol=1e-6, atol=1e-12while the report kept its own copy of the
earlier1e-3, 1e-2, 5e-4, 5e-3, 1e-4, so it certified equivalence at a
p-value tolerance 50x wider than the suite's. The report now has no
tolerance or path of its own to drift. -
The "Eigendecomp memory" log line names the planned driver and the real
reason for it. On a host with no vendor DSYEVD or DSYEVR and
JLINALG_NO_VENDOR_LAPACKunset it read(DSYEVD): ... (JLINALG_NO_VENDOR_LAPACK set, ...), contradicting the "No vendor LAPACK"
line printed just before it; it now reads(numpy): ... (no vendor DSYEVD or DSYEVR, using np.linalg.eigh; estimate is DSYEVD-sized)and the
duplicate line is gone. The forced-numpy line namesnumpytoo. The
planner records the reason where it chooses the driver (EigenDriverPlan
now carriesdriver,required_gbandreason;plan_eigen_driver
takesforced_numpyandinplace_blockerin place ofno_vendorand
inplace_eligible), so the runtime no longer re-derives it. The chosen
driver and every memory estimate are unchanged. -
run_lmm_association_numpy,run_lmm_association_numpy_streamingand
run_lmm_locoreject a phenotype vector containinginfwith
ValueError("prepared phenotypes must contain only finite values"), the
error the CLI already raised. Before,infpassed the missing-sample mask
and every statistic came back NaN with no error. The rule lives once, in
compute_valid_mask, which every entry runs first; the grouped runner's
separate copy is gone. -
LOCO yields a chromosome with no kinship SNPs (one that
-ksnpsor the
MAF filter left empty) in biological order instead of after every other
chromosome, sorun_lmm_locooutput and the reported PVE follow chromosome
order whether the eigenpairs are computed or read from the cache. The
stream writes that chromosome's full kinship into its shared buffer instead
of copyingS_fullonce per empty chromosome, and the "PVE computed from
chromosome N" log states the real reason: earlier chromosomes had no SNPs
to test. -
--mem-budgetnow reaches the-gkkinship accumulation gate and the
-gk -eigeneigendecomposition gate.jamma -gk 1 -eigen --mem-budget 0.001
previously exited 0 and wrote every file while-gk 1 -locoand-lmm
rejected the same budget.compute_kinship_streamingtakesmem_budget
and vetoes the run before the first genotype read, as
compute_loco_kinship_streamingdoes, and the eigendecomposition plans
its driver against the budget as the-lmmpath already did. -
jamma.jlinalg.eigh(K, driver="dsyevr")raisesRuntimeErrorwhen vendor
DSYEVR is not wired, as documented, instead of running DSYEVD. The memory
plan chooses DSYEVR for its O(N) workspace, so a DSYEVD run would touch
pages the plan never reserved.jlinalg_eigh_cnow calls DSYEVR through
one helper on both paths. The pipeline is unchanged:core/eigen_plan.py
only plans DSYEVR whenblas_has_dsyevris set. -
The raw
_jlinalg.dgemmand_jlinalg.dsyrkentry points raise
RuntimeErrorwhen their vendor routine is not wired, asjlinalg.h,
_jlinalg.pyi, andpymodule.cdocument, instead of callingabort().
The publicjamma.jlinalg.dgemmanddsyrkbind NumPy in that state and
were unaffected.JLINALG_NO_VENDOR_DSYRKandJLINALG_NO_VENDOR_DSYEVR
joinJLINALG_NO_VENDOR_DGEMMas test seams that leave one vendor routine
unwired. -
jamma -gkno longer applies-lmm's rules to options it never reads,
and no longer accepts them silently.-gk 1 -lmin -1and
-gk 1 -loco -hwe 0.01exited 2 with-lmm's validation messages, while
-gk 1 -k missing.cXX.txt -snps missing.txtran without a word.-gknow
leaves every-lmm-only option (-k,-d,-u,--eigen-dir,-hwe,
-lmin,-lmax,-snps,-widv,--backend) at its default and logs
one warning naming each one given on the command line,
matching GEMMA's tolerance of flags a mode does not use.-lmmis
unchanged. -
compare_assoc_resultsgates AF with an absolute tolerance of one printing
unit,ToleranceConfig.af_atol = 1e-3, in place ofaf_rtol = 0.05. The
relative gate passed a flipped allele whenever AF was within about 0.013 of
0.5 (0.49against0.51), while the docstring promised that a flip fails.
strict()andrelaxed()no longer widen AF, and the math-validation
wrapper's separateaf_orientationcheck is gone because the oracle now
enforces the same bound; its failure ids read<rs>:af.
Removed
- The PLINK-only genotype API, now that every consumer reads a
GenotypeDataset:PlinkMetadata,get_plink_metadata,PlinkData,
read_genotypes,load_plink_binary,partitions_from_metadataand
stream_genotype_chunksfromjamma.io.plink(and theirjamma.io
re-exports);jamma.genotype.snp_stats.collect_streamed_snp_stats;
SnpMeta.from_plink_meta; and theGenotypeSourceprotocol with its
MatrixSource,BedSourceandbed_chunk_sourceimplementations. Use
GenotypeDataset.open_plink(prefix)for metadata (variants,samples,
partitions),blocks()andstats()for streaming, andmaterialize()
for an in-memory copy. jamma.kinship.SnpStatsCacheand thesample_scopefield onSnpStats.
LocoKinshipStream.snp_statsis now always aSnpStatsover the rows the
kinship pass filtered on, neverNone.jamma.core.memory.MemoryLedger,estimate_lmm_memory,
estimate_streaming_memory,kinship_cost,eigen_costandlmm_cost,
and theirjamma.corere-exports.ExecutableAssociationPlan.price()
(plan_association(...).price(eigen=...)) is the one memory quote, the
figure every preflight gates on.jamma.core.eigen_plan.square_matrix_gb
is gone too;array_gbnow lives injamma.core.memory, and
array_gb(n, n)replacessquare_matrix_gb(n). The batch association
quote drops4 * nfloat64 it carried and streaming never did (3.2 MB at
100,000 samples). Streaming and LOCO quotes are unchanged.jamma.lmm.ioandjamma.lmm.resultsare merged into
jamma.lmm.assoc_output, which owns the.assoc.txtrow.
IncrementalAssocWriter(path, mode)takes aModeSpecin place of a
test_typestring, andwrite_arrays_batchno longer takeslmm_mode.
AssocResultmoves there fromjamma.lmm.stats, andModeSpec.header
replacesHEADERS.format_assoc_line,TEST_TYPE_MAP,ACCUM_KEYS,
RESULT_FIELDS,FORMAT_COLUMNSandLocoResultare gone;
run_lmm_locoreturnsLmmRunResult.PipelineTimingmoves to
jamma.pipeline_config. Output bytes are unchanged.jamma.core.cleanup_memory. It rangc.collect()twice and logged RSS,
and nothing in JAMMA or its scripts called it. Callgc.collect()and
jamma.core.log_memory_snapshot()directly instead.jamma.core.hardware.get_hardware_context()andHardwareContextnow
live inscripts/_hardware_context.py, beside the benchmark scripts that
are their only users.- Three one-shot developer scripts whose job is done, none of which a
workflow, hook, test or doc invokes.scripts/codemod_test_builders.py
migrated tests ontorotated_lmm_inputsin #266 and now finds nothing to
rewrite.scripts/verify_loco_thread_limits.pychecked LOCO BLAS scope
ownership for #361 and #362 through a private function.
tests/test_loco_worker_resources.pycovers the same completion, early
close and input-failure exits, plus a solve failure, with a fake BLAS
controller.
scripts/lmm_accel_sections.pyand its test were the section census for
splitting_lmm_accel.c, and that split is done. Git history keeps all
three. jamma.lmm._compile_accel.compile_extensionno longer takesdiagnose,
which added clang-Rpassor gcc-fopt-info-vec-allvectorization
reports to the dev rebuild. Nothing in the repository passed it. To see the
reports, run withverbose=True, copy a printed compile command, and add
the flag by hand.compile_extensioninjamma.lmm._compile_acceland
jamma.jlinalg._compile_jlinalgno longer takeson_retry. Nothing in the
repository passed it. The OpenMP retry notice still prints, on the same
stream as the rest of the build output.jamma.lmm.likelihood.reml_log_likelihoodandmle_log_likelihooddrop
thenc_totalkeyword and evaluate the null model only, the one form
production calls. The alternative-model REML moves to
tests/reference/likelihood.pyasreml_log_likelihood_alt, and the
alternative-model MLE fast path_mle_p_yy_scalar_ncvt1, which only tests
reached, is deleted.jamma.lmm.specialno longer exports the scalarbetaincandchi2_sf.
Production calls onlybetainc_batchandchi2_sf_batch; the scalar forms
move totests/reference/special.pyas their test oracles. No result
changes.