Fitting Functions

`fit()`

The primary estimation entry point. Runs FOCE, FOCEI, or SAEM depending on options.method.

pub fn fit(
    model: &CompiledModel,
    population: &Population,
    init_params: &ModelParameters,
    options: &FitOptions,
) -> Result<FitResult, String>

Parameters: - model: Compiled model from parse_model_file() or parse_full_model_file() - population: Population data from read_population_for() (or read_nonmem_csv() for a Gaussian-only model — see Parsing) - init_params: Initial parameter values - options: Estimation configuration

Returns: FitResult with parameter estimates, standard errors, and per-subject diagnostics.

Example:

use ferx_core::api::read_population_for;

let model = parse_model_file(Path::new("model.ferx"))?;
let (population, _) = read_population_for(&model, &None, "data.csv", None, None, None, &[])?;
let options = FitOptions::default();

let result = fit(&model, &population, &model.default_params, &options)?;
println!("OFV: {:.4}", result.ofv);

`fit_from_files()`

Convenience wrapper that handles parsing and data reading. data_path is None to rely solely on the model’s own [data] block (see data.qmd); when both are given, data_path overrides the model’s, with a warning recorded on the result.

pub fn fit_from_files(
    model_path: &str,
    data_path: Option<&str>,
    covariate_columns: Option<&[&str]>,
    options: Option<FitOptions>,
) -> Result<FitResult, String>

Example:

let result = fit_from_files(
    "model.ferx",
    Some("data.csv"),
    None,          // Auto-detect covariates
    None,          // Default options
)?;

`run_model_with_data()`

Full pipeline: parse model file, read data, fit. Returns both the fit result and the population. data_path is None to rely solely on the model’s own [data] block.

pub fn run_model_with_data(
    model_path: &str,
    data_path: Option<&str>,
) -> Result<(FitResult, Population), String>

Uses the [fit_options] from the model file.

`run_model_simulate()`

Simulation-estimation: parse model, generate data from [simulation] block, fit.

pub fn run_model_simulate(
    model_path: &str,
) -> Result<(FitResult, Population), String>

Requires a [simulation] block in the model file.

Running many fits

A tool that runs hundreds of fits (bootstrap, SCM, model search) has two things to get right that a single interactive fit does not: how the thread budget is split, and how much each fit prints.

Splitting the thread budget (`PoolPlan`)

fit() already parallelises over subjects, so a naive outer par_iter over replicates nests one Rayon pool on top of another and the two levels compete for the same workers. PoolPlan makes the split explicit and builds the outer pool correctly:

use ferx_core::{fit, FitOptions, PoolPlan};

// 200 replicates on an 8-thread budget: 8 fits at a time, each single-threaded.
let plan = PoolPlan::from_budget(8, 200);
assert_eq!((plan.replicates(), plan.threads_per_fit()), (8, 1));

let mut options = FitOptions::default().quiet();
plan.apply_to(&mut options);          // sets FitOptions::threads

let results: Vec<_> = plan.install(|| {
    use rayon::prelude::*;
    (0..200).into_par_iter().map(|_i| fit(&model, &population, &init, &options)).collect()
})?;

The budget is spent on the outer level first, because replicate-level parallelism is embarrassingly parallel while the per-subject level synchronises every iteration:

replicates      = min(total_threads, n_units)     (at least 1)
threads_per_fit = total_threads / replicates      (at least 1, floored)

So from_budget(8, 200) is 8 x 1, from_budget(8, 2) is 2 x 4, and from_budget(8, 3) is 3 x 2 — the division floors rather than oversubscribing. total_threads = 0 means the engine default (available cores minus one, capped at 8); both levels are always at least 1, so a plan never silently means “let Rayon decide”.

PoolPlan::install() is the part you should not hand-roll. ferx builds its own pools with a 32 MiB worker stack (FIT_RAYON_STACK_SIZE) because wide ODE+IOV analytic-gradient models overflow the platform-default 2 MiB Rayon stack. An outer pool built from a bare rayon::ThreadPoolBuilder::new() inherits that default and faults on exactly the models that need the bigger stack — and only under a tool, where nothing else looks.

With threads_per_fit == 1, apply_to() gives the documented single-threaded inner mode: the fit runs its per-subject loops on a one-worker pool and reports n_threads_used == 1. Estimates do not depend on the thread count — the per-subject likelihood is reduced in subject order — so a replicate fitted on one thread is bit-identical to the same replicate fitted on four.

