Metadata-Version: 2.4
Name: blupf90-wrapper
Version: 0.3.0
Summary: Unofficial Python wrapper for the BLUPF90 family of programs with R-style formula notation.
Author-email: Rajesh Neupane <rajeshneupane7@gmail.com>, "Juan P. Nani" <juanonani@gmail.com>, Umit Bilginer <umitbilginer44@gmail.com>
Maintainer-email: Rajesh Neupane <rajeshneupane7@gmail.com>
License: MIT
Keywords: animal-breeding,quantitative-genetics,reml,blupf90,livestock,variance-components,random-regression,heritability
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: POSIX :: Linux
Classifier: Operating System :: MacOS
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3 :: Only
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Typing :: Typed
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.20
Provides-Extra: dev
Requires-Dist: pytest>=8.0; extra == "dev"
Requires-Dist: pytest-cov>=5.0; extra == "dev"
Requires-Dist: ruff>=0.5; extra == "dev"
Requires-Dist: mypy>=1.10; extra == "dev"
Requires-Dist: build>=1.0; extra == "dev"
Provides-Extra: pandas
Requires-Dist: pandas>=2.0; extra == "pandas"
Provides-Extra: polars
Requires-Dist: polars>=1.0; extra == "polars"
Dynamic: license-file

# `blupf90` — a Python interface for the BLUPF90 suite

> **Disclaimer.** This package is an **independent, unofficial** wrapper.
> It is **not** affiliated with, endorsed by, or maintained by the
> University of Georgia Animal Breeding and Genetics group (the authors
> of the BLUPF90 family of programs). "BLUPF90" refers to the underlying
> Fortran programs, which are distributed separately by UGA under their
> own terms — see <https://nce.ads.uga.edu/wiki/doku.php>. This wrapper
> is distributed on PyPI as **`blupf90-wrapper`** to make the distinction
> explicit; the Python import name remains `blupf90` for brevity, in the
> same style as PyYAML (`import yaml`) or beautifulsoup4 (`import bs4`).

`blupf90` is a Python package that lets quantitative geneticists fit
animal models with the [BLUPF90 family of
programs](https://nce.ads.uga.edu/wiki/doku.php) (`renumf90`,
`blupf90+`) using a familiar **R / `statsmodels`-style formula
syntax**.  You write

```python
import blupf90 as bf

model = bf.univariate(
    "body_weight ~ contemporary_group + sex + cov(age_days)"
    " + ped(animal) + pe(animal)",
    data="phenotypes.txt",
    pedigree="pedigree.txt",
    columns="phenotypes.cols.txt",
)
result = model.fit()
print(result.h2)              # +0.40 +/- 0.02
```

…and the package handles `renum_trait.par` rendering, folder layout,
local execution or SLURM submission, log parsing, and a tidy result
object — so you spend your time on the science, not on shell
plumbing.

It is meant for everyday lab work *and* for reproducible publications:
every run leaves a self-contained folder on disk with the par file,
data symlinks, BLUPF90 logs, and a parsed summary.

---

## Table of contents

