[mhovd]

Stacked on #333 (chore/pharmsol-main). Retarget to main once that merges. Until then CI won't run, because rust.yml only triggers for PRs into main.

Why

NonParametricRunner had ~20 methods. Most were accessors (theta, psi, cycle, set_status, log_cycle_state, into_result, …) that each of the 4 algorithms reimplemented almost identically. NPMAP and NCNPAG stubbed out condensation, optimizations and expansion, then overrode fit to skip the cycle loop. Because the controller path didn't use that override, the two paths behaved differently. Nothing kept θ, Ψ and the weights in agreement, and every algorithm had to remember the cycle cap, the stop file and logging on its own.

Design

Solution      theta + psi + weights + log-likelihood, private fields. Built only by
              solving for optimal weights (IPM) or from given weights, and replaced
              as a whole, so the four always describe the same fit.

FitState<E>   the problem (equation, data, error models, prior) + the current Solution
              + cycle bookkeeping, all private. Runners change it only through
              set_solution / adopt_error_models. Shared steps: estimate, condense,
              compute_psi.

Runner<E>     trait Runner<E> {
                  fn step(&mut self, &mut FitState<E>) -> Result<CycleOutcome>;
              }                                   // CycleOutcome = Continue | Converged
              Each algorithm only holds its own state:
                Npag   { config, ranges, eps, last_fine_objective, error_optimizer }
                Npod   { config, error_optimizer }
                Npmap
                Ncnpag { config }

Fit<E>        FitState + Box<dyn Runner> + max_cycles. Owns everything around a cycle:
              the counter, tracing span, cycle cap, stop file, summary logging and
              cycle log. Hands its state over as the result without copying.

An NPAG cycle now reads top to bottom:

fn step(&mut self, state: &mut FitState<E>) -> Result<CycleOutcome> {
    let mut theta = state.solution().theta().clone();
    if !state.is_first_cycle() {
        adaptative_grid(&mut theta, self.eps, &self.ranges, self.config.grid_tolerance)?;
    }
    let estimate = state.estimate(theta, self.config.progress && state.is_first_cycle())?;
    let solution = state.condense(estimate, self.config.prune_threshold, self.config.qr_tolerance)?;
    state.set_solution(solution);
    self.error_optimizer.optimize(state)?;
    Ok(if self.converged(state) { CycleOutcome::Converged } else { CycleOutcome::Continue })
}

Other changes:

• Condensation (lambda pruning + QR) was duplicated in NPAG and NPOD. It now lives in FitState::condense.
• Error-model optimizer: ErrorModelOptimizer owns its adaptive steps and works on candidate Solutions. It previously took 9 &mut arguments.
• Stale weights: the optimizer has always kept the old weights when it adopts a better Ψ. That is now a single named call, Solution::with_stale_weights, documented as the one deliberate break of the invariant (see the open question below).
• Algorithm impls: one blanket impl<A: Into<NonParametricAlgorithm>> Algorithm<E, NonParametric> for A replaces the 5 hand-written impls. The parametric side gets the same treatment.
• Unwraps: the NPOD candidate optimizer and the zero-probability diagnostics no longer unwrap(). Errors are now propagated or logged instead.

Behaviour changes

NPAG and NPOD results are bit-identical to the base branch, including with error-model optimisation active (see Verification).

NPMAP and NCNPAG change:

1. objf is now −2LL. It used to be the log-likelihood (and cycles.csv labelled it neg2ll).
2. They report cycles() == 1, not 0. tests/onecomp.rs is updated for this.
3. fit_controller(NpmapConfig) now writes the cycle log entry. It used to be empty on that path.
4. They now log the cycle summary, like NPAG and NPOD.

Stop file:

• The stop file is still honoured by fit_with and fit_controller, and a stale one is still removed at start. This is now documented on the nonparametric module.
• Fit::new itself no longer touches the filesystem. The public entry points opt in with watch_stop_file().
• NCNPAG's nested per-point NPAG runs no longer delete or react to the stop file. Previously each nested run deleted it, so a stop created during NCNPAG stopped at most one point's refinement.

