Structural model search

Maturity: alpha — see Feature Maturity for what this means.

ferx modelsearch is Pharmpy’s modelsearch: a search over the structural PK model — absorption route, number of peripheral compartments, transit compartments and lag time — where every candidate is one pk template swap away from its parent. Where covsearch inserts and deletes [covariate_model] lines, this tool rewrites the pk line and, with it, the [parameters] and [individual_parameters] declarations the new template needs, pruning the ones it no longer reads.

ferx modelsearch warfarin.ferxsearch --directory warfarin-modelsearch --threads 8

What to search over comes from a .ferxsearch file. Its ABSORPTION, PERIPHERALS, TRANSITS and LAGTIME statements are the space, [rank] the criterion, and a [modelsearch] section picks the algorithm and how η is given to the parameters a candidate introduces:

base = "warfarin.ferx"
data = "warfarin.csv"

[space]
mfl = """
ABSORPTION([INST,FO]); PERIPHERALS(0..2)
TRANSITS([0,1,3], NODEPOT); TRANSITS(N)
LAGTIME([OFF,ON])
"""

[modelsearch]
algorithm    = "reduced_stepwise"     # or "exhaustive_stepwise", "exhaustive"
iiv_strategy = "absorption_delay"     # or "add_diagonal", "no_add"

[rank]
type   = "bic"                        # the mixed BIC, Pharmpy's default
cutoff = 3.84                         # optional: the improvement a candidate must show over the base

[strictness]
require_converged = true

[run]
retries   = 2
threads   = 8
cache_dir = "warfarin-modelsearch"

key	default	meaning
algorithm	reduced_stepwise	one of Pharmpy’s three, below
iiv_strategy	absorption_delay	η on a new lag time or mean transit time only; add_diagonal gives every new parameter an η; no_add none. Pharmpy’s fullblock is not offered — a block over the new and existing η is iivsearch’s move
[rank] type	bic	any of the rank types; every candidate is scored on it
[rank] cutoff	none	a candidate is selected only when it beats the base by at least this much; the base wins otherwise

The space

The four structural categories, and what each value becomes:

MFL feature	ferx	how
ABSORPTION(INST)	✓	the *_iv templates
ABSORPTION(FO)	✓	the *_oral templates; *_transit when a chain is present
ABSORPTION(ZO)	✓	the *_iv template with zero_order(dur = …) into central — an [odes] candidate
ABSORPTION(WEIBULL)	✓	the *_iv template with weibull(td = …, beta = …) into central — an [odes] candidate
ABSORPTION(SEQ-ZO-FO)	✗	not one input term but a depot of its own, filled at a constant rate and emptied by ka; no template, and the search does not generate one — refused by name when the file is loaded
ELIMINATION(FO)	✓	every template’s own CL/V
ELIMINATION(ZO / MM / MIX-FO-MM)	✓	an [odes] override of the flux out of central — see the [odes] candidates
PERIPHERALS(0 / 1 / 2)	✓	one_cpt_* / two_cpt_* / three_cpt_*
TRANSITS(n, NODEPOT)	✓	*_transit with a FIXed count, theta TVNTR(n, …) FIX; TRANSITS(0) is a move off a chain
TRANSITS(N)	✓	*_transit with the count estimated
TRANSITS(n, DEPOT)	✗	needs a separate depot with its own ka — an [odes] model
LAGTIME(ON / OFF)	✓	the lagtime= binding

Anything marked ✗ is a hard error naming the feature when the file is loaded — the coverage check of every .ferxsearch file — so a search never narrows itself silently. Two combinations the coverage check cannot see are refused per candidate, with a note in the report rather than a candidate: transit absorption with two peripheral compartments (there is no three_cpt_transit template), and the pairs Pharmpy itself declares meaningless (below).

The base model must be a pk model — one_cpt_oral, two_cpt_iv, and so on. An ode(...) model has no template to swap and is refused; so is the *_ig family, which has no MFL coordinate. A *_transit base is read at the count its n binding declares: a literal (n=3) or a parameter behind a FIXed θ is TRANSITS(3), a parameter behind a free θ is TRANSITS(N), and a non-integral fixed count is refused by name. A base’s f= (bioavailability) binding is not a coordinate and is carried through every swap — every template reads it — so F and its θ and η survive a change of compartment count.