1. [Features](#features)
2. [Why this package](#why-this-package)
3. [Installation](#installation)
4. [Package layout](#package-layout)
5. [Formula grammar](#formula-grammar)
6. [Quick start — univariate](#quick-start--univariate)
7. [Quick start — bivariate](#quick-start--bivariate)
8. [Batches — many traits, one SLURM array](#batches--many-traits-one-slurm-array)
9. [The `columns` mapping](#the-columns-mapping)
10. [Configuration knobs](#configuration-knobs)
11. [What gets written on disk](#what-gets-written-on-disk)
12. [Command-line interface](#command-line-interface)
13. [Random regression models](#random-regression-models)
14. [Testing](#testing)
15. [Compatibility](#compatibility)
16. [Contributing](#contributing)
17. [Citation](#citation)
18. [Licence](#licence)

---

## Features

- **R / `statsmodels`-style formula syntax** for univariate,
  bivariate, and general multi-trait animal models.
- **Automatic `renum_trait.par` rendering** with correct column
  indexing, effect ordering, and multi-trait `(CO)VARIANCES` blocks.
- **`OPTION se_covar_function`** blocks for heritability,
  repeatability, and genetic / permanent-environment / residual /
  phenotypic correlations emitted automatically.
- **Batch API** (`univariate_set`, `bivariate_set`, `model_set`) that
  submits many models as a single SLURM array job.
- **Log parser** that turns BLUPF90+'s `blup_vc.log` into a typed
  `Result` object with `.h2`, `.rg`, `.Va`, `.Vpe`, `.Ve`, ….
- **Random regression** (Legendre polynomial) support for
  longitudinal traits — growth curves, lactation curves, and other
  time-varying phenotypes.  See [`README_RRM.md`](README_RRM.md).
- **Zero runtime dependencies** for the parsing / rendering core
  (NumPy required only for random regression).

---

## Why this package

BLUPF90 (Misztal et al., 2002–present) is the de facto standard for
variance-component estimation and genetic evaluation in dairy and beef
cattle, pigs, poultry, aquaculture, and many other livestock species.
Its file-based interface — the `renum_*.par` parameter file with its
many `OPTION` keywords — is **flexible but unforgiving**:

- column numbers must match the data file exactly,
- random-effect blocks must be ordered correctly,
- `OPTION se_covar_function` formulas are written by hand,
- batch submission across many traits or many trait pairs typically
  ends up as ad-hoc bash + sed.

`blupf90` removes that friction by lifting the same workflow into a
small, well-typed Python API:

| Manual BLUPF90                   | `blupf90` Python                    |
| -------------------------------- | ----------------------------------- |
| Write `renum_trait.par`          | `bf.univariate(formula, ...)`       |
| `sbatch run.slurm` per model     | `models.submit()` (one array job)   |
| `grep`/`awk` the log             | `result.h2`, `result.rg`, `.Va`, …  |
| Per-trait folder + symlinks      | `model.setup()`                     |

The API is **data-agnostic** — there is no built-in column map for any
particular schema.  Bring your own data file, pedigree, and a tiny
column-position sidecar, and every formula name is resolved at build
time.

---

## Installation

The package is pure Python (no compiled extensions) but requires the
BLUPF90 executables (`renumf90`, `blupf90+`) on `PATH` for the
`run` / `submit` / `fit` actions to do real work.  Building the model,
inspecting it, and parsing existing logs needs only Python.

```bash
# From PyPI (once released):
pip install blupf90-wrapper

# From a checkout of this repository:
pip install .

# Or, for a development install:
pip install -e .[dev]
```

Note the **distribution name is `blupf90-wrapper`** but the **import
name is `blupf90`** — you `pip install blupf90-wrapper` once, then
write `import blupf90 as bf` in every script.

BLUPF90 itself is distributed by the [UGA Animal Breeding and Genetics
group](https://nce.ads.uga.edu) under its own (free for non-commercial
research) licence.  Install or `module load` it separately.

### Requirements

- Python ≥ 3.9
- `numpy` ≥ 1.20 (for the random-regression module only)
- BLUPF90 executables on `PATH` (for `run`, `submit`, `fit`)

The parsing / rendering core has no third-party runtime dependencies;
`polars` or `pandas` are only needed if you want to post-process the
TSV that `ResultSet.write_tsv()` produces.

---

## Package layout

```
blupf90/
├── __init__.py    public re-exports + version
├── config.py      Config (paths, columns, SLURM defaults), read_columns()
├── terms.py       Term ABC + FixedClass, Covariate, PedigreeRandom,
│                  DiagonalRandom, LegendreFixed, RRMPedigree,
│                  RRMPermanent, parse_formula()
├── legendre.py    orthonormal Legendre polynomial helpers
├── model.py       Model (formula -> par-file rendering -> run / collect)
├── runner.py      Runner (local), SlurmJob, SlurmArrayJob
├── result.py      VarComp, Result, ResultSet, ModelSet
├── api.py         statsmodels-style functions:
│                    univariate(), bivariate(), from_formula(), fit(),
│                    univariate_set(), bivariate_set(), model_set()
└── __main__.py    `python -m blupf90 <subcommand>` CLI
```

The split mirrors the conceptual layers of an analysis: **configure →
specify → render → run → parse**.  Every layer is independently usable
(you can `Result.from_log(model, path)` on its own to re-parse a
finished run, for example).

---

## Formula grammar

```
<response>[ + <response>]  ~  <term>[ + <term>]*
```

where each `<term>` is one of

| token                            | semantics                                              | par-file effect          |
| -------------------------------- | ------------------------------------------------------ | ------------------------ |
| `name`                           | cross-classified fixed effect                          | `cross` class            |
| `cov(name)`                      | linear covariate                                       | `cov` (regression)       |
| `ped(name)`                      | additive genetic via pedigree, *V<sub>a</sub>*         | random, NUMTYPE 2        |
| `pe(name)`                       | diagonal random, permanent environment, *V<sub>pe</sub>* | random, NUMTYPE 1      |
| `leg(cov, order)`                | fixed Legendre-polynomial regression                   | `order + 1` covariates   |
| `rrm_ped(cov, order, animal)`    | additive-genetic random regression                     | coupled random group     |
| `rrm_pe(cov, order, animal)`     | permanent-environment random regression                | coupled random group     |

- **Bivariate** models just put two responses on the LHS:
  `"weaning_weight + yearling_weight ~ ..."`.
- A `cov(X)` term whose name is also a response is silently dropped
  (perfect collinearity), and the drop is noted in `model.summary()`.
- At most one `ped(...)` and one `pe(...)` per model.  The package
  emits `OPTION se_covar_function` blocks for heritability, repeatability,
  and (for bivariate models) the genetic / permanent-environment /
  residual correlations automatically.
- Random regression (`leg`, `rrm_ped`, `rrm_pe`) is documented
  separately in [`README_RRM.md`](README_RRM.md).

---

## Quick start — univariate

### Example 1: Beef body weight (repeatability model)

A single-trait animal model for **body weight** measured at multiple
ages on the same animal:

```python
import blupf90 as bf

model = bf.univariate(
    "body_weight ~ contemporary_group + sex + age_group"
    " + cov(age_days) + ped(animal) + pe(animal)",
    data="phenotypes.txt",
    pedigree="pedigree.txt",
    columns="phenotypes.cols.txt",     # <-- column-position sidecar
)

print(model.summary())   # effects, function definitions, what will be reported
print(model.par())       # rendered renum_trait.par (no files written)

result = model.fit()     # setup folder + renumf90 + blupf90+ + parse log
print(result.summary())
print("h2:", result.h2)           # VarComp(h2 = +0.40 +/- 0.022)
print("Va:", result.Va)           # VarComp(Va = +95.21 +/- 4.10)
print("Vpe:", result.Vpe)
print("Ve:", result.Ve)
```

### Example 2: Dairy milk yield

Same syntax, different domain — a lactation-mean milk yield model:

```python
model = bf.univariate(
    "milk_yield ~ parity + herd_year_season + dim_class"
    " + cov(dim) + ped(animal) + pe(animal)",
    data="phenotypes.txt",
    pedigree="pedigree.txt",
    columns="phenotypes.cols.txt",
)
result = model.fit()
print(result.h2)                  # VarComp(h2 = +0.30 +/- 0.02)
```

### What's in a `Result`

| Attribute        | Type             | Meaning                                          |
| ---------------- | ---------------- | ------------------------------------------------ |
| `.converged`     | `bool`           | parsed from "Final Estimates" block              |
| `.rounds`        | `int`            | AI-REML iterations                               |
| `.n_records`     | `int`            | observations used                                |
| `.Va`            | `VarComp`        | additive genetic variance + SE                   |
| `.Vpe`           | `VarComp`        | permanent-environment variance + SE              |
| `.Ve`            | `VarComp`        | residual variance + SE                           |
| `.h2`            | `VarComp`        | heritability (univariate) from `se_covar_function` |
| `.repeatability` | `VarComp`        | repeatability (univariate)                       |
| `.rg`            | `VarComp`        | genetic correlation (bivariate)                  |
| `.rpe`           | `VarComp`        | permanent-environment correlation (bivariate)    |
| `.re`            | `VarComp`        | residual correlation (bivariate)                 |
| `.h2_t1`, `.h2_t2` | `VarComp`      | heritability of each trait (bivariate)           |

`VarComp` is a frozen dataclass with three fields — `mean`, `smean`
(sample mean), `se` — and a readable `repr` like `+0.40 +/- 0.022`.
Convenience attributes such as `result.Va_SE` are also provided.

---

## Quick start — bivariate

A two-trait genetic-correlation model between **weaning weight** and
**yearling weight** (beef selection index inputs):

```python
model = bf.bivariate(
    "weaning_weight + yearling_weight ~ contemporary_group + sex"
    " + cov(age_days) + ped(animal) + pe(animal)",
    data="phenotypes.txt", pedigree="pedigree.txt", columns="phenotypes.cols.txt",
)

result = model.fit()
print("rg :", result.rg)            # genetic correlation +/- SE
print("rpe:", result.rpe)
print("re :", result.re)
print("h2 :", result.h2_t1, result.h2_t2)
```

The same syntax applies to any two traits — e.g. milk yield /
somatic-cell score in dairy, or backfat / ribeye area in beef
carcass evaluation.

---

## Batches — many traits, one SLURM array

For a publication-grade analysis you typically fit many traits and many
trait pairs in parallel.  The batch helpers below build a `ModelSet`
that shares a single `Config` (one data file, one pedigree, one column
map) so the whole batch becomes **one `sbatch --array` submission**.

### `univariate_set` — same RHS, one model per trait

Heritability for a panel of production traits:

```python
TRAITS = [
    "birth_weight",     # early growth
    "weaning_weight",   # pre-weaning growth
    "yearling_weight",  # post-weaning growth
    "backfat",          # carcass quality
    "ribeye_area",      # carcass quality
]
RHS = ("contemporary_group + sex"
       " + cov(age_days) + ped(animal) + pe(animal)")

ms = bf.univariate_set(
    traits=TRAITS,
    rhs=RHS,
    data="phenotypes.txt", pedigree="pedigree.txt",
    columns="phenotypes.cols.txt",
    modules=["intel-oneapi", "mkl", "blupf90"],   # for `module load` on SLURM
    root="runs/",
)

job_id = ms.submit()         # one sbatch --array=1-N, one node per trait
# ...wait until `squeue -j <job_id>` is empty...
results = ms.collect()        # parses every blup_vc.log
print(results.summary())
results.write_tsv("heritabilities.tsv")
```

`heritabilities.tsv` is a long-format file with one row per
*(model, derived function)* pair plus `Va`, `Va_SE`, `Vpe`, `Vpe_SE`,
`Ve`, `Ve_SE`, `converged`, `rounds`, `n_records` repeated on every row
for easy pivoting.

### `bivariate_set` — every unordered trait pair

Genetic correlations across the whole trait panel:

```python
ms = bf.bivariate_set(
    traits=TRAITS,                # 5 traits -> C(5, 2) = 10 bivariate models
    rhs=RHS,
    data="phenotypes.txt", pedigree="pedigree.txt",
    columns="phenotypes.cols.txt",
    modules=["intel-oneapi", "mkl", "blupf90"],
    root="runs/",
    # pairs=[("weaning_weight", "yearling_weight"), ...]   # optional subset
)

ms.submit()
results = ms.collect()
results.write_tsv("genetic_correlations.tsv")
```

Each bivariate model lands in `<root>/pair_<a>__<b>/`, so univariate
and bivariate batches over the same traits never collide.

### `model_set` — arbitrary list of formulas

When models in the batch don't share an RHS, or when mixing univariate
and bivariate runs:

```python
ms = bf.model_set([
    "weaning_weight ~ contemporary_group + ped(animal) + pe(animal)",
    "weaning_weight + yearling_weight ~ contemporary_group + cov(age_days)"
    " + ped(animal) + pe(animal)",
    "calving_interval ~ parity + herd_year_season + ped(animal)",
], data="phenotypes.txt", pedigree="pedigree.txt", columns="phenotypes.cols.txt")
```

---

## The `columns` mapping

Every name used in a formula — responses, fixed effects, covariates,
the argument of `ped(...)` and `pe(...)` — must appear in the
`columns` mapping with its 1-based position in the data file.  Two
forms are accepted everywhere:

**1. A Python dict** (explicit):

```python
columns = {
    "animal":              1,
    "contemporary_group":  2,
    "sex":                 3,
    "age_days":            4,
    "age_group":           5,
    "birth_weight":        6,
    "weaning_weight":      7,
    "yearling_weight":     8,
    "backfat":             9,
    "ribeye_area":        10,
}
```

**2. A path to a sidecar file** — parsed by `bf.read_columns(...)`:

- `*.cols.txt` — two-column `name<TAB>index` file (recommended; mirrors
  `OPTION fields_passed_to_output` output)
- `*.csv` / `*.tsv` — header row of the data file itself (column index
  inferred from header position)
- `*.json` — `{name: int}` object

If a formula references a name not in `columns`, `Model.__post_init__`
raises a clear error at construction time — before any par file or
folder is created.

---

## Configuration knobs

Every API entry point accepts the same set of optional keyword
arguments, all forwarded to `Config.from_paths(...)`:

| keyword      | default                          | meaning                                     |
| ------------ | -------------------------------- | ------------------------------------------- |
| `modules`    | `[]`                             | `module load` line in the SLURM script      |
| `slurm`      | `{time=24:00:00, mem=64G, ...}`  | per-job SLURM resource overrides            |
| `missing`    | `-9999`                          | `OPTION missing N` in the par file          |
| `root`       | current working directory        | parent of every per-model folder            |
| `em_rounds`  | `50`                             | EM-REML warm-up rounds before AI-REML       |
| `maxrounds`  | `500`                            | maximum AI-REML iterations                  |
| `conv_crit`  | `"1d-10"`                        | convergence criterion (Fortran notation)    |
| `tol`        | `"1d-14"`                        | positive-definite tolerance                 |

---

## What gets written on disk

For each model in `<root>/<model_name>/`:

```
renum_trait.par     # generated by Python
data_clean.txt      # symlink to the data file
pedigree.txt        # symlink to the pedigree file
renum.log           # renumf90 stdout
blup_vc.log         # blupf90+ stdout (parsed by Result.from_log)
renf90.par          # produced by renumf90
renf90.dat          # renumbered data
renaddXX.ped        # renumbered pedigree
solutions           # EBVs / fixed-effect solutions
```

In `<root>/` itself:

```
logs/               # SLURM .out / .err
_job_<name>.slurm   # single-model SLURM script
_job_models.slurm   # batch / array SLURM script
```

For univariate batches the model folder is named after the trait
(`<root>/body_weight/`); for bivariate batches it is named
`<root>/pair_<a>__<b>/`.

---

## Command-line interface

```
python -m blupf90 show     <formula> --data ... --pedigree ... --columns ...
python -m blupf90 par      <formula> --data ... --pedigree ... --columns ...
python -m blupf90 setup    <formula> --data ... --pedigree ... --columns ...
python -m blupf90 run      <formula> --data ... --pedigree ... --columns ...
python -m blupf90 submit   <formula> --data ... --pedigree ... --columns ...
python -m blupf90 collect  <formula> --data ... --pedigree ... --columns ...
```

Every subcommand takes the formula as its single positional argument.
`--columns` accepts any of the sidecar formats listed above.  Add
`--kind univariate` (or `bivariate`) to enforce a runtime check.

For batches use the Python `ModelSet` API directly — the CLI handles
one model at a time on purpose.

---

## Random regression models

For longitudinal phenotypes — growth curves in beef, pigs, poultry,
and aquaculture; lactation curves in dairy; length-at-age in fish —
the classical repeatability model with `pe(animal)` can inflate
heritability estimates by absorbing within-animal autocorrelation
into the additive genetic variance.  The package provides three
random-regression tokens (`leg`, `rrm_ped`, `rrm_pe`) that fit smooth
per-animal curves using orthonormal Legendre polynomials:

```python
model = bf.univariate(
    "body_weight ~ contemporary_group + leg(age_days, 2)"
    " + rrm_ped(age_days, 1, animal) + rrm_pe(age_days, 1, animal)",
    data="phenotypes.txt", pedigree="pedigree.txt", columns="phenotypes.cols.txt",
)
result = model.fit()
# h2 is now a function of age_days; reconstruct h2(t) from the raw
# variance components in blup_vc.log (see README_RRM.md section 6).
```

Full documentation — mathematical background, model choice guidance,
worked example, and troubleshooting — is in
[`README_RRM.md`](README_RRM.md).

---

## How `submit()` works

- `Model.submit()` writes a one-shot SLURM script that `cd`s into
  `<root>/<model>/` and runs `renumf90 < renum_trait.par` followed by
  `blupf90+ < renf90.par`.  No Python re-entry is required because the
  par file is already on disk.
- `ModelSet.submit()` writes **one array script** for every model in
  the set; `$SLURM_ARRAY_TASK_ID` indexes into a bash array of model
  names baked into the script.  One node per task, all tasks run in
  parallel, results are gathered with `ModelSet.collect()` once
  `squeue` clears.

Both scripts honour the `modules` and `slurm` arguments passed in at
construction time.

---

## Hot-reload in Jupyter

If you edit `blupf90/*.py` from a long-lived kernel, hot-reload the
package without restarting:

```python
import importlib, sys
for name in list(sys.modules):
    if name == "blupf90" or name.startswith("blupf90."):
        importlib.reload(sys.modules[name])
```

The result parser is robust to this pattern — it never relies on
class-identity comparisons that would silently break after a reload.

---

## Testing

Smoke tests live under `tests/` and exercise the formula parser,
Legendre polynomial engine, par-file renderer, and data-file
augmentation logic — no BLUPF90 executables required.  Run with:

```bash
pip install -e .[dev]
pytest tests/ -v
```

---

## Compatibility

- Tested against RENUMF90 / BLUPF90+ version 1.169 (2024–).
- The generated par files use only the widely-supported `OPTION` set:
  `sol se`, `EM-REML`, `maxrounds`, `conv_crit`, `tol`, `method VCE`,
  `se_covar_function`, `missing`, `remove_all_missing`, `alpha_size`,
  and (for random regression) `RANDOM_REGRESSION data` / `RR_POSITION`
  / `(CO)VARIANCES_PE`.
- The wrapper does not currently emit ssGBLUP (`OPTION SNP_file`),
  maternal-genetic (`RANDOM animal maternal`), or `OPTION heterogeneous`
  variance blocks.  These can be added by subclassing `Term` and
  registering the subclass in `terms.TERM_FACTORIES`.

---

## Contributing

Bug reports, feature requests, and pull requests are welcome.  When
opening a PR please:

1. add a smoke test under `tests/` covering the new behaviour;
2. keep the parsing / rendering core free of runtime dependencies
   (NumPy is only imported by `legendre.py` and the RRM setup path);
3. update `README.md` / `README_RRM.md` if the public API changes.

The codebase follows PEP 8 with a 100-column line limit and uses
`ruff` for linting.

---

## Citation

If you use this package in a publication, please cite the underlying
BLUPF90 family of programs.  The standard reference is:

> Misztal, I., S. Tsuruta, D. A. L. Lourenco, Y. Masuda, I. Aguilar,
> A. Legarra, and Z. Vitezica (2002–present). *Manual for BLUPF90 family
> of programs.* University of Georgia.
> <https://nce.ads.uga.edu/wiki/doku.php>

For specific algorithms (e.g. AI-REML, ssGBLUP, single-step) cite the
appropriate primary literature as listed in the BLUPF90 manual.

---

## Licence

The Python wrapper (this package) is released under the MIT licence
(see [`LICENSE`](LICENSE)), Copyright (c) 2024-2026 Rajesh Neupane,
Juan P. Nani, and Umit Bilginer.

**BLUPF90 itself is not included** in this distribution and is
maintained separately by the University of Georgia Animal Breeding and
Genetics group under its own terms. This wrapper is not affiliated
with, endorsed by, or maintained by that group. All trademarks and
program names belong to their respective owners.
