Skip to content

Releases: bowo1698/masreml

MasRemL v0.5.0

Choose a tag to compare

@bowo1698 bowo1698 released this 18 May 13:47

Add a new dataset for forward prediction scenario

MasRemL v0.4.0

Choose a tag to compare

@bowo1698 bowo1698 released this 18 May 06:10

masreml 0.4.0

Repository: https://github.com/bowo1698/masreml
Documentation: https://bowo1698.github.io/masgenomics-docs/


RELEASE STATUS

Pre-1.0 maintenance release on top of 0.3.1. The public R-level
API is unchanged. The Rust kernel ships with extensive new
module-level and function-level documentation, inline math
comments at every REML / EMMAX / matrix algorithm, and a new
Rust unit-test layer (22 tests, all passing). One numerical bug
in the pedigree A-matrix builder is fixed; see CHANGES below.


CHANGES

Bug fix

  • matrix/pedigree.rs: build_a_ped_internal now uses the correct
    direction of the Henderson (1976) recursion for off-diagonal
    kinship entries. Previously the inner loop used the parents of
    j (the older individual); this collapsed parent-offspring
    kinship to zero whenever the older individual was a founder
    with unknown parents. The recursion now uses the parents of i
    (the younger individual), matching Mrode (2014) eq. 2.1:

    a_ii = 1 + F_i,   F_i = 0.5 * a[sire_i, dam_i]
    a_ij = 0.5 * (a[j, sire_i] + a[j, dam_i])   for j < i
    

    Diagonal entries (a_ii, inbreeding F_i) were correct in 0.3.1
    and are unchanged.

    Downstream impact: any analysis that fitted a mixed model with
    the pedigree A matrix as a relationship component on data
    containing founders should be re-run. SNP-only GBLUP, MH-only
    GBLUP, GWAS, and GWABLUP results from 0.3.1 are unaffected
    (those paths do not call build_a_ped).