The algorithms

All three are Pharmpy’s (tools/modelsearch/algorithms.py), and the candidate enumeration is checked against Pharmpy’s own — see Compared with Pharmpy.

exhaustive. Every combination of at most one feature per category, every one derived from the base in a single layer and fitted in parallel. With k candidate moves across c categories that is at most ∏(n_i + 1) − 1 models.

exhaustive_stepwise. One feature at a time, in every order. Layer 1 applies each allowed feature to the base; layer k applies each feature still allowed to every model of layer k − 1, each child seeded from its parent’s final estimates (Pharmpy’s update_initial_estimates). Two orders of the same features are two candidates, because they start from different parents — the runner’s canonical-hash dedup catches the case where two orders happen to render the same text. A variance the parent collapsed (an ω at the optimizer’s rail, 6 × 10⁻⁶) is seeded at the smallest value the engine will start from, 10⁻⁵, so the child can still be fitted; on the warfarin example a lag whose η collapsed would otherwise have left every two-compartment lag candidate seeded from it unfittable. A block_omega can collapse the same way through its correlations — every declared variance ordinary, one Cholesky diagonal at the rail (the vancomycin base’s ETA_V2) — and that block alone is nudged by 10⁻⁵ on its diagonal until its factor clears the rail; a standalone ω beside it, and a FIXed block, are left as they are, and a block that does not clear the rail after a bounded number of nudges goes through verbatim so the engine’s own message names it. SeedInits on its own stays faithful to the fit; the floor and the nudge belong to a search deriving a child.

reduced_stepwise — the default, as in Pharmpy. As above, but before a layer is extended the models that carry the same feature set are collapsed to the one with the lowest OFV, and only it continues. Every model is still ranked at the end; the table’s continued column says which one went on.

Which move may follow which

Pharmpy’s _is_allowed, kept verbatim so the enumeration anchors:

• a category is moved along once per path — LAGTIME(ON) after LAGTIME(ON), or ABSORPTION(INST) after ABSORPTION(FO), is never a step;
• the peripheral count starts at the smallest count in the space; a later peripheral step may take any other count not yet tried;
• TRANSITS(0, NODEPOT) is never a move (it is first-order absorption);
• the pairs Pharmpy declares meaningless are never combined: a lag time or a transit chain on a bolus, a lag time with a transit chain, and first-order absorption with a one-transit chain.

ferx applies the pair table in two places Pharmpy does not: to the exhaustive combinations (Pharmpy builds and fits INST + LAGTIME(ON) there), and to the candidate’s whole structure rather than only to the features a path applied — so a base that already carries a lag time is never given a transit chain, and a bolus base is never given a lag. Both are stated below.

The base model

The base is the input model when its structure lies in the space. When it does not — a one-compartment input in a PERIPHERALS(1..2) space — the input is fitted first and a base is derived from it with Pharmpy’s least number of transformations: for every category the space names whose values do not include the input’s, the space’s first value (its smallest count for PERIPHERALS). The base is fitted and used as the root; both appear in the table, and the base’s own features are removed from the candidate moves, so a feature the base already carries is never a step.

The practical consequence: a category the space names is a set of allowed values, and the base is moved onto it. LAGTIME(ON) alone means “every model has a lag time”; write LAGTIME([OFF,ON]) to explore it. A category the space does not name keeps the base’s value. (Pharmpy differs here — it fills an unnamed category with its default, so an oral base in a space without an ABSORPTION statement becomes a bolus.)

What a candidate becomes

A candidate is its parent’s text with the pk line replaced and the parameters reconciled. Roles the parent already binds keep their variable — by slot, so v1=V on a two-compartment line satisfies v=V on a one-compartment one, and q=Q becomes q2=Q under three compartments. Roles it does not bind get a default name and a declaration whose init is Pharmpy’s (modeling/odes.py), scaled from the parent’s seeded estimates:

role	name	init	bounds	η
first peripheral q, v2	Q, V2	Q = CL, V2 = 0.05 · Vc	(0, 10⁶)	add_diagonal only
second peripheral q3, v3	Q3, V3	Q3 = 0.9 · CL (and a first one added alongside it 0.1 · CL), V3 = 0.05 · Vc	(0, 10⁶)	add_diagonal only
ka	KA	1 / (2 · t₁), t₁ the first positive observation time	(0, 10⁶)	add_diagonal only
lagtime	ALAG	t₁ / 2	(0, 10⁶)	absorption_delay, add_diagonal
mtt	MTT	t₁ / 2	(0, 10⁶)	absorption_delay, add_diagonal
n	NTR	the count, FIX; or 2.0 for TRANSITS(N)	(0, 64)	never

