refactor: Remove likelihood calculation

.github/workflows/pr_benchmarks.yml

@@ -10,6 +10,7 @@ jobs:
     # DO NOT REMOVE: For handling Fork PRs see Pull Requests from Forks
     if: github.event_name == 'pull_request' && github.event.pull_request.head.repo.full_name == github.repository
     permissions:
+      checks: write
       pull-requests: write
     runs-on: self-hosted
     steps:

CHANGELOG.md

@@ -38,12 +38,26 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
-- Support time in the DSL ([#310](https://github.com/LAPKB/pharmsol/pull/310))
-- Rename DSL types and add additional safeguards ([#302](https://github.com/LAPKB/pharmsol/pull/302))
+- Add a caller-controlled, simulation-neutral SDE particle session that pauses
+  at observation boundaries and can resume with retained or replaced states.
+- Add focused Criterion coverage for standard SDE prediction and particle-session
+  retain/select paths.
 
-### Other
+### Removed
+
+- Remove equation-level and prediction-container observation evaluation APIs;
+  downstream fitting crates now generate predictions with pharmsol and evaluate
+  observations themselves.
+- Remove residual-distribution declarations, estimator caches, and parameter
+  optimization helpers from pharmsol's public API.
+
+### Changed
 
-- Make CI faster ([#303](https://github.com/LAPKB/pharmsol/pull/303))
+- Relocate the existing data-only `ErrorPoly` DTO while preserving its public
+  path and unchanged C0-C3 transport through observations and predictions.
+  As before, incomplete C0-C3 rows do not create an `ErrorPoly` value.
+- Keep analytical, ODE, and SDE execution focused exclusively on simulation and
+  prediction generation, including simulation-only examples and benchmarks.
 
 ## [0.28.2](https://github.com/LAPKB/pharmsol/compare/pharmsol-v0.28.1...pharmsol-v0.28.2) - 2026-07-20
 
@@ -133,7 +147,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 
-- Handle missing observations in log-likelihood ([#230](https://github.com/LAPKB/pharmsol/pull/230))
+- Handle missing observations during observation evaluation ([#230](https://github.com/LAPKB/pharmsol/pull/230))
 
 ## [0.24.0](https://github.com/LAPKB/pharmsol/compare/v0.23.0...v0.24.0) - 2026-03-23
 
@@ -146,7 +160,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
-- Remove LikelihoodMatrixOptions ([#222](https://github.com/LAPKB/pharmsol/pull/222))
+- Remove obsolete matrix configuration options ([#222](https://github.com/LAPKB/pharmsol/pull/222))
 - Use a sized cache with new subject hashing ([#214](https://github.com/LAPKB/pharmsol/pull/214))
 - Normalized rows ([#205](https://github.com/LAPKB/pharmsol/pull/205))
 
@@ -246,7 +260,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Exhaustive match instead of if-else block
 - Add setters
 - Make cens column optional
-- Censoring likelihood calculation
+- Add censored-observation evaluation
 - Update argmin requirement from 0.10.0 to 0.11.0 ([#135](https://github.com/LAPKB/pharmsol/pull/135))
 
 ## [0.17.2](https://github.com/LAPKB/pharmsol/compare/v0.17.1...v0.17.2) - 2025-09-30

Cargo.toml

@@ -31,11 +31,11 @@ documentation = "https://lapkb.github.io/pharmsol/"
 default = []
 dsl-core = []
 dsl-jit = [
-    "dsl-core",
-    "dep:cranelift",
-    "dep:cranelift-jit",
-    "dep:cranelift-module",
-    "dep:cranelift-native",
+  "dsl-core",
+  "dep:cranelift",
+  "dep:cranelift-jit",
+  "dep:cranelift-module",
+  "dep:cranelift-native",
 ]
 dsl-aot = ["dsl-core"]
 dsl-aot-load = ["dsl-core", "dep:libloading"]
@@ -87,6 +87,10 @@ bench = false
 name = "native_matrix"
 harness = false
 
+[[bench]]
+name = "sde_sessions"
+harness = false
+
 [[bench]]
 name = "dsl_matrix"
 harness = false

README.md

@@ -91,9 +91,7 @@ See [examples/analytical_readme.rs](examples/analytical_readme.rs),
 [examples/sde_readme.rs](examples/sde_readme.rs),
 [examples/dsl_jit_analytical_covariates.rs](examples/dsl_jit_analytical_covariates.rs),
 [examples/analytical_vs_ode.rs](examples/analytical_vs_ode.rs), and
-[examples/compare_solvers.rs](examples/compare_solvers.rs). For migration-oriented notes,
-see [docs/analytical-authoring-migration.md](docs/analytical-authoring-migration.md) and
-[docs/ode-authoring-migration.md](docs/ode-authoring-migration.md).
+[examples/compare_solvers.rs](examples/compare_solvers.rs).
 
 ### Built-In Analytical Kernels
 

benches/common/mod.rs

@@ -1,4 +1,4 @@
-//! Shared bench fixtures (workloads, subjects, params, error models, model factories)
+//! Shared bench fixtures (workloads, subjects, parameters, and 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.
 //!
@@ -11,15 +11,11 @@
 
 #![allow(dead_code)]
 
-use ndarray::Array2;
+use pharmsol::equation::{self, Analytical, Route, ODE, SDE};
 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 {
@@ -106,7 +102,7 @@ impl SolverKind {
         }
     }
 
-    // SDE bench cells temporarily disabled — too slow for the current matrix.
+    // Stochastic workloads use the dedicated SDE benchmark target.
     pub fn all() -> [SolverKind; 2] {
         [SolverKind::Ode, SolverKind::Analytical]
     }
@@ -117,18 +113,11 @@ impl SolverKind {
 /// 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());
@@ -153,60 +142,6 @@ pub fn subject_for_predictions(workload: Workload) -> Subject {
     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 named parameters per `(workload, kind)` in the source order shared by bench fixtures.
@@ -231,45 +166,6 @@ pub fn named_params(workload: Workload, kind: SolverKind) -> Vec<(&'static str,
     }
 }
 
-/// 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<f64> {
-    named_params(workload, kind)
-        .into_iter()
-        .map(|(_, value)| value)
-        .collect()
-}
-
-/// 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<f64> {
-    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 {