Documentation

  • Every Rust source file under src/rust/src/ now carries a
    module-level docstring (//!) covering purpose, algorithms,
    references, and design invariants. The previous TODO
    placeholders in masgenomics-docs/internals/ have been replaced
    with auto-extracted content; the docs site picks up changes to
    the Rust source on each render.
  • Inline math comments document every step of the BayesA-style
    conditional updates in reml/he_regression.rs, reml/ai_reml.rs,
    reml/em_reml.rs, and gwas/emmax.rs. Each sampling step or
    iteration includes the closed-form formula, the role of every
    helper quantity, and a one-line justification for the step.
  • matrix/snp_additive.rs and matrix/pedigree.rs gained function-
    level rustdoc covering the mathematical definition, GWABLUP
    weight option, and error semantics.

Testing

  • A Rust unit-test layer is now shipped under src/rust/src/
    (#[cfg(test)] blocks). 22 tests cover:
    • matrix::mh_additive — Da (2015) three-genotype encoding,
      frequency-weighted row shrinkage (partial damping factor,
      dimension preservation, zero-frequency skip).
    • matrix::snp_additive — VanRaden centering, allele-frequency
      computation, G symmetry, mean(diag(G)) under HWE,
      monomorphic-marker error path.
    • matrix::pedigree — diagonal correctness on a full-sib
      family, parent-offspring kinship (regression test for the
      bug above), inbreeding from a selfing event, dimension
      validation.
    • Existing gwas::emmax and gwas::smoother tests are retained
      and continue to pass.
      All tests run with cargo test inside src/rust/. They are
      excluded from the release binary via #[cfg(test)] and add no
      runtime cost.

Internal renaming (no math change)

  • The post-encoding step in matrix/mh_additive.rs was renamed
    from project_sum_to_zero to frequency_weighted_row_shrinkage
    to reflect what the operation actually does. The previous name
    was misleading: the function partially damps the
    frequency-weighted row sum by a factor (1 - Q/F), it does not
    drive it to zero. The math is unchanged and is identical to the
    0.3.1 implementation; outputs are bit-for-bit identical for any
    data that was processed by 0.3.1.

COMPATIBILITY

  • R API: unchanged. All public functions in 0.3.1 keep their
    signatures, argument names, and return types in 0.4.0.
  • FFI / extendr exports: unchanged.
  • Pedigree A-matrix consumers: numerical results will differ from
    0.3.1 because of the Henderson-recursion fix. The 0.4.0 values
    are the intended ones; pBLUP / A-based variance-component
    estimates should be re-run on data of interest.

DEPENDENCIES

No changes since 0.3.1. The Rust crate continues to use
extendr-api 0.7, ndarray 0.15 + ndarray-linalg 0.16 (OpenBLAS),
rayon 1.8, statrs 0.16, libm 0.2, thiserror 1.0.


SYSTEM REQUIREMENTS

  • R >= 4.2
  • Rust toolchain: rustc >= 1.70, cargo (installed via rustup)
  • OpenBLAS development headers on Linux/macOS (system); bundled
    static OpenBLAS on Windows.
  • Builds and tested on macOS (arm64) and Linux (x86_64) for this
    release. Windows build paths exist but have not been re-
    validated against the documentation changes.

MasRemL v0.3.1

Choose a tag to compare

@bowo1698 bowo1698 released this 16 May 05:44

masreml 0.3.1

Repository: https://github.com/bowo1698/masreml
Documentation: https://bowo1698.github.io/masgenomics-docs/


RELEASE STATUS

Pre-1.0 release. The public R-level API (masreml, cv_masreml,
run_gwas, gwablup, build_G_snp, build_D_snp, build_G_mh,
build_A_ped, varcomp, predict, evaluate_prediction, and the S3
summary/print methods) is considered usable for SNP and multi-
allelic GBLUP, pedigree BLUP, EMMAX-style GWAS, and GWABLUP.
Default algorithm choices, the PCG threshold, and the binary
Laplace step may still be tuned before 1.0; any change that
affects fitted variance components or BLUPs at working precision
will be called out in the release notes.

Validation in this release:

  • No testthat suite is shipped yet. End-to-end behaviour is
    exercised through examples/ and the bundled load_data() demo
    (n = 200, 500 SNPs, 250 MH blocks, n_qtl = 20, h^2 ~ 0.5) for all four
    relationship-matrix types, all REML algorithms, both
    continuous and binary traits, EMMAX GWAS, and GWABLUP.
  • A unit-test layer covering variance-component recovery on
    simulated data and parity checks against an established REML
    implementation is planned for 0.4.x.
    Users running production analyses should sanity-check the
    fitted variance components and prediction accuracy against an
    independent pipeline on a subset of their data before scaling
    up.

IMPLEMENTATION HIGHLIGHTS

  • Rust backend via extendr-api 0.7 (with the ndarray feature),
    organised into module directories: reml/, gwas/, solver/,
    matrix/, utils/. A build.rs handles platform-specific linkage
    to OpenBLAS (system OpenBLAS on Linux/macOS, static OpenBLAS
    on Windows).
  • The single Cholesky factorisation of V = sigma_g^2 * G +
    sigma_e^2 * I (held in the FactorizedV struct) is reused across
    the REML gradient computation, the BLUP solve, and the EMMAX
    test statistic. This is the main reason an AI iteration and a
    full genome-wide scan can share most of their O(n^3) cost.
  • Solver dispatch by sample size: dense Cholesky (via
    ndarray-linalg 0.16 / OpenBLAS) below the n ~ 10000 threshold,
    preconditioned conjugate gradient at or above it. The choice is
    exposed at R level via solver = "auto" / "cholesky" / "pcg".
  • Algorithm dispatch in method = "auto": continuous traits resolve
    to AI; binary traits resolve to HE on the Laplace-augmented
    working response. AI, EM, HE, and the HE+AI ("HI") combined
    routine are also selectable explicitly from R.
  • ID alignment is driven by R-side named vectors: y must be named,
    marker rownames must match y names, and the C-level routines
    assume this contract. Two marker input modes converge on the
    same solver path: raw markers + G built internally by Rust, or
    pre-built G matrices passed directly from R. For train/test
    splits the ref_W / ref_mh arguments compute allele frequencies
    from training only to prevent leakage.
  • rayon 1.8 is used inside Rust for embarrassingly parallel work
    (per-marker GWAS test statistics, per-fold cross-validation);
    BLAS-level parallelism is delegated to OpenBLAS. OMP_NUM_THREADS
    and RAYON_NUM_THREADS control the two layers independently.
  • Errors are returned through a typed thiserror chain so that R
    receives specific messages (non-PD relationship matrix, PCG
    non-convergence, etc.) instead of generic Rust panics. log 0.4
    is used internally for diagnostics.

NUMERICS AND REPRODUCIBILITY

  • Linear algebra: ndarray 0.15 + ndarray-linalg 0.16 backed by
    OpenBLAS. statrs 0.16 for distributions; libm 0.2 for the
    special functions used in the Laplace approximation.
  • Deterministic given identical inputs, identical BLAS, and an
    identical thread count. PCG iteration counts and therefore the
    working-precision residual may differ across thread counts;
    switch to solver = "cholesky" if exact reproducibility across
    machines with different core counts is required and n permits.
  • Binary-trait fits use a single-step Laplace approximation
    (algorithm string "Laplace-1step (HE)"); the Laplace step is
    not iterated to full convergence by design. The fit returns
    fitted probabilities and metrics on the observed scale; total
    GEBVs remain on the liability scale for compatibility with
    downstream genetic-evaluation tooling.
  • Release build uses opt-level = 3, LTO = true, codegen-units = 1.

SYSTEM REQUIREMENTS

  • R >= 4.2
  • Rust toolchain: rustc >= 1.70, cargo (installed via rustup)
  • OpenBLAS development headers on Linux/macOS (system); bundled
    static OpenBLAS on Windows. The configure script probes for
    these at install time.
  • Suggested R packages: MASS (used by some helpers), devtools and
    rextendr for the development workflow.
  • Builds and tested on Linux (x86_64) and macOS (arm64).
    Windows build paths exist (Makevars.win.in, masreml-win.def)
    but the Windows toolchain has not been validated for this
    release.