Minor:

• On the first cycle, NPCycle::delta_objf is inf instead of ~1e30.
• NPOD now runs the zero-probability check after condensation, as NPAG already did.
• The "converged" log line is the same for every algorithm.

API changes

• Removed: NonParametricRunner and the public structs NPAG, NPOD, NPMAP and NCNPAG. Their constructors were pub(crate), so nothing outside the crate could build one.
• NpodConfig gains objective_tolerance, min_distance, prune_threshold and qr_tolerance, with setters. These were hard-coded constants; the defaults are unchanged. #[serde(default)] keeps older serialised configs loading, and a test covers this.
• FitController gains weights(), log_likelihood() and objf(); before this there was no way to read the weights while a fit ran. likelihood() and n2ll() are removed. objf is now the name everywhere in the Rust API (result, cycle log, controller, FitResult) and always means −2LL. The neg2ll CSV column is unchanged.
• Looser bounds: FitController<E> and fit_controller go from E: Equation + Send + 'static to E: Equation. That is no looser in practice, since Equation: Clone + Send + Sync + 'static.

Public API (third commit, feat!)

These changes are breaking, and nothing is kept as a deprecated alias.

1. A stopped fit stays stopped. FitController::step used to run another full cycle even after the fit had stopped, despite its docs saying further steps changed nothing. It also overwrote a request_stop(), which is why the docs said to call request_stop "just before into_result". Now:
   • once a fit has stopped, step does nothing and returns the final status;
   • request_stop is final;
   • a fit that already stopped keeps its reason.
2. NonParametricResult::status() and stop_reason(). Telling MaxCycles, StopFile and Aborted apart no longer needs the cycle log.
3. FitController::run_until(|fit| -> bool) replaces fit_with_observer. FitObserver and CycleFlow are removed. The closure runs after every cycle, the last one included, so it also works for progress reporting:

   let result = problem
       .fit_controller(NpagConfig::new())?
       .run_until(|fit| {
           println!("cycle {} | -2LL {:.4}", fit.cycle(), fit.objf());
           fit.cycle() >= 50
       })?;

4. get_theta() → theta().
5. PredictionGrid replaces the positional idelta, tad floats in predictions, write_predictions, write_json, write_outputs and NPPredictions::calculate:

   result.write_outputs("outputs/run", PredictionGrid::new().idelta(0.5).tad(12.0))?;

   NPPredictions::calculate also stops cloning the whole dataset before expanding it.
6. ZeroProbabilityError. Subjects with zero probability now fail with a typed error that names them, which you can recover with err.downcast_ref::<ZeroProbabilityError>():
   • With log-scaled Ψ, this failure happens inside calculate_psi. It used to surface there as a generic "Each subject must have at least one finite log-likelihood" without naming anyone, and the per-subject likelihood diagnostics in FitState were effectively unreachable.
   • Both now happen on that path.
7. One way to pick an algorithm. Docs, examples, tests and benches use the config structs (NpagConfig::new()). NonParametricAlgorithm is documented for choosing at runtime.

Not changed on purpose:

• EstimationProblem::nonparametric(equation, data, prior, error_models). All four arguments have distinct types, so the compiler already rejects a mix-up. A builder would replace that with a runtime "missing prior" error.
• NonParametricAlgorithm::npag() etc. These are kept, but nothing in the repo uses them any more.

Verification

• Fingerprint. A throwaway harness (not committed) ran every algorithm on bimodal_ke before and after each commit: NPAG and NPOD, each with λ = 0 and λ = 2 (error-model optimisation active), NPMAP, NCNPAG with 50 and 0 cycles, 20 controller-stepped NPAG cycles, the controller NPMAP path, and a callback that stops NPOD after 5 cycles. Cycles, −2LL, nspp, weights, θ, cycle log and error-model factors are bit-identical. The only exceptions are the NPMAP/NCNPAG changes in the first commit; their objf is exactly −2× the old value. The last commit is bit-identical to the one before it for all four algorithms.
• Tests: cargo test --all-features: 95 passed. The new tests cover:
  • NPMAP objf = −2LL;
  • Solution::optimal agreeing with the marginal likelihood;
  • controller θ, weights and Ψ staying aligned between steps;
  • old NpodConfig JSON still loading;
  • stopped fits staying stopped, and request_stop keeping an earlier reason;
  • run_until (every cycle seen; both stop paths);
  • ZeroProbabilityError detection and downcasting.