A move that changes a chain the parent already has — another count, or fixed ↔︎ estimated — binds n to a fresh parameter (NTR2 when the parent’s is NTR) rather than rewriting the parent’s declaration, and the old one is pruned with its θ; MTT is kept. Names and inits are read off the parent’s own text after its edits and seeding, not the input’s, so a name an earlier step pruned (a three-compartment input narrowed to one, then widened again) is re-declared, and a θ an earlier step added is not declared twice.

A new η has variance 0.01, Pharmpy’s initial estimate. A generated name that collides with a θ or η the model already declares takes a _2 suffix; a parameter the model declares but does not bind (a leftover KA on an IV base) is bound as it is. A template that drops a parameter — two compartments back to one, or first-order absorption to a transit chain, which reads MTT in place of KA — drops its θ and η too; a [covariate_model] relation on a dropped parameter goes with it.

That a generated candidate is the model a pharmacometrician would have typed is pinned the way the [covariate_model] desugar is: the search’s two-compartment lagged candidate and its hand-written twin agree on predict() bit for bit and on the evaluated objective bit for bit (crates/ferx-tools/tests/modelsearch_end_to_end.rs).

The [odes] candidates

Three of the MFL values have no analytic pk template at all: ABSORPTION(ZO), ABSORPTION(WEIBULL) and every ELIMINATION other than FO. They are one family, because they are written the same way — an ode_template NAME(...) line, which generates the standard disposition (ode_template), plus one [odes] line replacing the central equation:

feature	the candidate
ABSORPTION(ZO)	ode_template one_cpt_iv(cl=CL, v=V) + d/dt(central) = zero_order(dur=DUR) - (CL/V) * central
ABSORPTION(WEIBULL)	ode_template one_cpt_iv(cl=CL, v=V) + d/dt(central) = weibull(td=TD, beta=BETA) - (CL/V) * central
ELIMINATION(MM)	the base’s own template + d/dt(central) = KA * depot - ((CL * KM / (KM + central / V)) / V) * central
ELIMINATION(ZO)	the same, with theta TVKM(…) FIX
ELIMINATION(MIX-FO-MM)	… - ((CL + CLMM * KM / (KM + central / V)) / V) * central

Everything else about the candidate is unchanged: the disposition, the peripheral compartments, a transit chain, the covariate model, and the parameters the base already had. Only central changes, so ELIMINATION composes with every absorption and every compartment count — it appears in none of Pharmpy’s incompatible pairs.

Four things are worth knowing before you put one in a space.

The clearance keeps its name

Pharmpy renames CL to CLMM when it switches to Michaelis–Menten; ferx binds the same parameter, under the name the base gave it. That is what carries its initial estimate and its η across the move — a renamed parameter would arrive with neither. MIX-FO-MM is the one case with two clearances, and there the saturable one is a new CLMM starting at half the first-order one, exactly as Pharmpy does it.

They cost an order of magnitude more per fit

A closed form is replaced by numerical integration, and a saturable right-hand side is stiffer again — one FOCEI evaluation of the warfarin Michaelis–Menten candidate below takes 0.064 s against the analytic base’s 0.008 s, and a fit widens that, because every outer iteration integrates the system again with sensitivities. The runner plans candidates of equal cost together, heaviest group first, so an ODE candidate runs with subject-level threads rather than alone on one worker while the rest of the machine idles — but the wall clock is still the wall clock, and the report’s seconds column is the measurement.

A saturable elimination is fitted with more starts

Michaelis–Menten models are the canonical local-minimum case, and a search that ranks a stalled MM fit against a converged first-order one rejects a correct model on the strength of the optimizer. Those candidates — and only those — are given at least 8 starts, whatever [run] retries says; a run configured with more keeps them.

Initial estimates, and a floor Pharmpy does not have