Quiet fits

FitOptions::verbose defaults to true, which suits one interactive fit and is unusable for 200 of them. FitOptions::quiet() turns it off:

let options = FitOptions { outer_maxiter: 200, ..Default::default() }.quiet();

A quiet fit writes nothing to stdout or stderr. Nothing is lost: every non-fatal message the engine produces goes into FitResult::warnings (and warnings_structured) regardless of verbose, so the tool reads them off each result and prints its own one-line-per-replicate progress instead.

Ranking candidates: BIC variants and strictness

A search (covariate, structural, variability — the #1175 tools) is generate a candidate → fit it → rank it. Two things the ranking step needs live in ferx-core (#1177): the BIC convention to rank on, and a gate that decides whether a candidate’s fit is trustworthy enough to be ranked at all.

BIC variants (`bic`)

FitResult::bic is the classical OFV + p · ln(n_obs). Pharmpy’s iivsearch and modelsearch rank on the Delattre et al. (2014) mixed BIC instead, which penalises the random-effects class of parameters on ln(n_subjects) and the rest on ln(n_obs) — ranking IIV structures on the observation-count BIC systematically favours the wrong model. bic() gives all four of pharmpy.modeling.calculate_bic’s variants from a finished FitResult:

use ferx_core::{bic, BicType};

let mixed  = bic(&result, BicType::Mixed);   // Pharmpy's default
let iiv    = bic(&result, BicType::Iiv);     // free Ω elements only, ln(n_subjects)
let random = bic(&result, BicType::Random);  // every parameter on ln(n_subjects)
let fixed  = bic(&result, BicType::Fixed);   // == result.bic

`BicType`	penalty
`Mixed`	`n_random · ln(n_subjects) + n_fixed · ln(n_obs)`
`Iiv`	`n_omega · ln(n_subjects)`
`Random`	`n_parameters · ln(n_subjects)`
`Fixed`	`n_parameters · ln(n_obs)`

Class membership follows Pharmpy: every free Ω / Ω_IOV element is random; a θ is random when it reaches — directly or through intermediate [individual_parameters] assignments — a parameter expression that carries an ETA or KAPPA, and fixed otherwise; σ (and a block_sigma correlation) is fixed unless the residual error itself carries an ETA (iiv_on_ruv). The [event_model], [binary_model] and [markov_model] expressions count too, since they accept ETA directly: scale = TVSCALE * exp(ETA_SCALE) in a TTE-only model with no [individual_parameters] makes TVSCALE random-class, the same as writing that line as an individual parameter. (A θ that meets its η only in a Form-C y = … readout or an [odes] hazard line is not seen and stays fixed-class.) A [covariate_nn]’s weights are θ too: an output head read through an η-bearing parameter (CL = TYPICAL_PK.CL * exp(ETA_CL)) makes the shared hidden-layer weights and that head’s row and bias random-class, while another head’s own weights stay fixed-class. The tally is recorded on FitResult::bic_inputs at fit time, so a .fitrx bundle can be ranked without re-parsing the model; a bundle saved before the tally existed gets NaN rather than a wrong penalty.

The four variants are pinned in unit tests to calculate_bic()’s own outputs on Pharmpy’s pheno example (59 subjects, 155 observations, three θ that all reach an η-bearing parameter, two Ω, one σ): mixed 611.707, fixed 616.537, random 610.741, iiv 594.431.

Strictness (`check_strictness`)

FitResult::converged is a bool. It is never true at an objective that is not a usable number — NaN, infinite, or the clamped divergence sentinel all report converged: false with a W_NONFINITE_OBJECTIVE warning (#1303) — but it does not distinguish a genuine optimum from an init stall (#751), a boundary estimate, an ill-conditioned covariance step or a near-singular correlation matrix. Under automation every one of those becomes a model-selection error: a candidate that never left its initial estimates is ranked on an OFV that says nothing about the model. Strictness is the pyDarwin-style gate, and the verdict carries the reasons, so a search report can say why a candidate was excluded instead of silently dropping it:

use ferx_core::{check_strictness, Strictness};

let gate = Strictness {
    require_covariance: true,
    ..Strictness::default()
};
let verdict = check_strictness(&result, &gate);
if !verdict.passed {
    for reason in &verdict.failures {
        eprintln!("excluded: {reason}");
    }
}

field	default	fails when
`require_converged`	`true`	`converged == false` — which also covers an internal runaway-guard hit, and a non-finite objective (`NaN`, infinite, or the clamped divergence sentinel), both of which demote `converged`
`require_covariance`	`false`	the covariance step did not run or failed. A SIR fallback that delivered credible intervals is accepted — unless its importance weights collapsed onto a single draw (Kish ESS below 2, every interval of zero width), which fails like a failed step; an accepted fallback stores no matrix, so the threshold gates below then report `skipped`
`max_condition_number`	`Some(1000)`	`cov_condition_number` exceeds it (∞ for a singular correlation matrix)
`max_correlation`	`Some(0.95)`	any off-diagonal of the covariance’s correlation form exceeds it in absolute value — the whole matrix, not just θ, read on the natural scale (a `block_omega` / block-κ segment is mapped off its Cholesky coordinates by the delta method first, so the number is the ω correlation NONMEM’s `.cor` would show — the layout comes from `FitResult::omega_is_diagonal` / `kappa_is_diagonal`, and a bundle saved before those were recorded is read as stored); the offending pair is named by its packed coordinate (`CL ~ log_chol_ETA_CL`)
`reject_on_boundary`	`true`	a θ is pinned to a declared bound — the same predicate `bootstrap`’s `skip_estimate_near_boundary` applies
`reject_init_stall`	`true`	the fit never left its initial estimates (the #751 signature): the outer optimizer’s own escape verdict when the result carries it (`FitResult::left_init`, recorded by the NLopt outer loop), otherwise no free θ/Ω/σ moved more than 1 % of its initial value; a fit with `n_parameters == 0` cannot stall, and one whose only free coordinates are κ / mixture / `block_sigma` entries with no optimizer verdict is `skipped`. This assumes a cold start — a candidate warm-started from its parent’s estimates legitimately converges where it began, so turn the gate off for such a search

While either threshold is enabled, a fit whose covariance step had to floor a Hessian eigenvalue — the Covariance step regularized: eigenvalue floor applied … covariance_regularized warning — fails as well (#1512); the non-finite cross-partial note that shares that code does not. The floor replaces a direction of negative or near-zero curvature with a finite one, so the matrix both thresholds read no longer shows the problem: a collapsed peripheral compartment (V2 → 0, Q free) on warfarin read a condition number of 2.98 next to a TVQ RSE of 293519 %, and passed. The floor fires only when the Hessian is indefinite or its smallest eigenvalue is below 1e-10 of its largest (or below 1e-12 outright, when the whole spectrum is that small), the analogue of NONMEM’s “R matrix algorithmically singular”, so the check applies at every severity the warning reports: that grade measures how much the floor moved the reported standard errors, not whether the model is identified.

The two thresholds are evaluated only when their input exists — a fit without a covariance matrix has no condition number to test. Such a case lands in verdict.skipped rather than failing or silently passing; pair the thresholds with require_covariance = true to make them mandatory. Strictness::none() turns every gate off for a “rank everything, report everything” run.

The predicates are also available on their own — estimate_near_boundary, stalled_at_init, max_abs_correlation — for a report that prints the numbers rather than the verdict.

`build_fit_inputs()`

Extract initial parameters and fit options from a parsed model, separating parsing from estimation for timing purposes.

pub fn build_fit_inputs(
    parsed: &ParsedModel,
) -> Result<(ModelParameters, FitOptions), String>

Example:

let parsed = parse_full_model_file(Path::new("model.ferx"))?;
let (init_params, options) = build_fit_inputs(&parsed)?;

let (population, _) = read_population_for(
    &parsed.model, &parsed.covariate_decls, "data.csv", None, None, None, &[],
)?;
let result = fit(&parsed.model, &population, &init_params, &options)?;

Fitting Functions

fit()

fit_from_files()

run_model_with_data()

run_model_simulate()

Running many fits

Splitting the thread budget (PoolPlan)

Quiet fits

Ranking candidates: BIC variants and strictness

BIC variants (bic)

Strictness (check_strictness)

build_fit_inputs()

`fit()`

`fit_from_files()`

`run_model_with_data()`

`run_model_simulate()`

Splitting the thread budget (`PoolPlan`)

BIC variants (`bic`)

Strictness (`check_strictness`)

`build_fit_inputs()`