diff --git a/.github/workflows/base_benchmarks.yml b/.github/workflows/base_benchmarks.yml index 0353e30b..9b31d3d4 100644 --- a/.github/workflows/base_benchmarks.yml +++ b/.github/workflows/base_benchmarks.yml @@ -28,4 +28,4 @@ jobs: --err \ --adapter rust_criterion \ --github-actions '${{ secrets.GITHUB_TOKEN }}' \ - cargo bench + cargo bench --all-features diff --git a/.github/workflows/pr_benchmarks.yml b/.github/workflows/pr_benchmarks.yml index e71bb67b..ce40f86b 100644 --- a/.github/workflows/pr_benchmarks.yml +++ b/.github/workflows/pr_benchmarks.yml @@ -29,4 +29,4 @@ jobs: --err \ --adapter rust_criterion \ --github-actions '${{ secrets.GITHUB_TOKEN }}' \ - cargo bench + cargo bench --all-features diff --git a/Cargo.toml b/Cargo.toml index 18698813..faa35d37 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -6,6 +6,7 @@ resolver = "2" name = "pharmsol" version = "0.26.1" edition = "2021" +autobenches = false authors = ["Julián D. Otálvaro ", "Markus Hovd"] description = "Rust library for solving analytic and ode-defined pharmacometric models." license = "GPL-3.0" @@ -14,7 +15,13 @@ documentation = "https://lapkb.github.io/pharmsol/" [features] default = [] dsl-core = [] -dsl-jit = ["dsl-core", "dep:cranelift", "dep:cranelift-jit", "dep:cranelift-module", "dep:cranelift-native"] +dsl-jit = [ + "dsl-core", + "dep:cranelift", + "dep:cranelift-jit", + "dep:cranelift-module", + "dep:cranelift-native", +] dsl-aot = ["dsl-core"] dsl-aot-load = ["dsl-core", "dep:libloading"] dsl-wasm-compile = ["dsl-core", "dep:wasm-encoder"] @@ -61,25 +68,10 @@ webbrowser = "1.2.1" bench = false [[bench]] -name = "performance" +name = "native_matrix" harness = false [[bench]] -name = "data" -harness = false - -[[bench]] -name = "ode" -harness = false - -[[bench]] -name = "nca" -harness = false - -[[bench]] -name = "runtime_matrix" -harness = false - -[[bench]] -name = "likelihood_matrix" +name = "dsl_matrix" harness = false +required-features = ["dsl-jit", "dsl-aot", "dsl-aot-load", "dsl-wasm"] diff --git a/benches/common/mod.rs b/benches/common/mod.rs new file mode 100644 index 00000000..997e6e9c --- /dev/null +++ b/benches/common/mod.rs @@ -0,0 +1,650 @@ +//! Shared bench fixtures (workloads, subjects, params, error models, model factories) +//! used by both `native_matrix.rs` and `dsl_matrix.rs`. Every backend measures the +//! same model + subject + params so only the engine varies between cells. +//! +//! Workloads: +//! - [`Workload::Short`]: 1-cpt, 100 mg PO at t=0, 9 obs over 12 h. +//! - [`Workload::Repeat`]: 2-cpt, 100 mg IV q12h × 10, 14 obs over 120 h. +//! +//! `#![allow(dead_code)]` since each bench binary compiles separately and may +//! not use every helper. + +#![allow(dead_code)] + +use ndarray::Array2; +use pharmsol::prelude::*; +use pharmsol::simulator::equation::analytical::{ + one_compartment_with_absorption, two_compartments, +}; +use pharmsol::{ + equation::{self, Route}, + Analytical, ResidualErrorModel, ResidualErrorModels, ODE, SDE, +}; + +/// `ModelMetadata` for handwritten factories so route/output labels resolve like the macro/DSL paths. +fn model_metadata(workload: Workload, kind: SolverKind) -> equation::ModelMetadata { + let name = match (workload, kind) { + (Workload::Short, SolverKind::Ode) => "bench_one_cpt_po_ode", + (Workload::Short, SolverKind::Analytical) => "bench_one_cpt_po_analytical", + (Workload::Short, SolverKind::Sde) => "bench_one_cpt_po_sde", + (Workload::Repeat, SolverKind::Ode) => "bench_two_cpt_iv_ode", + (Workload::Repeat, SolverKind::Analytical) => "bench_two_cpt_iv_analytical", + (Workload::Repeat, SolverKind::Sde) => "bench_two_cpt_iv_sde", + }; + match workload { + Workload::Short => { + let params: &[&str] = match kind { + SolverKind::Sde => &["ka", "ke", "v", "sigma_ke"], + _ => &["ka", "ke", "v"], + }; + equation::metadata::new(name) + .parameters(params.iter().copied()) + .states(["depot", "central"]) + .outputs(["plasma"]) + .route( + Route::bolus("po") + .to_state("depot") + .inject_input_to_destination(), + ) + } + Workload::Repeat => { + let params: &[&str] = match kind { + SolverKind::Sde => &["ke", "kcp", "kpc", "v", "sigma_ke"], + _ => &["ke", "kcp", "kpc", "v"], + }; + equation::metadata::new(name) + .parameters(params.iter().copied()) + .states(["central", "peripheral"]) + .outputs(["plasma"]) + .route( + Route::bolus("iv") + .to_state("central") + .inject_input_to_destination(), + ) + } + } +} + +/// SDE particle count. Kept low so SDE wall-clock stays in the same ballpark as ODE/Analytical. +pub const SDE_PARTICLES: usize = 16; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Workload { + /// 1-cpt, 100 mg PO, 12 h. + Short, + /// 2-cpt, 100 mg IV q12h × 10, 120 h. + Repeat, +} + +impl Workload { + pub fn label(self) -> &'static str { + match self { + Self::Short => "1cpt-12h-po", + Self::Repeat => "2cpt-120h-q12h", + } + } + + pub fn all() -> [Workload; 2] { + [Workload::Short, Workload::Repeat] + } +} + +/// Solver kinds. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum SolverKind { + Ode, + Analytical, + Sde, +} + +impl SolverKind { + pub fn label(self) -> &'static str { + match self { + Self::Ode => "ode", + Self::Analytical => "analytical", + Self::Sde => "sde", + } + } + + // SDE bench cells temporarily disabled — too slow for the current matrix. + pub fn all() -> [SolverKind; 2] { + [SolverKind::Ode, SolverKind::Analytical] + } +} + +// ───────────────────────────── Subjects ────────────────────────────── + +/// 9 sampling points for the short workload. +const SHORT_TIMES: &[f64] = &[0.25, 0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0]; + +/// Synthetic plasma concentrations — only need to exercise the likelihood path, not match any model. +const SHORT_OBS: &[f64] = &[0.50, 0.90, 1.60, 2.40, 2.10, 1.50, 1.05, 0.72, 0.48]; + +/// 14 sampling points for the repeat workload. +const REPEAT_TIMES: &[f64] = &[ + 0.5, 2.0, 6.0, 10.0, 14.0, 24.0, 36.0, 48.0, 60.0, 72.0, 84.0, 96.0, 108.0, 120.0, +]; + +const REPEAT_OBS: &[f64] = &[ + 1.80, 1.45, 1.10, 0.90, 1.30, 1.60, 1.55, 1.50, 1.48, 1.45, 1.43, 1.42, 1.41, 0.95, +]; + +/// Subject with `missing_observation` slots — used for prediction benches. +pub fn subject_for_predictions(workload: Workload) -> Subject { + let id = format!("{}-pred", workload.label()); + let mut builder = Subject::builder(id); + match workload { + Workload::Short => { + builder = builder.bolus(0.0, 100.0, "po"); + for &t in SHORT_TIMES { + builder = builder.missing_observation(t, "plasma"); + } + } + Workload::Repeat => { + for dose in 0..10 { + let t = dose as f64 * 12.0; + builder = builder.bolus(t, 100.0, "iv"); + } + for &t in REPEAT_TIMES { + builder = builder.missing_observation(t, "plasma"); + } + } + } + builder.build() +} + +/// Subject with real observations — used for log-likelihood / matrix benches. +pub fn subject_for_likelihood(workload: Workload) -> Subject { + let id = format!("{}-lik", workload.label()); + let mut builder = Subject::builder(id); + match workload { + Workload::Short => { + builder = builder.bolus(0.0, 100.0, "po"); + for (&t, &y) in SHORT_TIMES.iter().zip(SHORT_OBS.iter()) { + builder = builder.observation(t, y, "plasma"); + } + } + Workload::Repeat => { + for dose in 0..10 { + let t = dose as f64 * 12.0; + builder = builder.bolus(t, 100.0, "iv"); + } + for (&t, &y) in REPEAT_TIMES.iter().zip(REPEAT_OBS.iter()) { + builder = builder.observation(t, y, "plasma"); + } + } + } + builder.build() +} + +/// Population dataset for `log_likelihood_matrix`: `n` subjects, same template with small obs perturbations. +pub fn matrix_data(workload: Workload, n: usize) -> Data { + let subjects = (0..n) + .map(|i| { + let offset = i as f64 * 0.01; + let id = format!("{}-{i:03}", workload.label()); + let mut builder = Subject::builder(id); + match workload { + Workload::Short => { + builder = builder.bolus(0.0, 100.0, "po"); + for (&t, &y) in SHORT_TIMES.iter().zip(SHORT_OBS.iter()) { + builder = builder.observation(t, y + offset, "plasma"); + } + } + Workload::Repeat => { + for dose in 0..10 { + let t = dose as f64 * 12.0; + builder = builder.bolus(t, 100.0, "iv"); + } + for (&t, &y) in REPEAT_TIMES.iter().zip(REPEAT_OBS.iter()) { + builder = builder.observation(t, y + offset, "plasma"); + } + } + } + builder.build() + }) + .collect(); + Data::new(subjects) +} + +// ───────────────────────────── Parameters ──────────────────────────── + +/// Reference parameter vector per `(workload, kind)`: +/// Short ODE/Analytical `[ka, ke, v]`, Short SDE adds `sigma_ke`; +/// Repeat ODE/Analytical `[ke, kcp, kpc, v]`, Repeat SDE adds `sigma_ke`. +pub fn params(workload: Workload, kind: SolverKind) -> Vec { + match (workload, kind) { + (Workload::Short, SolverKind::Ode | SolverKind::Analytical) => vec![1.0, 0.2, 50.0], + (Workload::Short, SolverKind::Sde) => vec![1.0, 0.2, 50.0, 0.05], + (Workload::Repeat, SolverKind::Ode | SolverKind::Analytical) => { + vec![0.10, 0.05, 0.04, 50.0] + } + (Workload::Repeat, SolverKind::Sde) => vec![0.10, 0.05, 0.04, 50.0, 0.01], + } +} + +/// Support-point grid for `log_likelihood_matrix`, shape `(n, nparams)`. +/// Rows are small perturbations of [`params`]. +pub fn support_points(workload: Workload, kind: SolverKind, n: usize) -> Array2 { + let base = params(workload, kind); + let nparams = base.len(); + Array2::from_shape_fn((n, nparams), |(row, col)| { + let p = base[col]; + let perturbation = (row as f64) * 0.001 * p.abs().max(1e-3); + p + perturbation + }) +} + +// ───────────────────────────── Error models ────────────────────────── + +/// Assay error model for `estimate_log_likelihood` / `log_likelihood_matrix`. +pub fn assay_error_models() -> AssayErrorModels { + AssayErrorModels::new() + .add( + "plasma", + AssayErrorModel::additive(ErrorPoly::new(0.1, 0.1, 0.0, 0.0), 0.0), + ) + .expect("plasma assay error model") +} + +/// Residual error model. Indexes outputs by dense `usize` — assumes a single output at index 0. +pub fn residual_error_models() -> ResidualErrorModels { + ResidualErrorModels::new().add(0, ResidualErrorModel::constant(0.2)) +} + +// ───────────────────────────── Handwritten factories ───────────────── + +pub fn handwritten_ode(workload: Workload) -> ODE { + match workload { + Workload::Short => { + // 1-cpt, oral absorption. Params = [ka, ke, v]. + ODE::new( + |x, p, _t, dx, b, _rateiv, _cov| { + let ka = p[0]; + let ke = p[1]; + dx[0] = -ka * x[0] + b[0]; + dx[1] = ka * x[0] - ke * x[1]; + }, + |_p, _t, _cov| lag! {}, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, p, _t, _cov, y| { + let v = p[2]; + y[0] = x[1] / v; + }, + ) + .with_nstates(2) + .with_ndrugs(1) + .with_nout(1) + .with_metadata(model_metadata(workload, SolverKind::Ode)) + .expect("short ODE metadata validates") + } + Workload::Repeat => { + // 2-cpt IV bolus. Params = [ke, kcp, kpc, v]. + ODE::new( + |x, p, _t, dx, b, _rateiv, _cov| { + let ke = p[0]; + let kcp = p[1]; + let kpc = p[2]; + dx[0] = -ke * x[0] - kcp * x[0] + kpc * x[1] + b[0]; + dx[1] = kcp * x[0] - kpc * x[1]; + }, + |_p, _t, _cov| lag! {}, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, p, _t, _cov, y| { + let v = p[3]; + y[0] = x[0] / v; + }, + ) + .with_nstates(2) + .with_ndrugs(1) + .with_nout(1) + .with_metadata(model_metadata(workload, SolverKind::Ode)) + .expect("repeat ODE metadata validates") + } + } +} + +pub fn handwritten_analytical(workload: Workload) -> Analytical { + match workload { + Workload::Short => Analytical::new( + one_compartment_with_absorption, + |_x, _t, _cov| {}, + |_p, _t, _cov| lag! {}, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, p, _t, _cov, y| { + let v = p[2]; + y[0] = x[1] / v; + }, + ) + .with_nstates(2) + .with_ndrugs(1) + .with_nout(1) + .with_metadata(model_metadata(workload, SolverKind::Analytical)) + .expect("short Analytical metadata validates"), + Workload::Repeat => Analytical::new( + two_compartments, + |_x, _t, _cov| {}, + |_p, _t, _cov| lag! {}, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, p, _t, _cov, y| { + let v = p[3]; + y[0] = x[0] / v; + }, + ) + .with_nstates(2) + .with_ndrugs(1) + .with_nout(1) + .with_metadata(model_metadata(workload, SolverKind::Analytical)) + .expect("repeat Analytical metadata validates"), + } +} + +pub fn handwritten_sde(workload: Workload) -> SDE { + match workload { + Workload::Short => { + // Drift mirrors handwritten_ode short. Diffusion = sigma_ke on central. + // Params = [ka, ke, v, sigma_ke]. + SDE::new( + |x, p, _t, dx, b, _cov| { + let ka = p[0]; + let ke = p[1]; + dx[0] = -ka * x[0] + b[0]; + dx[1] = ka * x[0] - ke * x[1]; + }, + |p, sigma| { + let sigma_ke = p[3]; + sigma[0] = 0.0; + sigma[1] = sigma_ke; + }, + |_p, _t, _cov| lag! {}, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, p, _t, _cov, y| { + let v = p[2]; + y[0] = x[1] / v; + }, + SDE_PARTICLES, + ) + .with_nstates(2) + .with_ndrugs(1) + .with_nout(1) + .with_metadata(model_metadata(workload, SolverKind::Sde)) + .expect("short SDE metadata validates") + } + Workload::Repeat => { + // 2-cpt IV bolus SDE. Params = [ke, kcp, kpc, v, sigma_ke]. + SDE::new( + |x, p, _t, dx, b, _cov| { + let ke = p[0]; + let kcp = p[1]; + let kpc = p[2]; + dx[0] = -ke * x[0] - kcp * x[0] + kpc * x[1] + b[0]; + dx[1] = kcp * x[0] - kpc * x[1]; + }, + |p, sigma| { + let sigma_ke = p[4]; + sigma[0] = sigma_ke; + sigma[1] = 0.0; + }, + |_p, _t, _cov| lag! {}, + |_p, _t, _cov| fa! {}, + |_p, _t, _cov, _x| {}, + |x, p, _t, _cov, y| { + let v = p[3]; + y[0] = x[0] / v; + }, + SDE_PARTICLES, + ) + .with_nstates(2) + .with_ndrugs(1) + .with_nout(1) + .with_metadata(model_metadata(workload, SolverKind::Sde)) + .expect("repeat SDE metadata validates") + } + } +} + +// ───────────────────────────── Macro factories ─────────────────────── + +pub fn macro_ode(workload: Workload) -> ODE { + match workload { + Workload::Short => ode! { + name: "bench_one_cpt_po_ode", + params: [ka, ke, v], + states: [depot, central], + outputs: [plasma], + routes: [ + bolus(po) -> depot, + ], + diffeq: |x, _p, _t, dx, _cov| { + dx[depot] = -ka * x[depot]; + dx[central] = ka * x[depot] - ke * x[central]; + }, + out: |x, _p, _t, _cov, y| { + y[plasma] = x[central] / v; + }, + }, + Workload::Repeat => ode! { + name: "bench_two_cpt_iv_ode", + params: [ke, kcp, kpc, v], + states: [central, peripheral], + outputs: [plasma], + routes: [ + bolus(iv) -> central, + ], + diffeq: |x, _p, _t, dx, _cov| { + dx[central] = -ke * x[central] - kcp * x[central] + kpc * x[peripheral]; + dx[peripheral] = kcp * x[central] - kpc * x[peripheral]; + }, + out: |x, _p, _t, _cov, y| { + y[plasma] = x[central] / v; + }, + }, + } +} + +pub fn macro_analytical(workload: Workload) -> Analytical { + match workload { + Workload::Short => analytical! { + name: "bench_one_cpt_po_analytical", + params: [ka, ke, v], + states: [depot, central], + outputs: [plasma], + routes: [ + bolus(po) -> depot, + ], + structure: one_compartment_with_absorption, + out: |x, _p, _t, _cov, y| { + y[plasma] = x[central] / v; + }, + }, + Workload::Repeat => analytical! { + name: "bench_two_cpt_iv_analytical", + params: [ke, kcp, kpc, v], + states: [central, peripheral], + outputs: [plasma], + routes: [ + bolus(iv) -> central, + ], + structure: two_compartments, + out: |x, _p, _t, _cov, y| { + y[plasma] = x[central] / v; + }, + }, + } +} + +pub fn macro_sde(workload: Workload) -> SDE { + match workload { + Workload::Short => sde! { + name: "bench_one_cpt_po_sde", + params: [ka, ke, v, sigma_ke], + states: [depot, central], + outputs: [plasma], + particles: 16, + routes: [ + bolus(po) -> depot, + ], + drift: |x, _p, _t, dx, _cov| { + dx[depot] = -ka * x[depot]; + dx[central] = ka * x[depot] - ke * x[central]; + }, + diffusion: |_p, sigma| { + sigma[central] = sigma_ke; + }, + out: |x, _p, _t, _cov, y| { + y[plasma] = x[central] / v; + }, + }, + Workload::Repeat => sde! { + name: "bench_two_cpt_iv_sde", + params: [ke, kcp, kpc, v, sigma_ke], + states: [central, peripheral], + outputs: [plasma], + particles: 16, + routes: [ + bolus(iv) -> central, + ], + drift: |x, _p, _t, dx, _cov| { + dx[central] = -ke * x[central] - kcp * x[central] + kpc * x[peripheral]; + dx[peripheral] = kcp * x[central] - kpc * x[peripheral]; + }, + diffusion: |_p, sigma| { + sigma[central] = sigma_ke; + }, + out: |x, _p, _t, _cov, y| { + y[plasma] = x[central] / v; + }, + }, + } +} + +// ───────────────────────────── DSL source strings ──────────────────── + +/// DSL source for `(workload, kind)`. Compiled via `compile_module_source_to_runtime`. +pub fn dsl_source(workload: Workload, kind: SolverKind) -> &'static str { + match (workload, kind) { + (Workload::Short, SolverKind::Ode) => SHORT_ODE_DSL, + (Workload::Short, SolverKind::Analytical) => SHORT_ANALYTICAL_DSL, + (Workload::Short, SolverKind::Sde) => SHORT_SDE_DSL, + (Workload::Repeat, SolverKind::Ode) => REPEAT_ODE_DSL, + (Workload::Repeat, SolverKind::Analytical) => REPEAT_ANALYTICAL_DSL, + (Workload::Repeat, SolverKind::Sde) => REPEAT_SDE_DSL, + } +} + +/// DSL `name = ...` field for `(workload, kind)`. +pub fn dsl_model_name(workload: Workload, kind: SolverKind) -> &'static str { + match (workload, kind) { + (Workload::Short, SolverKind::Ode) => "bench_one_cpt_po_ode", + (Workload::Short, SolverKind::Analytical) => "bench_one_cpt_po_analytical", + (Workload::Short, SolverKind::Sde) => "bench_one_cpt_po_sde", + (Workload::Repeat, SolverKind::Ode) => "bench_two_cpt_iv_ode", + (Workload::Repeat, SolverKind::Analytical) => "bench_two_cpt_iv_analytical", + (Workload::Repeat, SolverKind::Sde) => "bench_two_cpt_iv_sde", + } +} + +const SHORT_ODE_DSL: &str = r#" +name = bench_one_cpt_po_ode +kind = ode + +params = ka, ke, v +states = depot, central +outputs = plasma + +bolus(po) -> depot + +dx(depot) = -ka * depot +dx(central) = ka * depot - ke * central + +out(plasma) = central / v ~ continuous() +"#; + +const SHORT_ANALYTICAL_DSL: &str = r#" +name = bench_one_cpt_po_analytical +kind = analytical + +params = ka, ke, v +states = depot, central +outputs = plasma + +bolus(po) -> depot + +structure = one_compartment_with_absorption + +out(plasma) = central / v ~ continuous() +"#; + +const SHORT_SDE_DSL: &str = r#" +name = bench_one_cpt_po_sde +kind = sde + +params = ka, ke, v, sigma_ke +states = depot, central +particles = 16 +outputs = plasma + +bolus(po) -> depot + +dx(depot) = -ka * depot +dx(central) = ka * depot - ke * central + +noise(central) = sigma_ke + +out(plasma) = central / v ~ continuous() +"#; + +const REPEAT_ODE_DSL: &str = r#" +name = bench_two_cpt_iv_ode +kind = ode + +params = ke, kcp, kpc, v +states = central, peripheral +outputs = plasma + +bolus(iv) -> central + +dx(central) = -ke * central - kcp * central + kpc * peripheral +dx(peripheral) = kcp * central - kpc * peripheral + +out(plasma) = central / v ~ continuous() +"#; + +const REPEAT_ANALYTICAL_DSL: &str = r#" +name = bench_two_cpt_iv_analytical +kind = analytical + +params = ke, kcp, kpc, v +states = central, peripheral +outputs = plasma + +bolus(iv) -> central + +structure = two_compartments + +out(plasma) = central / v ~ continuous() +"#; + +const REPEAT_SDE_DSL: &str = r#" +name = bench_two_cpt_iv_sde +kind = sde + +params = ke, kcp, kpc, v, sigma_ke +states = central, peripheral +particles = 16 +outputs = plasma + +bolus(iv) -> central + +dx(central) = -ke * central - kcp * central + kpc * peripheral +dx(peripheral) = kcp * central - kpc * peripheral + +noise(central) = sigma_ke + +out(plasma) = central / v ~ continuous() +"#; diff --git a/benches/data.rs b/benches/data.rs deleted file mode 100644 index 1202972c..00000000 --- a/benches/data.rs +++ /dev/null @@ -1,247 +0,0 @@ -use criterion::{criterion_group, criterion_main, Criterion}; -use pharmsol::prelude::*; -use pharmsol::SubjectBuilderExt; // Ensure this trait is in scope -use std::hint::black_box; - -fn subject_builder_benchmark(c: &mut Criterion) { - // Simple subject creation - c.bench_function("SubjectBuilder simple", |b| { - b.iter(|| { - let subject = Subject::builder("patient1") - .bolus(0.0, 100.0, 0) - .observation(3.0, 100.0, 0) - .observation(4.0, 200.0, 0) - .observation(5.0, 300.0, 0) - .build(); - black_box(subject); - }) - }); - - // Subject with covariates - c.bench_function("SubjectBuilder with covariates", |b| { - b.iter(|| { - let subject = Subject::builder("patient1") - .bolus(0.0, 100.0, 0) - .observation(3.0, 100.0, 0) - .observation(4.0, 200.0, 0) - .observation(5.0, 300.0, 0) - .observation(12.0, 300.0, 0) - .covariate("age", 0.0, 30.0) - .covariate("weight", 0.0, 70.0) - .build(); - black_box(subject); - }) - }); - - // Multi-occasion subject - c.bench_function("SubjectBuilder multi-occasion", |b| { - b.iter(|| { - let subject = Subject::builder("patient1") - .bolus(0.0, 100.0, 0) - .observation(1.0, 50.0, 0) - .observation(6.0, 25.0, 0) - .covariate("age", 0.0, 30.0) - .reset() - .bolus(24.0, 120.0, 0) - .observation(25.0, 55.0, 0) - .observation(30.0, 30.0, 0) - .covariate("age", 24.0, 30.0) - .build(); - black_box(subject); - }) - }); -} - -fn data_expand_benchmark(c: &mut Criterion) { - // Create datasets for expansion testing - let simple_data = { - let subject = Subject::builder("patient1") - .bolus(0.0, 100.0, 0) - .observation(1.0, 50.0, 0) - .observation(6.0, 25.0, 0) - .observation(12.0, 15.0, 0) - .build(); - Data::new(vec![subject]) - }; - - let complex_data = { - let subjects = (1..=10) - .map(|i| { - Subject::builder(format!("patient{}", i)) - .bolus(0.0, 100.0 * i as f64, 0) - .observation(1.0, 50.0, 0) - .observation(2.0, 45.0, 1) // Different output equation - .observation(6.0, 25.0, 0) - .observation(24.0, 8.0, 1) - .infusion(48.0, 200.0, 0, 2.0) - .observation(50.0, 30.0, 0) - .build() - }) - .collect(); - Data::new(subjects) - }; - - // Test different expansion intervals - c.bench_function("Data expand simple (1h intervals)", |b| { - b.iter(|| { - let expanded = simple_data.expand(1.0, 0.0); - black_box(expanded); - }) - }); - - c.bench_function("Data expand complex (1h intervals)", |b| { - b.iter(|| { - let expanded = complex_data.expand(1.0, 0.0); - black_box(expanded); - }) - }); - - c.bench_function("Data expand with additional time", |b| { - b.iter(|| { - let expanded = complex_data.expand(1.0, 24.0); - black_box(expanded); - }) - }); -} - -fn dose_modification_benchmark(c: &mut Criterion) { - let create_dosing_data = || { - let subjects = (1..=5) - .map(|i| { - Subject::builder(format!("patient{}", i)) - .bolus(0.0, 100.0, 0) - .bolus(12.0, 100.0, 0) - .bolus(24.0, 100.0, 0) - .infusion(48.0, 200.0, 0, 2.0) - .infusion(72.0, 250.0, 0, 3.0) - .observation(1.0, 50.0, 0) - .observation(25.0, 45.0, 0) - .observation(50.0, 30.0, 0) - .build() - }) - .collect(); - Data::new(subjects) - }; - - c.bench_function("Modify all bolus doses", |b| { - b.iter(|| { - let mut data = create_dosing_data(); - for subject in data.iter_mut() { - for occasion in subject.occasions_iter_mut() { - for event in occasion.events_iter_mut() { - if let Event::Bolus(bolus) = event { - bolus.set_amount(bolus.amount() * 2.0); - } - } - } - } - black_box(data); - }) - }); - - c.bench_function("Modify all infusion doses", |b| { - b.iter(|| { - let mut data = create_dosing_data(); - for subject in data.iter_mut() { - for occasion in subject.occasions_iter_mut() { - for event in occasion.events_iter_mut() { - if let Event::Infusion(infusion) = event { - infusion.set_amount(infusion.amount() * 1.5); - infusion.set_duration(infusion.duration() * 0.8); - } - } - } - } - black_box(data); - }) - }); - - c.bench_function("Conditional dose modification", |b| { - b.iter(|| { - let mut data = create_dosing_data(); - for subject in data.iter_mut() { - for occasion in subject.occasions_iter_mut() { - for event in occasion.events_iter_mut() { - match event { - Event::Bolus(bolus) if bolus.time() >= 24.0 => { - bolus.set_amount(bolus.amount() * 1.2); - } - Event::Infusion(infusion) if infusion.time() >= 48.0 => { - infusion.set_amount(infusion.amount() * 0.9); - } - _ => {} - } - } - } - } - black_box(data); - }) - }); -} - -fn data_operations_benchmark(c: &mut Criterion) { - // Combined benchmark for data filtering and large dataset creation - let large_data = { - let subjects = (1..=100) - .map(|i| { - Subject::builder(format!("patient{:03}", i)) - .bolus(0.0, 100.0, 0) - .observation(1.0, 50.0, 0) - .observation(6.0, 25.0, 0) - .observation(12.0, 15.0, 0) - .covariate("age", 0.0, 20.0 + (i as f64 % 50.0)) - .covariate("weight", 0.0, 60.0 + (i as f64 % 40.0)) - .build() - }) - .collect(); - Data::new(subjects) - }; - - c.bench_function("Create large dataset (100 subjects)", |b| { - b.iter(|| { - let subjects: Vec = (1..=100) - .map(|i| { - Subject::builder(format!("patient{:03}", i)) - .bolus(0.0, 100.0, 0) - .observation(0.5, 80.0, 0) - .observation(1.0, 60.0, 0) - .observation(2.0, 45.0, 0) - .observation(4.0, 30.0, 0) - .observation(8.0, 20.0, 0) - .observation(12.0, 15.0, 0) - .observation(24.0, 8.0, 0) - .covariate("age", 0.0, 25.0 + (i as f64 % 50.0)) - .covariate("weight", 0.0, 60.0 + (i as f64 % 40.0)) - .build() - }) - .collect(); - let data = Data::new(subjects); - black_box(data); - }) - }); - - c.bench_function("Filter include subjects", |b| { - b.iter(|| { - let include_ids: Vec = (1..=20).map(|i| format!("patient{:03}", i)).collect(); - let filtered = large_data.filter_include(&include_ids); - black_box(filtered); - }) - }); - - c.bench_function("Filter exclude subjects", |b| { - b.iter(|| { - let exclude_ids: Vec = (80..=100).map(|i| format!("patient{:03}", i)).collect(); - let filtered = large_data.filter_exclude(exclude_ids); - black_box(filtered); - }) - }); -} - -criterion_group!( - benches, - subject_builder_benchmark, - data_expand_benchmark, - dose_modification_benchmark, - data_operations_benchmark -); -criterion_main!(benches); diff --git a/benches/dsl_matrix.rs b/benches/dsl_matrix.rs new file mode 100644 index 00000000..50072072 --- /dev/null +++ b/benches/dsl_matrix.rs @@ -0,0 +1,485 @@ +//! DSL bench matrix (feature-gated): JIT, native AoT, WASM across all workloads + solvers. +//! Mirrors `native_matrix.rs` but compiles models from DSL source. +//! +//! IDs: +//! - `dsl/compile` → `{workload}/{kind}/{backend}` +//! - `dsl/predictions` → `{workload}/{kind}/{backend}/{cache}` +//! - `dsl/log-likelihood` → `{workload}/{kind}/{backend}/{cache}` +//! - `dsl/likelihood-matrix` → `{workload}/{kind}/{backend}` + +use std::hint::black_box; +use std::path::PathBuf; +use std::time::Duration; + +use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion, SamplingMode}; +use tempfile::TempDir; + +use pharmsol::dsl::{ + compile_module_source_to_runtime, CompiledRuntimeModel, NativeAnalyticalModel, + NativeAotCompileOptions, NativeOdeModel, NativeSdeModel, RuntimeCompilationTarget, +}; +use pharmsol::prelude::*; +use pharmsol::Cache; + +mod common; +use common::{ + assay_error_models, dsl_model_name, dsl_source, matrix_data, params, subject_for_likelihood, + subject_for_predictions, support_points, SolverKind, Workload, +}; + +const MATRIX_N_SUBJECTS: usize = 32; +const MATRIX_N_SUPPORT: usize = 64; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +#[allow(dead_code)] // Aot/Wasm temporarily disabled in `Backend::all` +enum Backend { + Jit, + Aot, + Wasm, +} + +impl Backend { + fn label(self) -> &'static str { + match self { + Self::Jit => "dsl-jit", + Self::Aot => "dsl-aot", + Self::Wasm => "dsl-wasm", + } + } + + // AoT and WASM backends temporarily disabled — too slow for the current matrix. + fn all() -> [Backend; 1] { + [Backend::Jit] + } +} + +#[derive(Debug, Clone, Copy)] +enum CacheState { + Hot, + Cold, +} + +impl CacheState { + fn label(self) -> &'static str { + match self { + Self::Hot => "hot", + Self::Cold => "cold", + } + } + + fn all() -> [CacheState; 2] { + [CacheState::Hot, CacheState::Cold] + } +} + +/// One `TempDir` shared across the bench binary; each compile gets a fresh subdir. +struct AotWorkspace { + root: TempDir, + counter: std::cell::Cell, +} + +impl AotWorkspace { + fn new() -> Self { + Self { + root: tempfile::Builder::new() + .prefix("pharmsol-bench-dsl-aot-") + .tempdir() + .expect("create AoT workspace tempdir"), + counter: std::cell::Cell::new(0), + } + } + + fn fresh(&self, stem: &str) -> PathBuf { + let n = self.counter.get(); + self.counter.set(n + 1); + self.root.path().join(format!("{stem}-{n:04}")) + } +} + +/// Compile `(workload, kind)` with `backend` and return the full `CompiledRuntimeModel`. +fn compile_runtime( + workload: Workload, + kind: SolverKind, + backend: Backend, + aot: &AotWorkspace, +) -> CompiledRuntimeModel { + let source = dsl_source(workload, kind); + let name = dsl_model_name(workload, kind); + let target = match backend { + Backend::Jit => RuntimeCompilationTarget::Jit, + Backend::Aot => { + let dir = aot.fresh(&format!("{}-{}", workload.label(), kind.label())); + RuntimeCompilationTarget::NativeAot(NativeAotCompileOptions::new(dir)) + } + Backend::Wasm => RuntimeCompilationTarget::Wasm, + }; + compile_module_source_to_runtime(source, Some(name), target, |_, _| {}) + .unwrap_or_else(|e| panic!("compile {} via {} failed: {e:?}", name, backend.label())) +} + +fn compile_ode(workload: Workload, backend: Backend, aot: &AotWorkspace) -> NativeOdeModel { + match compile_runtime(workload, SolverKind::Ode, backend, aot) { + CompiledRuntimeModel::Ode(model) => model, + other => panic!( + "expected Ode model for {}, got {:?}", + workload.label(), + other.backend() + ), + } +} + +fn compile_analytical( + workload: Workload, + backend: Backend, + aot: &AotWorkspace, +) -> NativeAnalyticalModel { + match compile_runtime(workload, SolverKind::Analytical, backend, aot) { + CompiledRuntimeModel::Analytical(model) => model, + other => panic!( + "expected Analytical model for {}, got {:?}", + workload.label(), + other.backend() + ), + } +} + +fn compile_sde(workload: Workload, backend: Backend, aot: &AotWorkspace) -> NativeSdeModel { + match compile_runtime(workload, SolverKind::Sde, backend, aot) { + CompiledRuntimeModel::Sde(model) => model, + other => panic!( + "expected Sde model for {}, got {:?}", + workload.label(), + other.backend() + ), + } +} + +// ───────────────────────────── compile group ───────────────────────── + +fn compile_group(c: &mut Criterion) { + use std::time::Instant; + + let mut group = c.benchmark_group("dsl/compile"); + group.sampling_mode(SamplingMode::Flat); + group.sample_size(10); + group.measurement_time(Duration::from_secs(5)); + // Each compile leaks an executable mmap (JIT) or runs rustc (AoT). Without + // a cap, a fast JIT compile (~60 µs) lets Criterion request hundreds of + // thousands of iterations per cell and exhausts the runner's executable + // memory pool / `vm.max_map_count`. We hard-cap real iterations per + // Criterion batch to `MAX_ITERS_PER_BATCH` and scale the reported elapsed + // time linearly so per-iteration timings stay accurate. + const MAX_ITERS_PER_BATCH: u64 = 25; + + let aot = AotWorkspace::new(); + + for workload in Workload::all() { + for kind in SolverKind::all() { + for backend in Backend::all() { + let bench_id = BenchmarkId::from_parameter(format!( + "{}/{}/{}", + workload.label(), + kind.label(), + backend.label() + )); + group.bench_function(bench_id, |b| { + b.iter_custom(|iters| { + let actual = iters.min(MAX_ITERS_PER_BATCH).max(1); + let start = Instant::now(); + for _ in 0..actual { + black_box(compile_runtime( + black_box(workload), + black_box(kind), + black_box(backend), + &aot, + )); + } + let elapsed = start.elapsed(); + elapsed.mul_f64(iters as f64 / actual as f64) + }); + }); + } + } + } + + group.finish(); +} + +// ───────────────────────────── predictions group ───────────────────── + +fn predictions_group(c: &mut Criterion) { + let mut group = c.benchmark_group("dsl/predictions"); + group.sampling_mode(SamplingMode::Flat); + + let aot = AotWorkspace::new(); + + for workload in Workload::all() { + let subject = subject_for_predictions(workload); + for kind in SolverKind::all() { + let theta = params(workload, kind); + for backend in Backend::all() { + for cache in CacheState::all() { + let bench_id = BenchmarkId::from_parameter(format!( + "{}/{}/{}/{}", + workload.label(), + kind.label(), + backend.label(), + cache.label() + )); + match kind { + SolverKind::Ode => { + let model = match cache { + CacheState::Hot => compile_ode(workload, backend, &aot), + CacheState::Cold => { + compile_ode(workload, backend, &aot).disable_cache() + } + }; + group.bench_function(bench_id, |b| { + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + }); + } + SolverKind::Analytical => { + let model = match cache { + CacheState::Hot => compile_analytical(workload, backend, &aot), + CacheState::Cold => { + compile_analytical(workload, backend, &aot).disable_cache() + } + }; + group.bench_function(bench_id, |b| { + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + }); + } + SolverKind::Sde => { + let model = match cache { + CacheState::Hot => compile_sde(workload, backend, &aot), + CacheState::Cold => { + compile_sde(workload, backend, &aot).disable_cache() + } + }; + group.bench_function(bench_id, |b| { + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + }); + } + } + } + } + } + } + + group.finish(); +} + +// ───────────────────────────── log-likelihood group ────────────────── + +fn log_likelihood_group(c: &mut Criterion) { + let mut group = c.benchmark_group("dsl/log-likelihood"); + group.sampling_mode(SamplingMode::Flat); + + let aot = AotWorkspace::new(); + let error_models = assay_error_models(); + + for workload in Workload::all() { + let subject = subject_for_likelihood(workload); + for kind in SolverKind::all() { + let theta = params(workload, kind); + for backend in Backend::all() { + for cache in CacheState::all() { + let bench_id = BenchmarkId::from_parameter(format!( + "{}/{}/{}/{}", + workload.label(), + kind.label(), + backend.label(), + cache.label() + )); + match kind { + SolverKind::Ode => { + let model = match cache { + CacheState::Hot => compile_ode(workload, backend, &aot), + CacheState::Cold => { + compile_ode(workload, backend, &aot).disable_cache() + } + }; + group.bench_function(bench_id, |b| { + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + }); + } + SolverKind::Analytical => { + let model = match cache { + CacheState::Hot => compile_analytical(workload, backend, &aot), + CacheState::Cold => { + compile_analytical(workload, backend, &aot).disable_cache() + } + }; + group.bench_function(bench_id, |b| { + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + }); + } + SolverKind::Sde => { + let model = match cache { + CacheState::Hot => compile_sde(workload, backend, &aot), + CacheState::Cold => { + compile_sde(workload, backend, &aot).disable_cache() + } + }; + group.bench_function(bench_id, |b| { + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + }); + } + } + } + } + } + } + + group.finish(); +} + +// ───────────────────────────── likelihood-matrix group ─────────────── + +fn likelihood_matrix_group(c: &mut Criterion) { + use pharmsol::prelude::simulator::log_likelihood_matrix; + + let mut group = c.benchmark_group("dsl/likelihood-matrix"); + group.sampling_mode(SamplingMode::Flat); + group.sample_size(10); + group.measurement_time(Duration::from_secs(20)); + + let aot = AotWorkspace::new(); + let error_models = assay_error_models(); + + for workload in Workload::all() { + let data = matrix_data(workload, MATRIX_N_SUBJECTS); + for kind in SolverKind::all() { + let theta = support_points(workload, kind, MATRIX_N_SUPPORT); + for backend in Backend::all() { + let bench_id = BenchmarkId::from_parameter(format!( + "{}/{}/{}", + workload.label(), + kind.label(), + backend.label() + )); + match kind { + SolverKind::Ode => { + let model = compile_ode(workload, backend, &aot); + group.bench_function(bench_id, |b| { + b.iter(|| { + black_box( + log_likelihood_matrix( + black_box(&model), + black_box(&data), + black_box(&theta), + black_box(&error_models), + false, + ) + .unwrap(), + ) + }); + }); + } + SolverKind::Analytical => { + let model = compile_analytical(workload, backend, &aot); + group.bench_function(bench_id, |b| { + b.iter(|| { + black_box( + log_likelihood_matrix( + black_box(&model), + black_box(&data), + black_box(&theta), + black_box(&error_models), + false, + ) + .unwrap(), + ) + }); + }); + } + SolverKind::Sde => { + let model = compile_sde(workload, backend, &aot); + group.bench_function(bench_id, |b| { + b.iter(|| { + black_box( + log_likelihood_matrix( + black_box(&model), + black_box(&data), + black_box(&theta), + black_box(&error_models), + false, + ) + .unwrap(), + ) + }); + }); + } + } + } + } + } + + group.finish(); +} + +criterion_group!( + dsl_matrix, + compile_group, + predictions_group, + log_likelihood_group, + likelihood_matrix_group +); +criterion_main!(dsl_matrix); diff --git a/benches/likelihood_matrix.rs b/benches/likelihood_matrix.rs deleted file mode 100644 index b8f94d54..00000000 --- a/benches/likelihood_matrix.rs +++ /dev/null @@ -1,247 +0,0 @@ -use criterion::{ - criterion_group, criterion_main, BenchmarkId, Criterion, SamplingMode, Throughput, -}; -use ndarray::Array2; -use pharmsol::prelude::simulator::{log_likelihood_batch, log_likelihood_matrix}; -use pharmsol::prelude::*; -use pharmsol::{Cache, ResidualErrorModel, ResidualErrorModels, ODE}; -use std::hint::black_box; -use std::time::Duration; - -fn example_equation() -> ODE { - equation::ODE::new( - |x, p, _t, dx, _b, _rateiv, _cov| { - fetch_params!(p, ka, ke, _tlag, _v); - dx[0] = -ka * x[0]; - dx[1] = ka * x[0] - ke * x[1]; - }, - |p, _t, _cov| { - fetch_params!(p, _ka, _ke, tlag, _v); - lag! {0=>tlag} - }, - |_p, _t, _cov| fa! {}, - |_p, _t, _cov, _x| {}, - |x, p, _t, _cov, y| { - fetch_params!(p, _ka, _ke, _tlag, v); - y[0] = x[1] / v; - }, - ) - .with_nstates(2) - .with_nout(1) -} - -fn example_data(n_subjects: usize) -> Data { - let subjects = (0..n_subjects) - .map(|index| { - let dose = 100.0 + index as f64; - let offset = index as f64 * 0.02; - Subject::builder(format!("subject_{index}")) - .bolus(0.0, dose, 0) - .observation(0.5, 4.2 + offset, 0) - .observation(1.0, 5.9 + offset, 0) - .observation(2.0, 7.4 + offset, 0) - .observation(4.0, 6.8 + offset, 0) - .observation(8.0, 4.7 + offset, 0) - .observation(12.0, 3.3 + offset, 0) - .build() - }) - .collect(); - - Data::new(subjects) -} - -fn example_subject() -> Subject { - example_data(1) - .get_subject("subject_0") - .expect("example subject exists") - .clone() -} - -fn example_params() -> [f64; 4] { - [0.60, 0.09, 0.05, 20.0] -} - -fn support_points(n_support_points: usize) -> Array2 { - Array2::from_shape_fn((n_support_points, 4), |(row, column)| match column { - 0 => 0.55 + row as f64 * 0.002, - 1 => 0.08 + row as f64 * 0.0004, - 2 => 0.05, - 3 => 18.0 + row as f64 * 0.15, - _ => unreachable!(), - }) -} - -fn batch_parameters(n_subjects: usize) -> Array2 { - Array2::from_shape_fn((n_subjects, 4), |(row, column)| match column { - 0 => 0.60 + row as f64 * 0.001, - 1 => 0.09 + row as f64 * 0.0002, - 2 => 0.05, - 3 => 20.0 + row as f64 * 0.05, - _ => unreachable!(), - }) -} - -fn assay_error_models() -> AssayErrorModels { - AssayErrorModels::new() - .add( - 0, - AssayErrorModel::additive(ErrorPoly::new(0.0, 0.1, 0.0, 0.0), 0.0), - ) - .unwrap() -} - -fn residual_error_models() -> ResidualErrorModels { - ResidualErrorModels::new().add(0, ResidualErrorModel::constant(0.2)) -} - -fn criterion_benchmark(c: &mut Criterion) { - let assay_error_models = assay_error_models(); - let residual_error_models = residual_error_models(); - - let subject = example_subject(); - let params = example_params(); - let equation_cold = example_equation().disable_cache(); - let equation_hot = example_equation(); - let _ = equation_hot - .estimate_predictions(&subject, params.as_slice()) - .unwrap(); - let predictions = equation_cold - .estimate_predictions(&subject, params.as_slice()) - .unwrap(); - - let mut breakdown_group = c.benchmark_group("ode/runtime-breakdown"); - breakdown_group.sample_size(10); - breakdown_group.warm_up_time(Duration::from_millis(250)); - breakdown_group.measurement_time(Duration::from_secs(1)); - breakdown_group.bench_function("predict-cold", |b| { - b.iter(|| { - black_box( - equation_cold - .estimate_predictions(&subject, params.as_slice()) - .unwrap(), - ) - }) - }); - breakdown_group.bench_function("predict-hot", |b| { - b.iter(|| { - black_box( - equation_hot - .estimate_predictions(&subject, params.as_slice()) - .unwrap(), - ) - }) - }); - breakdown_group.bench_function("score-only", |b| { - b.iter(|| black_box(predictions.log_likelihood(&assay_error_models)).unwrap()) - }); - breakdown_group.bench_function("loglik-cold", |b| { - b.iter(|| { - black_box( - equation_cold - .estimate_log_likelihood(&subject, params.as_slice(), &assay_error_models) - .unwrap(), - ) - }) - }); - breakdown_group.bench_function("loglik-hot", |b| { - b.iter(|| { - black_box( - equation_hot - .estimate_log_likelihood(&subject, params.as_slice(), &assay_error_models) - .unwrap(), - ) - }) - }); - breakdown_group.finish(); - - let mut matrix_group = c.benchmark_group("likelihood/matrix"); - matrix_group.sampling_mode(SamplingMode::Flat); - matrix_group.sample_size(10); - matrix_group.warm_up_time(Duration::from_millis(250)); - matrix_group.measurement_time(Duration::from_secs(2)); - - for (n_subjects, n_support_points) in [(16usize, 32usize), (64usize, 128usize)] { - let case = (example_data(n_subjects), support_points(n_support_points)); - let equation_cold = example_equation().disable_cache(); - let equation_hot = example_equation(); - let _ = log_likelihood_matrix(&equation_hot, &case.0, &case.1, &assay_error_models, false) - .unwrap(); - - matrix_group.throughput(Throughput::Elements((n_subjects * n_support_points) as u64)); - matrix_group.bench_with_input( - BenchmarkId::new("cold", format!("{n_subjects}x{n_support_points}")), - &case, - |b, (data, theta)| { - b.iter(|| { - black_box(log_likelihood_matrix( - &equation_cold, - data, - theta, - &assay_error_models, - false, - )) - .unwrap() - }) - }, - ); - matrix_group.bench_with_input( - BenchmarkId::new("hot", format!("{n_subjects}x{n_support_points}")), - &case, - |b, (data, theta)| { - b.iter(|| { - black_box(log_likelihood_matrix( - &equation_hot, - data, - theta, - &assay_error_models, - false, - )) - .unwrap() - }) - }, - ); - } - matrix_group.finish(); - - let mut batch_group = c.benchmark_group("likelihood/batch"); - batch_group.sample_size(10); - batch_group.warm_up_time(Duration::from_millis(250)); - batch_group.measurement_time(Duration::from_secs(1)); - - for n_subjects in [16usize, 64usize, 256usize] { - let data = example_data(n_subjects); - let parameters = batch_parameters(n_subjects); - let equation_cold = example_equation().disable_cache(); - let equation_hot = example_equation(); - let _ = log_likelihood_batch(&equation_hot, &data, ¶meters, &residual_error_models) - .unwrap(); - - batch_group.throughput(Throughput::Elements(n_subjects as u64)); - batch_group.bench_with_input(BenchmarkId::new("cold", n_subjects), &data, |b, data| { - b.iter(|| { - black_box(log_likelihood_batch( - &equation_cold, - data, - ¶meters, - &residual_error_models, - )) - .unwrap() - }) - }); - batch_group.bench_with_input(BenchmarkId::new("hot", n_subjects), &data, |b, data| { - b.iter(|| { - black_box(log_likelihood_batch( - &equation_hot, - data, - ¶meters, - &residual_error_models, - )) - .unwrap() - }) - }); - } - batch_group.finish(); -} - -criterion_group!(benches, criterion_benchmark); -criterion_main!(benches); diff --git a/benches/native_matrix.rs b/benches/native_matrix.rs new file mode 100644 index 00000000..4ad12e39 --- /dev/null +++ b/benches/native_matrix.rs @@ -0,0 +1,570 @@ +//! Native bench matrix: handwritten + macro models, all solvers, both workloads. +//! Default features only — this is the CI regression signal. DSL backends are in `dsl_matrix.rs`. +//! +//! IDs: +//! - `native/predictions` → `{workload}/{solver}/{authoring}/{cache}` +//! - `native/log-likelihood` → `{workload}/{solver}/{authoring}/{cache}` +//! - `native/likelihood-matrix` → `{workload}/{solver}/{authoring}` + +use std::hint::black_box; +use std::time::Duration; + +use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion, SamplingMode}; +use pharmsol::prelude::*; +use pharmsol::Cache; + +mod common; +use common::{ + assay_error_models, handwritten_analytical, handwritten_ode, handwritten_sde, macro_analytical, + macro_ode, macro_sde, matrix_data, params, subject_for_likelihood, subject_for_predictions, + support_points, SolverKind, Workload, +}; + +const MATRIX_N_SUBJECTS: usize = 32; +const MATRIX_N_SUPPORT: usize = 64; + +#[derive(Debug, Clone, Copy)] +enum Authoring { + Handwritten, + Macro, +} + +impl Authoring { + fn label(self) -> &'static str { + match self { + Self::Handwritten => "handwritten", + Self::Macro => "macro", + } + } + + fn all() -> [Authoring; 2] { + [Authoring::Handwritten, Authoring::Macro] + } +} + +#[derive(Debug, Clone, Copy)] +enum CacheState { + Hot, + Cold, +} + +impl CacheState { + fn label(self) -> &'static str { + match self { + Self::Hot => "hot", + Self::Cold => "cold", + } + } + + fn all() -> [CacheState; 2] { + [CacheState::Hot, CacheState::Cold] + } +} + +fn id(workload: Workload, kind: SolverKind, authoring: Authoring, cache: CacheState) -> String { + format!( + "{}/{}/{}/{}", + workload.label(), + kind.label(), + authoring.label(), + cache.label() + ) +} + +fn predictions_group(c: &mut Criterion) { + let mut group = c.benchmark_group("native/predictions"); + group.sampling_mode(SamplingMode::Flat); + + for workload in Workload::all() { + let subject = subject_for_predictions(workload); + for kind in SolverKind::all() { + let theta = params(workload, kind); + for authoring in Authoring::all() { + for cache in CacheState::all() { + let bench_id = + BenchmarkId::from_parameter(id(workload, kind, authoring, cache)); + group.bench_function(bench_id, |b| match (kind, authoring, cache) { + (SolverKind::Ode, Authoring::Handwritten, CacheState::Hot) => { + let model = handwritten_ode(workload); + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + } + (SolverKind::Ode, Authoring::Handwritten, CacheState::Cold) => { + let model = handwritten_ode(workload).disable_cache(); + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + } + (SolverKind::Ode, Authoring::Macro, CacheState::Hot) => { + let model = macro_ode(workload); + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + } + (SolverKind::Ode, Authoring::Macro, CacheState::Cold) => { + let model = macro_ode(workload).disable_cache(); + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + } + (SolverKind::Analytical, Authoring::Handwritten, CacheState::Hot) => { + let model = handwritten_analytical(workload); + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + } + (SolverKind::Analytical, Authoring::Handwritten, CacheState::Cold) => { + let model = handwritten_analytical(workload).disable_cache(); + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + } + (SolverKind::Analytical, Authoring::Macro, CacheState::Hot) => { + let model = macro_analytical(workload); + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + } + (SolverKind::Analytical, Authoring::Macro, CacheState::Cold) => { + let model = macro_analytical(workload).disable_cache(); + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + } + (SolverKind::Sde, Authoring::Handwritten, CacheState::Hot) => { + let model = handwritten_sde(workload); + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + } + (SolverKind::Sde, Authoring::Handwritten, CacheState::Cold) => { + let model = handwritten_sde(workload).disable_cache(); + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + } + (SolverKind::Sde, Authoring::Macro, CacheState::Hot) => { + let model = macro_sde(workload); + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + } + (SolverKind::Sde, Authoring::Macro, CacheState::Cold) => { + let model = macro_sde(workload).disable_cache(); + b.iter(|| { + black_box( + model + .estimate_predictions( + black_box(&subject), + black_box(&theta), + ) + .unwrap(), + ) + }); + } + }); + } + } + } + } + + group.finish(); +} + +fn log_likelihood_group(c: &mut Criterion) { + let mut group = c.benchmark_group("native/log-likelihood"); + group.sampling_mode(SamplingMode::Flat); + let error_models = assay_error_models(); + + for workload in Workload::all() { + let subject = subject_for_likelihood(workload); + for kind in SolverKind::all() { + let theta = params(workload, kind); + for authoring in Authoring::all() { + for cache in CacheState::all() { + let bench_id = + BenchmarkId::from_parameter(id(workload, kind, authoring, cache)); + group.bench_function(bench_id, |b| match (kind, authoring, cache) { + (SolverKind::Ode, Authoring::Handwritten, CacheState::Hot) => { + let model = handwritten_ode(workload); + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + } + (SolverKind::Ode, Authoring::Handwritten, CacheState::Cold) => { + let model = handwritten_ode(workload).disable_cache(); + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + } + (SolverKind::Ode, Authoring::Macro, CacheState::Hot) => { + let model = macro_ode(workload); + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + } + (SolverKind::Ode, Authoring::Macro, CacheState::Cold) => { + let model = macro_ode(workload).disable_cache(); + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + } + (SolverKind::Analytical, Authoring::Handwritten, CacheState::Hot) => { + let model = handwritten_analytical(workload); + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + } + (SolverKind::Analytical, Authoring::Handwritten, CacheState::Cold) => { + let model = handwritten_analytical(workload).disable_cache(); + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + } + (SolverKind::Analytical, Authoring::Macro, CacheState::Hot) => { + let model = macro_analytical(workload); + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + } + (SolverKind::Analytical, Authoring::Macro, CacheState::Cold) => { + let model = macro_analytical(workload).disable_cache(); + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + } + (SolverKind::Sde, Authoring::Handwritten, CacheState::Hot) => { + let model = handwritten_sde(workload); + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + } + (SolverKind::Sde, Authoring::Handwritten, CacheState::Cold) => { + let model = handwritten_sde(workload).disable_cache(); + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + } + (SolverKind::Sde, Authoring::Macro, CacheState::Hot) => { + let model = macro_sde(workload); + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + } + (SolverKind::Sde, Authoring::Macro, CacheState::Cold) => { + let model = macro_sde(workload).disable_cache(); + b.iter(|| { + black_box( + model + .estimate_log_likelihood( + black_box(&subject), + black_box(&theta), + black_box(&error_models), + ) + .unwrap(), + ) + }); + } + }); + } + } + } + } + + group.finish(); +} + +fn likelihood_matrix_group(c: &mut Criterion) { + use pharmsol::prelude::simulator::log_likelihood_matrix; + + let mut group = c.benchmark_group("native/likelihood-matrix"); + group.sampling_mode(SamplingMode::Flat); + // One matrix iteration is multi-second; small sample budget keeps wall-clock sane. + group.sample_size(10); + group.measurement_time(Duration::from_secs(20)); + + let error_models = assay_error_models(); + + for workload in Workload::all() { + let data = matrix_data(workload, MATRIX_N_SUBJECTS); + for kind in SolverKind::all() { + let theta = support_points(workload, kind, MATRIX_N_SUPPORT); + for authoring in Authoring::all() { + let bench_id = BenchmarkId::from_parameter(format!( + "{}/{}/{}", + workload.label(), + kind.label(), + authoring.label() + )); + group.bench_function(bench_id, |b| match (kind, authoring) { + (SolverKind::Ode, Authoring::Handwritten) => { + let model = handwritten_ode(workload); + b.iter(|| { + black_box( + log_likelihood_matrix( + black_box(&model), + black_box(&data), + black_box(&theta), + black_box(&error_models), + false, + ) + .unwrap(), + ) + }); + } + (SolverKind::Ode, Authoring::Macro) => { + let model = macro_ode(workload); + b.iter(|| { + black_box( + log_likelihood_matrix( + black_box(&model), + black_box(&data), + black_box(&theta), + black_box(&error_models), + false, + ) + .unwrap(), + ) + }); + } + (SolverKind::Analytical, Authoring::Handwritten) => { + let model = handwritten_analytical(workload); + b.iter(|| { + black_box( + log_likelihood_matrix( + black_box(&model), + black_box(&data), + black_box(&theta), + black_box(&error_models), + false, + ) + .unwrap(), + ) + }); + } + (SolverKind::Analytical, Authoring::Macro) => { + let model = macro_analytical(workload); + b.iter(|| { + black_box( + log_likelihood_matrix( + black_box(&model), + black_box(&data), + black_box(&theta), + black_box(&error_models), + false, + ) + .unwrap(), + ) + }); + } + (SolverKind::Sde, Authoring::Handwritten) => { + let model = handwritten_sde(workload); + b.iter(|| { + black_box( + log_likelihood_matrix( + black_box(&model), + black_box(&data), + black_box(&theta), + black_box(&error_models), + false, + ) + .unwrap(), + ) + }); + } + (SolverKind::Sde, Authoring::Macro) => { + let model = macro_sde(workload); + b.iter(|| { + black_box( + log_likelihood_matrix( + black_box(&model), + black_box(&data), + black_box(&theta), + black_box(&error_models), + false, + ) + .unwrap(), + ) + }); + } + }); + } + } + } + + group.finish(); +} + +criterion_group!( + native_matrix, + predictions_group, + log_likelihood_group, + likelihood_matrix_group +); +criterion_main!(native_matrix); diff --git a/benches/nca.rs b/benches/nca.rs deleted file mode 100644 index d53413e1..00000000 --- a/benches/nca.rs +++ /dev/null @@ -1,100 +0,0 @@ -use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion}; -use pharmsol::nca::{NCAOptions, NCA}; -use pharmsol::prelude::*; -use std::hint::black_box; - -/// Build a typical PK subject with 12 time points (oral dose) -fn typical_oral_subject(id: &str) -> Subject { - Subject::builder(id) - .bolus(0.0, 100.0, 0) - .observation(0.0, 0.0, 0) - .observation(0.25, 2.5, 0) - .observation(0.5, 5.0, 0) - .observation(1.0, 8.0, 0) - .observation(2.0, 10.0, 0) - .observation(4.0, 7.5, 0) - .observation(6.0, 5.0, 0) - .observation(8.0, 3.5, 0) - .observation(12.0, 1.5, 0) - .observation(16.0, 0.8, 0) - .observation(24.0, 0.2, 0) - .observation(36.0, 0.05, 0) - .build() -} - -/// Build a population of n subjects with slight variation -fn build_population(n: usize) -> Data { - let subjects: Vec = (0..n) - .map(|i| { - let scale = 1.0 + (i as f64 % 7.0) * 0.05; // slight variation - Subject::builder(format!("subj_{}", i)) - .bolus(0.0, 100.0, 0) - .observation(0.0, 0.0, 0) - .observation(0.25, 2.5 * scale, 0) - .observation(0.5, 5.0 * scale, 0) - .observation(1.0, 8.0 * scale, 0) - .observation(2.0, 10.0 * scale, 0) - .observation(4.0, 7.5 * scale, 0) - .observation(6.0, 5.0 * scale, 0) - .observation(8.0, 3.5 * scale, 0) - .observation(12.0, 1.5 * scale, 0) - .observation(16.0, 0.8 * scale, 0) - .observation(24.0, 0.2 * scale, 0) - .observation(36.0, 0.05 * scale, 0) - .build() - }) - .collect(); - Data::new(subjects) -} - -fn bench_single_subject_nca(c: &mut Criterion) { - let subject = typical_oral_subject("bench_subj"); - let opts = NCAOptions::default(); - - c.bench_function("nca_single_subject", |b| { - b.iter(|| { - let result = black_box(&subject).nca(black_box(&opts)); - let _ = black_box(result); - }); - }); -} - -fn bench_population_nca(c: &mut Criterion) { - let mut group = c.benchmark_group("nca_population"); - - for size in [10, 100, 500] { - let data = build_population(size); - let opts = NCAOptions::default(); - - group.bench_with_input(BenchmarkId::from_parameter(size), &size, |b, _| { - b.iter(|| { - let results = black_box(&data).nca_all(black_box(&opts)); - black_box(results); - }); - }); - } - - group.finish(); -} - -fn bench_observation_metrics(c: &mut Criterion) { - use pharmsol::data::event::AUCMethod; - - let subject = typical_oral_subject("bench_subj"); - - c.bench_function("nca_auc_cmax_metrics", |b| { - b.iter(|| { - let auc = black_box(&subject).auc(0, &AUCMethod::Linear); - let cmax = black_box(&subject).cmax(0); - black_box((auc, cmax)); - }); - }); -} - -criterion_group!( - benches, - bench_single_subject_nca, - bench_population_nca, - bench_observation_metrics, -); -criterion_main!(benches); diff --git a/benches/ode.rs b/benches/ode.rs deleted file mode 100644 index c677fe94..00000000 --- a/benches/ode.rs +++ /dev/null @@ -1,123 +0,0 @@ -use criterion::{criterion_group, criterion_main, Criterion}; -use pharmsol::*; -use std::hint::black_box; - -fn example_subject() -> Subject { - Subject::builder("1") - .bolus(0.0, 100.0, 0) - .observation(3.0, 0.1, 0) - .observation(6.0, 0.4, 0) - .observation(12.0, 1.0, 0) - .observation(24.0, 1.1, 0) - .build() -} - -fn example_subject_covariates() -> Subject { - Subject::builder("1") - .bolus(0.0, 100.0, 0) - .covariate("iron", 0.0, 10.0) - .observation(3.0, 0.1, 0) - .covariate("iron", 3.0, 40.0) - .observation(6.0, 0.4, 0) - .covariate("iron", 6.0, 40.0) - .observation(12.0, 1.0, 0) - .covariate("iron", 12.0, 50.0) - .observation(24.0, 1.1, 0) - .covariate("iron", 24.0, 100.0) - .build() -} - -fn one_compartment() { - let subject = example_subject(); - let ode = equation::ODE::new( - |x, p, _t, dx, _b, _rateiv, _cov| { - fetch_params!(p, ka, ke, _tlag, _v); - dx[0] = -ka * x[0]; - dx[1] = ka * x[0] - ke * x[1]; - }, - |p, _t, _cov| { - fetch_params!(p, _ka, _ke, tlag, _v); - lag! {0=>tlag} - }, - |_p, _t, _cov| fa! {}, - |_p, _t, _cov, _x| {}, - |x, p, _t, _cov, y| { - fetch_params!(p, _ka, _ke, _tlag, v); - y[0] = x[1] / v; - }, - ) - .with_nstates(2) - .with_nout(1); - black_box( - ode.estimate_predictions(&subject, &[0.3, 0.5, 0.1, 70.0]) - .unwrap(), - ); -} - -fn one_compartment_covariates() { - let subject = example_subject_covariates(); - let ode = equation::ODE::new( - |x, p, t, dx, _b, _rateiv, cov| { - fetch_cov!(cov, t, iron); - fetch_params!(p, ka, ke, _tlag, _v); - - ke = ke * iron.powf(0.75) / 50.0; - dx[0] = -ka * x[0]; - dx[1] = ka * x[0] - ke * x[1]; - }, - |p, _t, _cov| { - fetch_params!(p, _ka, _ke, tlag, _v); - lag! {0=>tlag} - }, - |_p, _t, _cov| fa! {}, - |_p, _t, _cov, _x| {}, - |x, p, _t, _cov, y| { - fetch_params!(p, _ka, _ke, _tlag, v); - y[0] = x[1] / v; - }, - ) - .with_nstates(2) - .with_nout(1); - black_box( - ode.estimate_predictions(&subject, &[0.3, 0.5, 0.1, 70.0]) - .unwrap(), - ); -} - -fn two_compartment() { - let subject = example_subject(); - let ode = equation::ODE::new( - |x, p, _t, dx, _b, _rateiv, _cov| { - fetch_params!(p, ka, ke, k12, k21, _tlag, _v); - dx[0] = -ka * x[0] - k12 * x[0] + k21 * x[1]; - dx[1] = k12 * x[0] - k21 * x[1] - ke * x[1]; - }, - |p, _t, _cov| { - fetch_params!(p, _ka, _ke, _k12, _k21, tlag, _v); - lag! {0=>tlag} - }, - |_p, _t, _cov| fa! {}, - |_p, _t, _cov, _x| {}, - |x, p, _t, _cov, y| { - fetch_params!(p, _ka, _ke, _k12, _k21, _tlag, v); - y[0] = x[1] / v; - }, - ) - .with_nstates(2) - .with_nout(1); - black_box( - ode.estimate_predictions(&subject, &[0.3, 0.5, 0.1, 0.04, 0.08, 70.0]) - .unwrap(), - ); -} - -fn criterion_benchmark(c: &mut Criterion) { - c.bench_function("one_compartment", |b| b.iter(one_compartment)); - c.bench_function("one_compartment_covariates", |b| { - b.iter(one_compartment_covariates) - }); - c.bench_function("two_compartment", |b| b.iter(two_compartment)); -} - -criterion_group!(benches, criterion_benchmark); -criterion_main!(benches); diff --git a/benches/performance.rs b/benches/performance.rs deleted file mode 100644 index 2ab17403..00000000 --- a/benches/performance.rs +++ /dev/null @@ -1,55 +0,0 @@ -use criterion::{criterion_group, criterion_main, Criterion}; -use pharmsol::*; -use std::hint::black_box; - -fn example_subject() -> Subject { - Subject::builder("id1") - .bolus(0.0, 100.0, 0) - .repeat(2, 0.5) - .observation(0.5, 0.1, 0) - .observation(1.0, 0.4, 0) - .observation(2.0, 1.0, 0) - .observation(2.5, 1.1, 0) - .covariate("wt", 0.0, 80.0) - .covariate("wt", 1.0, 83.0) - .covariate("age", 0.0, 25.0) - .build() -} - -fn readme(n: usize) { - let subject = example_subject(); - let ode = equation::ODE::new( - |x, p, _t, dx, _b, _rateiv, _cov| { - // fetch_cov!(cov, t, _wt, _age); - fetch_params!(p, ka, ke, _tlag, _v); - //Struct - dx[0] = -ka * x[0]; - dx[1] = ka * x[0] - ke * x[1]; - }, - |p, _t, _cov| { - fetch_params!(p, _ka, _ke, tlag, _v); - lag! {0=>tlag} - }, - |_p, _t, _cov| fa! {}, - |_p, _t, _cov, _x| {}, - |x, p, _t, _cov, y| { - fetch_params!(p, _ka, _ke, _tlag, v); - y[0] = x[1] / v; - }, - ) - .with_nstates(2) - .with_nout(1); - for _ in 0..n { - let op = ode - .estimate_predictions(&subject, &[0.3, 0.5, 0.1, 70.0]) - .unwrap(); - black_box(op); - } -} - -fn criterion_benchmark(c: &mut Criterion) { - c.bench_function("readme 20", |b| b.iter(|| readme(black_box(20)))); -} - -criterion_group!(benches, criterion_benchmark); -criterion_main!(benches); diff --git a/benches/runtime_matrix.rs b/benches/runtime_matrix.rs deleted file mode 100644 index 81648b8d..00000000 --- a/benches/runtime_matrix.rs +++ /dev/null @@ -1,100 +0,0 @@ -#[path = "../tests/support/runtime_corpus.rs"] -mod runtime_corpus; - -use criterion::{criterion_group, criterion_main, Criterion}; - -#[cfg(all( - feature = "dsl-jit", - feature = "dsl-aot", - feature = "dsl-aot-load", - feature = "dsl-wasm" -))] -fn runtime_matrix_benchmark(c: &mut Criterion) { - use std::hint::black_box; - use std::time::Duration; - - use criterion::{BatchSize, BenchmarkId, SamplingMode}; - use pharmsol::dsl::WasmCompileCache; - use runtime_corpus as corpus; - use runtime_corpus::{ArtifactWorkspace, CorpusCase}; - - let cases = [CorpusCase::Ode, CorpusCase::Analytical, CorpusCase::Sde]; - - let mut compile_group = c.benchmark_group("dsl/runtime-compile"); - compile_group.sampling_mode(SamplingMode::Flat); - compile_group.sample_size(10); - compile_group.warm_up_time(Duration::from_millis(500)); - compile_group.measurement_time(Duration::from_secs(1)); - - for case in cases { - compile_group.bench_with_input(BenchmarkId::new("jit", case.label()), &case, |b, case| { - b.iter(|| black_box(corpus::compile_runtime_jit_model(*case).unwrap())) - }); - compile_group.bench_with_input( - BenchmarkId::new("native-aot", case.label()), - &case, - |b, case| { - b.iter_batched( - || ArtifactWorkspace::new().unwrap(), - |workspace| { - black_box( - corpus::compile_runtime_native_aot_model(*case, &workspace).unwrap(), - ) - }, - BatchSize::PerIteration, - ) - }, - ); - compile_group.bench_with_input(BenchmarkId::new("wasm", case.label()), &case, |b, case| { - b.iter(|| black_box(corpus::compile_runtime_wasm_model(*case).unwrap())) - }); - compile_group.bench_with_input( - BenchmarkId::new("wasm-module", case.label()), - &case, - |b, case| b.iter(|| black_box(corpus::compile_wasm_module(*case).unwrap())), - ); - - let cache = WasmCompileCache::default(); - corpus::compile_wasm_module_with_cache(case, &cache).unwrap(); - compile_group.bench_function(BenchmarkId::new("wasm-module-cached", case.label()), |b| { - b.iter(|| black_box(corpus::compile_wasm_module_with_cache(case, &cache).unwrap())) - }); - } - compile_group.finish(); - - let mut runtime_group = c.benchmark_group("dsl/runtime-predict"); - runtime_group.sampling_mode(SamplingMode::Flat); - runtime_group.sample_size(10); - runtime_group.warm_up_time(Duration::from_millis(500)); - runtime_group.measurement_time(Duration::from_secs(1)); - - for case in cases { - let jit = corpus::compile_runtime_jit_model(case).unwrap(); - runtime_group.bench_function(BenchmarkId::new("jit", case.label()), |b| { - b.iter(|| black_box(corpus::estimate_runtime_predictions(case, &jit).unwrap())) - }); - - let aot_workspace = ArtifactWorkspace::new().unwrap(); - let aot = corpus::compile_runtime_native_aot_model(case, &aot_workspace).unwrap(); - runtime_group.bench_function(BenchmarkId::new("native-aot", case.label()), |b| { - b.iter(|| black_box(corpus::estimate_runtime_predictions(case, &aot).unwrap())) - }); - - let wasm = corpus::compile_runtime_wasm_model(case).unwrap(); - runtime_group.bench_function(BenchmarkId::new("wasm", case.label()), |b| { - b.iter(|| black_box(corpus::estimate_runtime_predictions(case, &wasm).unwrap())) - }); - } - runtime_group.finish(); -} - -#[cfg(not(all( - feature = "dsl-jit", - feature = "dsl-aot", - feature = "dsl-aot-load", - feature = "dsl-wasm" -)))] -fn runtime_matrix_benchmark(_: &mut Criterion) {} - -criterion_group!(benches, runtime_matrix_benchmark); -criterion_main!(benches); diff --git a/src/dsl/native.rs b/src/dsl/native.rs index 7197084e..7c2e327c 100644 --- a/src/dsl/native.rs +++ b/src/dsl/native.rs @@ -24,11 +24,11 @@ use crate::{ data::error_model::AssayErrorModels, data::{Covariates, Infusion, InputLabel, OutputLabel}, simulator::{ - cache::{PredictionCache, DEFAULT_CACHE_SIZE}, + cache::{PredictionCache, SdeLikelihoodCache, DEFAULT_CACHE_SIZE}, equation::{ ode::{closure_helpers::PMProblem, ExplicitRkTableau, OdeSolver, SdirkTableau}, sde::simulate_sde_event_with, - EqnKind, Equation, EquationPriv, EquationTypes, + EqnKind, Equation, EquationPriv, EquationTypes, Predictions, }, likelihood::{Prediction, SubjectPredictions}, Fa, Lag, M, T, V, @@ -745,12 +745,14 @@ pub struct NativeOdeModel { #[derive(Clone, Debug)] pub struct NativeAnalyticalModel { shared: Arc, + cache: Option, } #[derive(Clone, Debug)] pub struct NativeSdeModel { shared: Arc, nparticles: usize, + cache: Option, } #[derive(Clone, Debug)] @@ -1224,6 +1226,7 @@ impl NativeAnalyticalModel { pub(crate) fn new(info: NativeModelInfo, artifact: impl RuntimeArtifact + 'static) -> Self { Self { shared: Arc::new(SharedNativeModel::new(info, artifact)), + cache: Some(PredictionCache::new(DEFAULT_CACHE_SIZE)), } } @@ -1377,12 +1380,189 @@ impl NativeAnalyticalModel { } } +#[inline(always)] +fn runtime_analytical_predictions( + model: &NativeAnalyticalModel, + subject: &Subject, + support_point: &[f64], +) -> Result { + if let Some(cache) = &model.cache { + let key = ( + subject.hash(), + crate::simulator::equation::spphash(support_point), + ); + if let Some(cached) = cache.get(&key) { + return Ok(cached); + } + + let result = model.estimate_predictions(subject, support_point)?; + cache.insert(key, result.clone()); + Ok(result) + } else { + model.estimate_predictions(subject, support_point) + } +} + +impl crate::simulator::equation::Cache for NativeAnalyticalModel { + fn with_cache_capacity(mut self, size: u64) -> Self { + self.cache = Some(PredictionCache::new(size)); + self + } + + fn enable_cache(mut self) -> Self { + self.cache = Some(PredictionCache::new(DEFAULT_CACHE_SIZE)); + self + } + + fn clear_cache(&self) { + if let Some(cache) = &self.cache { + cache.invalidate_all(); + } + } + + fn disable_cache(mut self) -> Self { + self.cache = None; + self + } +} + +impl EquationTypes for NativeAnalyticalModel { + type S = V; + type P = SubjectPredictions; +} + +impl EquationPriv for NativeAnalyticalModel { + fn lag(&self) -> &Lag { + &(runtime_no_lag as Lag) + } + + fn fa(&self) -> &Fa { + &(runtime_no_fa as Fa) + } + + fn get_nstates(&self) -> usize { + self.shared.info.state_len + } + + fn get_ndrugs(&self) -> usize { + self.shared.info.route_len + } + + fn get_nouteqs(&self) -> usize { + self.shared.info.output_len + } + + fn metadata(&self) -> Option<&crate::ValidatedModelMetadata> { + None + } + + fn solve( + &self, + _state: &mut Self::S, + _support_point: &[f64], + _covariates: &Covariates, + _infusions: &[Infusion], + _start_time: f64, + _end_time: f64, + ) -> Result<(), PharmsolError> { + unimplemented!("solve is not used for runtime analytical models") + } + + fn process_observation( + &self, + _support_point: &[f64], + _observation: &Observation, + _error_models: Option<&AssayErrorModels>, + _time: f64, + _covariates: &Covariates, + _x: &mut Self::S, + _likelihood: &mut Vec, + _output: &mut Self::P, + ) -> Result<(), PharmsolError> { + unimplemented!("process_observation is not used for runtime analytical models") + } + + fn initial_state( + &self, + _support_point: &[f64], + _covariates: &Covariates, + _occasion_index: usize, + ) -> Self::S { + V::zeros(self.shared.info.state_len, NalgebraContext) + } +} + +impl Equation for NativeAnalyticalModel { + fn estimate_likelihood( + &self, + subject: &Subject, + support_point: &[f64], + error_models: &AssayErrorModels, + ) -> Result { + Ok(self + .estimate_log_likelihood(subject, support_point, error_models)? + .exp()) + } + + fn estimate_log_likelihood( + &self, + subject: &Subject, + support_point: &[f64], + error_models: &AssayErrorModels, + ) -> Result { + let bound_error_models = self.bind_error_models(error_models)?; + let predictions = runtime_analytical_predictions(self, subject, support_point)?; + predictions.log_likelihood(&bound_error_models) + } + + fn kind() -> EqnKind { + EqnKind::Analytical + } + + fn assay_error_models(&self) -> AssayErrorModels { + AssayErrorModels::with_output_names( + self.info() + .outputs + .iter() + .map(|output| output.name.as_str()), + ) + } + + fn estimate_predictions( + &self, + subject: &Subject, + support_point: &[f64], + ) -> Result { + runtime_analytical_predictions(self, subject, support_point) + } + + fn simulate_subject( + &self, + subject: &Subject, + support_point: &[f64], + error_models: Option<&AssayErrorModels>, + ) -> Result<(Self::P, Option), PharmsolError> { + let bound_error_models = match error_models { + Some(error_models) => Some(self.bind_error_models(error_models)?), + None => None, + }; + + let predictions = runtime_analytical_predictions(self, subject, support_point)?; + let likelihood = match bound_error_models.as_ref() { + Some(error_models) => Some(predictions.log_likelihood(error_models)?.exp()), + None => None, + }; + Ok((predictions, likelihood)) + } +} + impl NativeSdeModel { pub(crate) fn new(info: NativeModelInfo, artifact: impl RuntimeArtifact + 'static) -> Self { let nparticles = info.particles.unwrap_or(1); Self { shared: Arc::new(SharedNativeModel::new(info, artifact)), nparticles, + cache: Some(SdeLikelihoodCache::new(DEFAULT_CACHE_SIZE)), } } @@ -1633,6 +1813,192 @@ impl NativeSdeModel { } } +#[inline(always)] +fn runtime_sde_log_likelihood( + model: &NativeSdeModel, + subject: &Subject, + support_point: &[f64], + error_models: &AssayErrorModels, +) -> Result { + if let Some(cache) = &model.cache { + let key = ( + subject.hash(), + crate::simulator::equation::spphash(support_point), + error_models.hash(), + ); + if let Some(cached) = cache.get(&key) { + return Ok(cached); + } + + let predictions = model.estimate_predictions(subject, support_point)?; + let log_lik = predictions.log_likelihood(error_models)?; + cache.insert(key, log_lik); + Ok(log_lik) + } else { + let predictions = model.estimate_predictions(subject, support_point)?; + predictions.log_likelihood(error_models) + } +} + +impl crate::simulator::equation::Cache for NativeSdeModel { + fn with_cache_capacity(mut self, size: u64) -> Self { + self.cache = Some(SdeLikelihoodCache::new(size)); + self + } + + fn enable_cache(mut self) -> Self { + self.cache = Some(SdeLikelihoodCache::new(DEFAULT_CACHE_SIZE)); + self + } + + fn clear_cache(&self) { + if let Some(cache) = &self.cache { + cache.invalidate_all(); + } + } + + fn disable_cache(mut self) -> Self { + self.cache = None; + self + } +} + +impl EquationTypes for NativeSdeModel { + type S = Vec>; + type P = Array2; +} + +impl EquationPriv for NativeSdeModel { + fn lag(&self) -> &Lag { + &(runtime_no_lag as Lag) + } + + fn fa(&self) -> &Fa { + &(runtime_no_fa as Fa) + } + + fn get_nstates(&self) -> usize { + self.shared.info.state_len + } + + fn get_ndrugs(&self) -> usize { + self.shared.info.route_len + } + + fn get_nouteqs(&self) -> usize { + self.shared.info.output_len + } + + fn nparticles(&self) -> usize { + self.nparticles + } + + fn is_sde(&self) -> bool { + true + } + + fn metadata(&self) -> Option<&crate::ValidatedModelMetadata> { + None + } + + fn solve( + &self, + _state: &mut Self::S, + _support_point: &[f64], + _covariates: &Covariates, + _infusions: &[Infusion], + _start_time: f64, + _end_time: f64, + ) -> Result<(), PharmsolError> { + unimplemented!("solve is not used for runtime SDE models") + } + + fn process_observation( + &self, + _support_point: &[f64], + _observation: &Observation, + _error_models: Option<&AssayErrorModels>, + _time: f64, + _covariates: &Covariates, + _x: &mut Self::S, + _likelihood: &mut Vec, + _output: &mut Self::P, + ) -> Result<(), PharmsolError> { + unimplemented!("process_observation is not used for runtime SDE models") + } + + fn initial_state( + &self, + _support_point: &[f64], + _covariates: &Covariates, + _occasion_index: usize, + ) -> Self::S { + vec![DVector::zeros(self.shared.info.state_len); self.nparticles] + } +} + +impl Equation for NativeSdeModel { + fn estimate_likelihood( + &self, + subject: &Subject, + support_point: &[f64], + error_models: &AssayErrorModels, + ) -> Result { + let log_lik = self.estimate_log_likelihood(subject, support_point, error_models)?; + Ok(log_lik.exp()) + } + + fn estimate_log_likelihood( + &self, + subject: &Subject, + support_point: &[f64], + error_models: &AssayErrorModels, + ) -> Result { + let bound_error_models = self.bind_error_models(error_models)?; + runtime_sde_log_likelihood(self, subject, support_point, &bound_error_models) + } + + fn kind() -> EqnKind { + EqnKind::SDE + } + + fn assay_error_models(&self) -> AssayErrorModels { + AssayErrorModels::with_output_names( + self.info() + .outputs + .iter() + .map(|output| output.name.as_str()), + ) + } + + fn estimate_predictions( + &self, + subject: &Subject, + support_point: &[f64], + ) -> Result { + NativeSdeModel::estimate_predictions(self, subject, support_point) + } + + fn simulate_subject( + &self, + subject: &Subject, + support_point: &[f64], + error_models: Option<&AssayErrorModels>, + ) -> Result<(Self::P, Option), PharmsolError> { + let bound_error_models = match error_models { + Some(error_models) => Some(self.bind_error_models(error_models)?), + None => None, + }; + + let predictions = NativeSdeModel::estimate_predictions(self, subject, support_point)?; + let likelihood = match bound_error_models.as_ref() { + Some(error_models) => Some(predictions.log_likelihood(error_models)?.exp()), + None => None, + }; + Ok((predictions, likelihood)) + } +} + fn active_route_inputs(infusions: &[Infusion], time: f64, route_len: usize) -> Vec { let mut values = vec![0.0; route_len]; for infusion in infusions {