ABSORPTION(ZO) drops the depot, as Pharmpy’s set_zero_order_absorption does: the dose is delivered into central over a modeled duration, so KA and its η are pruned with it. The duration starts at Pharmpy’s 2·MAT with MAT = 2·t_first; the Weibull scale at MAT / Γ(1 + 1/β) with β = 1.5; the Michaelis constant at max(DV)/2, bounded above by 1.5·max(DV), or fixed at min(DV)/100 for zero-order elimination, which is how Pharmpy writes a zero-order elimination too.

The observation range those last two read is floored positive, which Pharmpy does only half of. A Michaelis constant at or below zero is not a model — the saturable term CL·KM/(KM + C) is singular at C = −KM, and at KM = 0 the generated right-hand side evaluates 0/0 when central is zero, which is every subject’s first record — and an additive residual error puts negative observations in perfectly ordinary data. Pharmpy resets a negative min(DV)/100 to 0.01; ferx applies the same fallback to a zero minimum and to a non-positive max(DV), which would otherwise give KM a negative upper bound. The fallback is the value Pharmpy uses when a model carries no dataset at all.

Bioavailability and the lag time survive every move

Including a second one. They are not coordinates of the search, and an ode_template line cannot carry them — it takes the analytic signature, which has no f or lagtime slot — so an ODE candidate states them in [individual_parameters] instead: under a reserved spelling (F, ALAG, LAGTIME) when the model already uses one, and otherwise as the alias lagtime = TLAG that the engine reads the slot off. A later step recovers the binding from there rather than from the template line, so ELIMINATION(MM) followed by PERIPHERALS(1) leaves both exactly where the base put them.

The base must be analytic

A base model whose disposition is already an ode_template line is refused: the search reads its starting coordinates off a pk NAME(...) line and cannot recover an absorption or an elimination from an arbitrary right-hand side. Write the base analytically and let the search generate the ODE candidates.

Ranking, and convergence beside it