• Clippy: cargo clippy --all-targets --all-features reports 1 warning, which predates this branch. The other old one, a redundant + Sync bound, is fixed.
• Examples: bimodal_ke gives the same result as chore: use pharmsol from main #333 (320 cycles, −2LL −429.2039). Its five output files (theta, posterior, pred, covs, cycles) are byte-identical before and after the third commit. theophylline, two_eq_lag, neely, chain_algorithms and bestdose (with --features dsl) all run cleanly.
• Fingerprint across commits: the third commit is bit-identical to the second. The harness's observer run was ported to run_until with the same stopping rule.

Performance

An earlier version of this description reported −8.9 % for NPAG. That number was wrong. The base had been built in a separate worktree, and since the hot ODE-solver code is instantiated inside the bench crate, that comparison was between two differently compiled solvers. The numbers below come from bench binaries built in the same checkout and run interleaved. They were measured with a since-dropped commit included (see below), which had no measurable effect.

Bench	base	this PR	change
bimodal_ke_npag (6 interleaved rounds, order rotated)	2.88 s	2.93–3.00 s	none measurable (see below)
bimodal_ke_npod (2 rounds)	0.997 s	0.946 s	~−5 %
bimodal_ke_postprob (2 rounds)	461 ms	443 ms	~−4 %

• NPAG: single-threaded, a full fit after the second commit retires 505.5 vs 504.9 billion instructions compared with the first commit (+0.13 %), with the same CPU cycles and wall time. Multi-threaded runs of the same binaries, alternating order, agree within ~1 %. Criterion runs on this laptop drift by ±5 %, which is the noise floor.
• NPOD and postprob: the gains come from turning the result into a consuming move instead of cloning data, Ψ and θ, and from dropping redundant clones.
• Dropped commit: a commit that ran the error-model optimiser's "up" and "down" trials concurrently (rayon::join) has been removed. It showed no measurable gain, because the Ψ computation already keeps every core busy, so the code stays sequential.

Open question (not changed here)

When the error-model optimiser adopts a better Ψ, it keeps the weights from before (Solution::with_stale_weights). From then on, Ψ and the weights no longer describe the same fit, and those weights flow into the cycle log, the result, NPAG's Σln(Ψ·w) convergence test and NPOD's expansion.

Dropping that call adopts the IPM weights instead. That leaves NPAG unaffected, but npod with λ = 2 then converges in 15 cycles at −2LL −431.75, where it currently takes 18 cycles to reach −435.40. Because the change alters numerical results, it is left for a separate PR.

Not addressed

• Psi::nspp() returns the number of subjects (rows) and nsub() the number of support points: the two are swapped. That's a public-API fix for a separate PR.

[github-actions]

Bencher Report

Project	PMcore
Branch	refactor/algorithm-trait
Testbed	mhovd-pgx

Click to view all benchmark results

Benchmark	Latency	Benchmark Result milliseconds (ms) (Result Δ%)	Upper Boundary milliseconds (ms) (Limit %)
bimodal_ke_npag	📈 view plot 🚷 view threshold	3,525.10 ms (-18.59%) Baseline: 4,330.08 ms	6,428.47 ms (54.84%)
bimodal_ke_npod	📈 view plot 🚷 view threshold	801.45 ms (-38.39%) Baseline: 1,300.76 ms	2,053.72 ms (39.02%)
bimodal_ke_postprob	📈 view plot 🚷 view threshold	345.02 ms (-11.30%) Baseline: 388.99 ms	831.47 ms (41.50%)

🐰 View full continuous benchmarking report in Bencher