feat: Non-compartmental analysis

Cargo.toml

@@ -52,3 +52,7 @@ harness = false
 [[bench]]
 name = "analytical_vs_ode"
 harness = false
+
+[[bench]]
+name = "nca"
+harness = false

README.md

@@ -87,7 +87,42 @@ let ode = equation::ODE::new(
 
 Analytical solutions provide 20-33× speedups compared to equivalent ODE formulations. See [benchmarks](benches/) for details.
 
+## Non-Compartmental Analysis (NCA)
 
+pharmsol includes a complete NCA module for calculating standard pharmacokinetic parameters.
+
+```rust
+use pharmsol::prelude::*;
+use pharmsol::nca::NCAOptions;
+
+let subject = Subject::builder("patient_001")
+    .bolus(0.0, 100.0, 0)  // 100 mg oral dose
+    .observation(0.5, 5.0, 0)
+    .observation(1.0, 10.0, 0)
+    .observation(2.0, 8.0, 0)
+    .observation(4.0, 4.0, 0)
+    .observation(8.0, 2.0, 0)
+    .build();
+
+let results = subject.nca(&NCAOptions::default(), 0);
+let result = results[0].as_ref().expect("NCA failed");
+
+println!("Cmax: {:.2}", result.exposure.cmax);
+println!("Tmax: {:.2} h", result.exposure.tmax);
+println!("AUClast: {:.2}", result.exposure.auc_last);
+
+if let Some(ref term) = result.terminal {
+    println!("Half-life: {:.2} h", term.half_life);
+}
+```
+
+**Supported NCA Parameters:**
+
+- Exposure: Cmax, Tmax, Clast, Tlast, AUClast, AUCinf, tlag
+- Terminal: λz, t½, MRT
+- Clearance: CL/F, Vz/F, Vss
+- IV-specific: C0 (back-extrapolation), Vd
+- Steady-state: AUCtau, Cmin, Cavg, fluctuation, swing
 
 # Links
 

benches/nca.rs

@@ -0,0 +1,127 @@
+use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion};
+use pharmsol::nca::{lambda_z_candidates, 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<Subject> = (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_lambda_z_candidates(c: &mut Criterion) {
+    use pharmsol::data::event::{AUCMethod, BLQRule};
+    use pharmsol::data::observation::ObservationProfile;
+    use pharmsol::nca::LambdaZOptions;
+
+    let subject = typical_oral_subject("bench_subj");
+    let occ = &subject.occasions()[0];
+    let profile = ObservationProfile::from_occasion(occ, 0, &BLQRule::Exclude).unwrap();
+    let lz_opts = LambdaZOptions::default();
+
+    // Get AUClast for the candidate scoring
+    let auc_results = subject.auc(0, &AUCMethod::Linear, &BLQRule::Exclude);
+    let auc_last = auc_results[0].as_ref().copied().unwrap_or(50.0);
+
+    c.bench_function("nca_lambda_z_candidates", |b| {
+        b.iter(|| {
+            let candidates = lambda_z_candidates(
+                black_box(&profile),
+                black_box(&lz_opts),
+                black_box(auc_last),
+            );
+            black_box(candidates);
+        });
+    });
+}
+
+fn bench_observation_metrics(c: &mut Criterion) {
+    use pharmsol::data::event::{AUCMethod, BLQRule};
+
+    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, &BLQRule::Exclude);
+            let cmax = black_box(&subject).cmax(0, &BLQRule::Exclude);
+            black_box((auc, cmax));
+        });
+    });
+}
+
+criterion_group!(
+    benches,
+    bench_single_subject_nca,
+    bench_population_nca,
+    bench_lambda_z_candidates,
+    bench_observation_metrics,
+);
+criterion_main!(benches);

examples/exa.rs

@@ -13,14 +13,20 @@ fn main() {
     use std::path::PathBuf;
 
     // Create test subject with infusion and observations
+    // Including missing observations to verify predictions work without observed values
     let subject = Subject::builder("1")
         .infusion(0.0, 500.0, 0, 0.5)
         .observation(0.5, 1.645776, 0)
+        .missing_observation(0.75, 0) // Missing observation
         .observation(1.0, 1.216442, 0)
+        .missing_observation(1.5, 0) // Missing observation
         .observation(2.0, 0.4622729, 0)
+        .missing_observation(2.5, 0) // Missing observation
         .observation(3.0, 0.1697458, 0)
         .observation(4.0, 0.06382178, 0)
+        .missing_observation(5.0, 0) // Missing observation
         .observation(6.0, 0.009099384, 0)
+        .missing_observation(7.0, 0) // Missing observation
         .observation(8.0, 0.001017932, 0)
         .build();
 
@@ -138,22 +144,32 @@ fn main() {
     let dynamic_ode_flat = dynamic_ode_preds.flat_predictions();
     let dynamic_analytical_flat = dynamic_analytical_preds.flat_predictions();
 
+    let static_times = static_ode_preds.flat_times();
+    let static_obs = static_ode_preds.flat_observations();
+
     println!(
-        "\n{:<12} {:>15} {:>15} {:>15}",
-        "Time", "Static ODE", "Dynamic ODE", "Analytical"
+        "\n{:<12} {:>12} {:>15} {:>15} {:>15}",
+        "Time", "Obs", "Static ODE", "Dynamic ODE", "Analytical"
     );
-    println!("{}", "-".repeat(60));
+    println!("{}", "-".repeat(75));
 
-    let times = [0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0];
-    for (i, &time) in times.iter().enumerate() {
+    for i in 0..static_times.len() {
+        let obs_str = match static_obs[i] {
+            Some(v) => format!("{:.4}", v),
+            None => "MISSING".to_string(),
+        };
         println!(
-            "{:<12.1} {:>15.6} {:>15.6} {:>15.6}",
-            time, static_flat[i], dynamic_ode_flat[i], dynamic_analytical_flat[i]
+            "{:<12.2} {:>12} {:>15.6} {:>15.6} {:>15.6}",
+            static_times[i],
+            obs_str,
+            static_flat[i],
+            dynamic_ode_flat[i],
+            dynamic_analytical_flat[i]
         );
     }
 
     // Verify predictions match
-    println!("\n{}", "=".repeat(60));
+    println!("\n{}", "=".repeat(75));
     println!("Verification:");
 
     let ode_match = static_flat
@@ -182,6 +198,10 @@ fn main() {
         }
     );
 
+    // Count zero predictions for missing observations
+    let zero_count = static_flat.iter().filter(|&&v| v == 0.0).count();
+    println!("  Zero predictions count: {} (should be 0)", zero_count);
+
     // =========================================================================
     // 5. Clean up compiled model files
     // =========================================================================