SpectraREML/README.md

# SpectraREML

SpectraREML is a reusable batch AI-REML engine for many continuous traits that share one genomic relationship matrix (GRM).

For each task, it fits

\[
y = X\beta + Zu + e, \qquad
u \sim N(0,\sigma_g^2G), \qquad
e \sim N(0,\sigma_e^2I).
\]

The engine diagonalizes the GRM once, rotates the common design and every unique task-specific covariate once, rotates phenotypes in blocks, and fits independent REML tasks in parallel.

## Features

- C++17 numerical core using oneMKL or OpenBLAS/LAPACKE.
- GRM eigendecomposition once per task set.
- AI-REML with a strong-Wolfe line search; no EM updates.
- Per-task OLS-residual phenotype scaling inside the numerical core, with all
  estimates and the restricted likelihood restored to the input phenotype units.
- Safeguarded quadratic zoom interpolation (central 96% of the bracket),
  bisection fallback, and a last-valid-improving-point fallback when strict Wolfe
  curvature cannot be reached because of numerical roundoff.
- Signed standard-deviation parameterization with a separate one-sided KKT check for the \(\sigma_g^2=0\) boundary.
- No explicit dense \(P\) matrix.
- Variable numbers of extra fixed-effect covariates per task through a CSR index.
- Satterthwaite t/F fixed-effect inference by default, with optional
  Kenward-Roger F tests and a joint test of task-specific covariates.
- Named arbitrary linear contrasts over the common fixed-effect design, using
  the Satterthwaite or Kenward-Roger method selected for the run.
- Atomic block output and safe resume/force semantics.
- Generic Python CLI for manifest creation, validation, execution, provenance signatures, status, and result export.

The numerical core has no knowledge of cohorts, molecular assay types, or domain-specific variable names.

The line search starts at 1, expands by 1.618 up to the configured maximum
step, and uses quadratic interpolation only when its stationary point lies at
least 2% away from both bracket endpoints. An invalid interpolation falls back
to bisection. If the evaluation limit is reached, the search accepts the last
finite covariance-valid point that improved the likelihood; it reports
`line_search_failed` only when no such point exists.

## Repository layout

```text
include/spectra_reml/  public C++ API
src/                   numerical core, batch I/O, and CLI
python/spectra_reml.py generic Python CLI
examples/example.py    fully synthetic input example
scripts/run_server.sh  generic Linux build/run wrapper
tests/                 numerical and file-contract tests
docs/FORMAT.md         binary and tabular file contract
```

## Build with Intel oneMKL

Load the oneAPI environment first:

```bash
source /opt/intel/oneapi/setvars.sh

cmake -S . -B build-mkl \
  -DCMAKE_BUILD_TYPE=Release \
  -DREML_BLAS=MKL \
  -DMKL_INTERFACE=lp64 \
  -DMKL_LINK=dynamic \
  -DMKL_THREADING=sequential
cmake --build build-mkl --parallel
ctest --test-dir build-mkl --output-on-failure
```

For a nonstandard installation, locate `MKLConfig.cmake` and set its directory explicitly:

```bash
cmake -S . -B build-mkl \
  -DREML_BLAS=MKL \
  -DMKL_DIR=/path/to/mkl/latest/lib/cmake/mkl \
  -DMKL_THREADING=sequential
```

Sequential BLAS is recommended because independent tasks are already parallelized by OpenMP.

## Build with OpenBLAS

```bash
cmake -S . -B build-openblas \
  -DCMAKE_BUILD_TYPE=Release \
  -DREML_BLAS=OPENBLAS \
  -DOpenBLAS_ROOT=/path/to/openblas \
  -DLAPACKE_ROOT=/path/to/lapacke
cmake --build build-openblas --parallel
ctest --test-dir build-openblas --output-on-failure
```

The configuration performs a real CBLAS/LAPACKE link check. It supports LAPACKE either inside OpenBLAS or in a separate library.

## Quick start

Create a synthetic bundle:

```bash
python examples/example.py
```

Create a manifest for your own data:

```bash
python python/spectra_reml.py make-manifest \
  --manifest work/manifest.json \
  --grm-bin data/example.grm.bin \
  --grm-id data/example.grm.id \
  --base-x data/base_x.f64.bin \
  --phenotypes data/phenotypes.f64.bin \
  --extra-covariates data/extra_covariates.f32.bin \
  --tasks data/tasks.tsv \
  --extra-offsets data/extra_offsets.i64.bin \
  --extra-indices data/extra_indices.i32.bin \
  --contrast-matrix data/contrasts.f64.bin \
  --contrast-metadata data/contrasts.tsv \
  --output-dir work/blocks \
  --n-samples 3523 \
  --n-base-covariates 8 \
  --n-phenotype-rows 10000 \
  --n-extra-covariate-rows 2000 \
  --n-contrasts 5
```

Validate, run, and export results:

```bash
python python/spectra_reml.py validate --manifest work/manifest.json

python python/spectra_reml.py run \
  --manifest work/manifest.json \
  --engine build-mkl/spectra_reml \
  --threads 28 \
  --blas-threads 1 \
  --block-size 256 \
  --resume

python python/spectra_reml.py finalize \
  --manifest work/manifest.json \
  --output work/results.tsv.gz
```

`results.tsv.gz` retains the complete per-task summary and stores the fixed-effect vector and row-wise packed lower covariance as JSON arrays.

Optional named contrasts are supplied by `--contrast-matrix`,
`--contrast-metadata`, and `--n-contrasts`. Contrast rows act on `base_x`;
task-specific covariate columns receive zero weights automatically. Finalized
results contain contrast estimates, adjusted standard errors, statistics,
numerator and denominator degrees of freedom, and p-values as JSON arrays.
This release uses manifest v2, run-signature v3, block v4, and finalized-output
v3 contracts. It intentionally does not read older block, run-signature, or
finalized-output contracts.

Fixed-effect inference defaults to Satterthwaite. Select Kenward-Roger or turn
inference off with `--fixed-effect-test kenward-roger` or
`--fixed-effect-test none`. Coefficient-wise standard errors, statistics,
denominator degrees of freedom, and p-values are exported as JSON arrays; tasks
with extra covariates also report their joint F test in the summary columns.
An accepted `sigma_g2=0` boundary is refitted as ordinary least squares without
the GRM term and is reported as `boundary_ols`. A retained
`line_search_failed` iterate continues through the selected fixed-effect
inference and is reported as `line_search_conditional` while preserving the
optimizer status.

## Recovery and provenance

Each block is written as four files, with `.complete` renamed last. The Python layer adds `run.signature.json`, which binds the canonical manifest, engine SHA-256, numerical options, thread settings, and block size.

- `--resume` reuses complete blocks only when the signature matches exactly.
- `--force` invalidates the old signature before deleting old blocks and starting a new generation.
- `finalize` refuses blocks that are not bound to the current manifest.
- `--dry-run` does not mutate output state.

## Threading

Use one BLAS thread with multiple outer task threads unless benchmarking shows otherwise:

```bash
export MKL_NUM_THREADS=1
export OPENBLAS_NUM_THREADS=1
export OMP_DYNAMIC=FALSE
```

Then set `--threads` to the physical cores allocated to the process.

See [docs/FORMAT.md](docs/FORMAT.md) for the exact file contract.

## License

BSD 3-Clause. See [LICENSE](LICENSE).