Every fitted model — input, base, every candidate — is ranked on [rank] type among the models that pass the [strictness] gate. Under automation a fit that stalled at its initial estimates (#751) or landed in the wrong inner-EBE mode (#864) is not a slow fit, it is a model-selection error, so every candidate gets [run] retries perturbed restarts first, then the gate, then a rank. A model that fails the gate, does not compile or does not fit keeps its row with the reason and no rank:

Base model (base): FO, 0 peripherals — OFV -286.004, bic_mixed -262.985
Algorithm: exhaustive, 4 models fitted

  model    structure                                 npar          OFV    bic_mixed         d rank seconds  status        decision
  base     FO, 0 peripherals                            7     -286.004     -262.985    +0.000    1     0.9  ok            SELECTED
  run1     FO, 1 peripheral                             9     -289.073     -256.653    +6.332    2     2.1  ok
  run3     FO, 1 peripheral, lag                       11     -287.483     -252.760   +10.225    3     3.4  ok
  run2     FO, 0 peripherals, lag                       9     -287.053     -257.427    +5.558    -     1.7  not converged stalled at the initial estimates (#751)

With a cutoff, a candidate is selected only when it beats the base by at least that much on the criterion; the base stays otherwise, and the ranks still show what came closest.

The input model is ranked but never selected. It only has a row of its own when a base had to be derived from it, and that happens exactly when its structure lies outside the space the MFL declares — so selecting it would write a final.ferx the space excluded. Read its row to see what the move onto the space cost; if you want the input’s structure to be selectable, name its value in the category concerned (ABSORPTION([INST,FO]) rather than ABSORPTION(FO)) and it becomes the base. When the base itself fails the gate the candidates are ranked among themselves and the cutoff cannot apply, which the notes say.

What a child is seeded from is never silent. A candidate the runner deduplicated (the same canonical text as an earlier one) is seeded from its representative’s fit, and carries that fit out if it wins; a candidate that produced no fit at all extends with its own initial estimates, and the notes say so; and a --resume row whose cached fits/<hash>.json is missing or unreadable is an error naming the row and the fix (refit without --resume) rather than a search that quietly continues from the file’s initials.

The seconds column is deliberate. The analytic templates cost about the same per fit, but the column is what shows when one does not — and it is the first thing an [odes] candidate would need, should the elimination gap close.

Output

Into --directory (default {search}-modelsearch next to the .ferxsearch file, or [run] cache_dir):

file	contents
models.csv	one row per fitted model: id, parent, layer, path, absorption, peripherals, transits, lagtime, n_parameters, ofv, criterion, d_criterion, rank, converged, passed, failures, error, seconds, selected, continued, reused
final.ferx	the selected model, with its own estimates written into the initial values
final-fit.yaml	its estimates
models/<id>.ferx	every fitted model as it was fitted, so the runner-up can be read and refitted by hand
base/, layer-1/, layer-2/, … (or candidates/ for exhaustive)	one runner directory per layer: candidates.csv, the journal and the cached fits

The per-layer directories are what --resume reads. A search interrupted at layer 3 restarts from the same base, generates the same layers, and refits only what its journals do not already hold; a reused row is marked in the table.

Compared with Pharmpy

The anchor for the enumeration is Pharmpy 2.0.0’s own workflow. tools/pharmpy-modelsearch-anchor/dump.py builds, for nine (base, space, algorithm) cases, the base Pharmpy derives, the candidate moves it filters and every create_candidate task with its feature set and parents — without fitting anything — and a Tier-1 test replays the cases through ferx. The candidate moves agree in Pharmpy’s order, the (feature set, parent feature set) edges agree as multisets through the reduced-stepwise collectors, and the exhaustive combinations agree in order, so run{n} is Pharmpy’s modelsearch_run{n}. Where ferx differs the case is named in the test with its reason, and the test asserts the two sides really do differ:

case	Pharmpy 2.0.0	ferx
a space with no ABSORPTION statement, oral base	fills the category with its default and derives a bolus base	keeps the base’s absorption
LAGTIME([OFF,ON]) on an IV base	a lag time on the bolus is a candidate (the pair table sees only applied features)	never builds a lagged bolus
reduced_stepwise with one duplicate group in a layer	does not collapse it (len(groups) > 1) and runs the full stepwise tree	collapses every group
exhaustive with INST and LAGTIME(ON) in the space	builds INST + LAGTIME(ON)	applies the pair table to combinations too

Transit cases are absent from the anchor on purpose: Pharmpy’s base derivation raises on a space that lists TRANSITS(0, NODEPOT) against a first-order model, and its model-feature reader counts compartments rather than transits, so a base derived with TRANSITS(1) re-applies its own chain as a candidate. The README in the anchor directory records both; the transit rules of the pair table are pinned by unit tests transcribed from _is_allowed.

Compared with NONMEM

The anchor for the numbers is NONMEM 7.5.1 on the warfarin dataset (data/warfarin.csv, 10 subjects), FOCEI, over the space PERIPHERALS(0..1); LAGTIME([OFF,ON]) run exhaustively. The base was fitted in NONMEM (base.ctl); the three candidates were written by hand the way the search writes them — seeded from the base’s .ext, with Q = CL, V2 = 0.05 · V, ALAG = 0.25 and an η of 0.01 on the lag — and fitted (lag.ctl, p1.ctl, p1_lag.ctl). The control streams, .lst and .ext files are under crates/ferx-tools/tests/nonmem/modelsearch_anchor/, and tests/modelsearch_nonmem_anchor.rs replays the comparison.

model	NONMEM, own run	NONMEM at ferx’s estimates	ferx	mixed BIC	rank
base, one_cpt_oral	−286.0042 (successful)	—	−286.0042	−267.488	1
LAGTIME(ON)	−287.0526 (terminated)	−287.6292	−287.6292	−264.508	2
PERIPHERALS(1)	−289.0727 (terminated)	−289.0986	−289.0986	−261.182	3
LAGTIME(ON); PERIPHERALS(1)	−287.4830 (terminated)	−287.6292	−287.6292	−255.107	4

NONMEM’s own minimiser reports MINIMIZATION TERMINATED on all three candidates: warfarin’s first sample is at 0.5 h, so the lag is barely identified and its η collapses (ω → 10⁻⁶), and the two-compartment lag model’s peripheral collapses with it. On those flat directions NONMEM stopped short, so the comparison is made at the same point: *_eval.ctl evaluates NONMEM’s objective (MAXEVAL=0) at ferx’s estimates, and the two engines agree to 1.2 × 10⁻⁹ (lag), 3.7 × 10⁻⁹ (two-compartment) and 5.4 × 10⁻⁵ (two-compartment lag, on a corner with Q = 3069 and V2 = 9 × 10⁻⁶). And lag_refit.ctl re-minimises the lag model from ferx’s estimates: NONMEM declares MINIMIZATION SUCCESSFUL there, at −287.62917388 against ferx’s −287.62917388 — ferx’s optimum is one NONMEM confirms, and NONMEM’s own run had stopped above it.

The ranking on the mixed BIC is the same under either NONMEM reading — base, lag, two compartments, both — and the base is selected on both sides: no candidate earns its extra parameters on ten subjects. The smallest BIC gap the ranking depends on (2.98, base against lag) is four orders of magnitude above the worst OFV disagreement, so it is not a ranking decided by numbers finer than two engines agree on. The tolerance in the test is measured, not assumed: 10⁻³, twenty times the worst realised error.

The [odes] candidates, anchored on an evaluation

The generated [odes] equation is a second numerical object, and NONMEM can state it exactly ($DES with a CLMM·KM/(KM + C) clearance), so it has an anchor of its own — mm_base.ctl / mm_eval.ctl beside the four above. Here both sides evaluate (MAXEVAL=0 / maxiter = 0) at one stated parameter vector, the base’s own θ plus KM = max(DV)/2 = 6.9409:

model	ferx	NONMEM 7.5.1	|Δ|
first-order elimination (analytic, the control)	−182.66670809	−182.66670808	9.7 × 10⁻⁹
ELIMINATION(MM) (ode_template + one override)	929.63056169	929.63055241	9.3 × 10⁻⁶

An evaluation rather than a fit, on purpose: the object under test is the equation, and comparing two minimisations of a Michaelis–Menten model would compare two optimizers on the canonical multimodal surface instead. The first-order row is the control — it says the agreement on the second row belongs to the saturable term and is not an accident of the dataset — and the 1112-unit gap between the two rows is what makes the anchor unable to pass on a candidate whose override the engine ignored. Both sides integrate to a stated tolerance (TOL=9 ATOL=12 against ode_reltol = 1e-10, ode_abstol = 1e-12), so the 9.3 × 10⁻⁶ residual is solver error rather than model error; at ferx’s default ode_reltol = 1e-4 the same comparison is 0.017 out, which is the FOCEI objective amplifying integration error and not a disagreement about the model.

Deviations from Pharmpy, stated

• The pair table is applied to the whole structure, and to the exhaustive combinations. Pharmpy checks its incompatible pairs against the features a path applied, so a lag time on a bolus base, or a transit chain on a base that already has a lag, is a candidate; and exhaustive skips the check entirely. ferx never generates a model the table calls meaningless, whichever way it would be reached.
• A category the space does not name keeps the base’s value. Pharmpy fills it with the category’s default (INST, no transits, no peripherals, no lag) and moves the base there.
• reduced_stepwise collapses every duplicate group. Pharmpy 2.0.0 collapses only when a layer holds two or more such groups.
• The collapse prefers a model that passed the gate. Pharmpy keeps the lowest OFV among the group’s models that have any result; ferx keeps the lowest OFV among those that passed strictness, and only when none did the lowest of the rest.
• TRANSITS(n) is a FIXed θ, not a literal binding, so the count keeps a name the estimates file and the ODE twin can see.
• fullblock is not offered as an IIV strategy; that move belongs to the variability search.
• The input model is ranked but not selectable. Pharmpy ranks the base and the candidates; ferx additionally shows the input’s row when a base had to be derived, and excludes that row from selection so the space binds the result.

From Rust

use ferx_tools::modelsearch::{run_modelsearch, ModelsearchRun};
use ferx_tools::search::SearchConfig;

let cfg = SearchConfig::load("warfarin.ferxsearch")?;
let base = cfg.load_base()?;
let result = run_modelsearch(&cfg, &base, ModelsearchRun {
    dir: Some("warfarin-modelsearch".into()),
    ..ModelsearchRun::default()
})?;
for row in result.ranked() {
    println!("{} {:?} {:.3} rank {:?} {:.1}s",
        row.id, row.structure, row.criterion, row.rank, row.seconds);
}
println!("{}", result.final_model.render());
let runner_up = &result.models["run2"];   // every fitted model's text, by id

ModelsearchRun also takes a threads override, a CancelFlag that stops the search between layers with what it has, and a progress callback (ModelsearchEvent).