Skip to content

feat: Remove API for regular likelihood #188

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mhovd wants to merge 1 commit into
base: main
Choose a base branch
from
Open

feat: Remove API for regular likelihood #188

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 15 additions & 15 deletions examples/one_compartment.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,28 +67,28 @@ fn main() -> Result<(), pharmsol::PharmsolError> {
let ke = 1.022; // Elimination rate constant
let v = 194.0; // Volume of distribution

// Compute likelihoods and predictions for both models
let analytical_likelihoods = an.estimate_likelihood(&subject, &vec![ke, v], &ems, false)?;
// Compute log-likelihoods and predictions for both models
let analytical_log_lik = an.estimate_log_likelihood(&subject, &vec![ke, v], &ems, false)?;

let analytical_predictions = an.estimate_predictions(&subject, &vec![ke, v])?;

let ode_likelihoods = ode.estimate_likelihood(&subject, &vec![ke, v], &ems, false)?;
let ode_log_lik = ode.estimate_log_likelihood(&subject, &vec![ke, v], &ems, false)?;

let ode_predictions = ode.estimate_predictions(&subject, &vec![ke, v])?;

// Print comparison table
println!("\n┌───────────┬─────────────────┬─────────────────┬─────────────────────┐");
println!("│ │ Analytical │ ODE │ Difference │");
println!("├───────────┼─────────────────┼─────────────────┼─────────────────────┤");
println!("\n┌─────────────────┬─────────────────┬─────────────────┬─────────────────────┐");
println!("│ │ Analytical │ ODE │ Difference │");
println!("├─────────────────┼─────────────────┼─────────────────┼─────────────────────┤");
println!(
"│ Likelihood│ {:>15.6} │ {:>15.6} │ {:>19.2e} │",
analytical_likelihoods,
ode_likelihoods,
analytical_likelihoods - ode_likelihoods
"│ Log-Likelihood │ {:>15.6} │ {:>15.6} │ {:>19.2e} │",
analytical_log_lik,
ode_log_lik,
analytical_log_lik - ode_log_lik
);
println!("├───────────┼─────────────────┼─────────────────┼─────────────────────┤");
println!("│ Time │ Prediction │ Prediction │ │");
println!("├───────────┼─────────────────┼─────────────────┼─────────────────────┤");
println!("├─────────────────┼─────────────────┼─────────────────┼─────────────────────┤");
println!("│ Time │ Prediction │ Prediction │ │");
println!("├─────────────────┼─────────────────┼─────────────────┼─────────────────────┤");

let times = analytical_predictions.flat_times();
let analytical_preds = analytical_predictions.flat_predictions();
Expand All @@ -101,12 +101,12 @@ fn main() -> Result<(), pharmsol::PharmsolError> {
{
let diff = a - b;
println!(
"│ {:>9.2} │ {:>15.9} │ {:>15.9} │ {:>19.2e} │",
"│ {:>15.2} │ {:>15.9} │ {:>15.9} │ {:>19.2e} │",
t, a, b, diff
);
}

println!("└───────────┴─────────────────┴─────────────────┴─────────────────────┘\n");
println!("└─────────────────┴─────────────────┴─────────────────┴─────────────────────┘\n");

Ok(())
}
4 changes: 2 additions & 2 deletions examples/pf.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,11 +37,11 @@ fn main() {
.unwrap();

let ll = sde
.estimate_likelihood(&subject, &vec![1.0], &ems, false)
.estimate_log_likelihood(&subject, &vec![1.0], &ems, false)
.unwrap();

dbg!(sde
.estimate_likelihood(&subject, &vec![1.0], &ems, false)
.estimate_log_likelihood(&subject, &vec![1.0], &ems, false)
.unwrap());
println!("{ll:#?}");
}
2 changes: 1 addition & 1 deletion src/lib.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ pub mod prelude {
pub use crate::simulator::{
equation,
equation::Equation,
likelihood::{log_psi, psi, PopulationPredictions, Prediction, SubjectPredictions},
likelihood::{psi, PopulationPredictions, Prediction, SubjectPredictions},
};
}
pub mod models {
Expand Down
27 changes: 14 additions & 13 deletions src/optimize/spp.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ pub struct SppOptimizer<'a, E: Equation> {
equation: &'a E,
data: &'a Data,
sig: &'a ErrorModels,
pyl: &'a Array1<f64>,
log_pyl: &'a Array1<f64>,
}

impl<E: Equation> CostFunction for SppOptimizer<'_, E> {
Expand All @@ -20,22 +20,23 @@ impl<E: Equation> CostFunction for SppOptimizer<'_, E> {
fn cost(&self, spp: &Self::Param) -> Result<Self::Output, Error> {
let theta = Array1::from(spp.clone()).insert_axis(Axis(0));

let psi = psi(self.equation, self.data, &theta, self.sig, false, false)?;
let log_psi = psi(self.equation, self.data, &theta, self.sig, false, false)?;

if psi.ncols() > 1 {
tracing::error!("Psi in SppOptimizer has more than one column");
if log_psi.ncols() > 1 {
tracing::error!("log_psi in SppOptimizer has more than one column");
}
if psi.nrows() != self.pyl.len() {
if log_psi.nrows() != self.log_pyl.len() {
tracing::error!(
"Psi in SppOptimizer has {} rows, but spp has {}",
psi.nrows(),
self.pyl.len()
"log_psi in SppOptimizer has {} rows, but spp has {}",
log_psi.nrows(),
self.log_pyl.len()
);
}
let nsub = psi.nrows() as f64;
let nsub = log_psi.nrows() as f64;
let mut sum = -nsub;
for (p_i, pyl_i) in psi.iter().zip(self.pyl.iter()) {
sum += p_i / pyl_i;
// Convert log-likelihoods back to likelihood ratio: exp(log_psi - log_pyl)
for (log_p_i, log_pyl_i) in log_psi.iter().zip(self.log_pyl.iter()) {
sum += (log_p_i - log_pyl_i).exp();
}
Ok(-sum)
}
Expand All @@ -46,13 +47,13 @@ impl<'a, E: Equation> SppOptimizer<'a, E> {
equation: &'a E,
data: &'a Data,
sig: &'a ErrorModels,
pyl: &'a Array1<f64>,
log_pyl: &'a Array1<f64>,
) -> Self {
Self {
equation,
data,
sig,
pyl,
log_pyl,
}
}
pub fn optimize_point(self, spp: Array1<f64>) -> Result<Array1<f64>, Error> {
Expand Down
27 changes: 1 addition & 26 deletions src/simulator/equation/analytical/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -191,7 +191,7 @@ impl EquationPriv for Analytical {
let pred = y[observation.outeq()];
let pred = observation.to_prediction(pred, x.as_slice().to_vec());
if let Some(error_models) = error_models {
likelihood.push(pred.likelihood(error_models)?);
likelihood.push(pred.log_likelihood(error_models)?);
}
output.add_prediction(pred);
Ok(())
Expand Down Expand Up @@ -283,16 +283,6 @@ mod tests {
}
}
impl Equation for Analytical {
fn estimate_likelihood(
&self,
subject: &Subject,
support_point: &Vec<f64>,
error_models: &ErrorModels,
cache: bool,
) -> Result<f64, PharmsolError> {
_estimate_likelihood(self, subject, support_point, error_models, cache)
}

fn estimate_log_likelihood(
&self,
subject: &Subject,
Expand Down Expand Up @@ -341,18 +331,3 @@ fn _subject_predictions(
) -> Result<SubjectPredictions, PharmsolError> {
Ok(ode.simulate_subject(subject, support_point, None)?.0)
}

fn _estimate_likelihood(
ode: &Analytical,
subject: &Subject,
support_point: &Vec<f64>,
error_models: &ErrorModels,
cache: bool,
) -> Result<f64, PharmsolError> {
let ypred = if cache {
_subject_predictions(ode, subject, support_point)
} else {
_subject_predictions_no_cache(ode, subject, support_point)
}?;
ypred.likelihood(error_models)
}
27 changes: 3 additions & 24 deletions src/simulator/equation/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -171,32 +171,10 @@ pub(crate) trait EquationPriv: EquationTypes {
/// and estimate parameters.
#[allow(private_bounds)]
pub trait Equation: EquationPriv + 'static + Clone + Sync {
/// Estimate the likelihood of the subject given the support point and error model.
///
/// This function calculates how likely the observed data is given the model
/// parameters and error model. It may use caching for performance.
///
/// # Parameters
/// - `subject`: The subject data
/// - `support_point`: The parameter values
/// - `error_model`: The error model
/// - `cache`: Whether to use caching
///
/// # Returns
/// The likelihood value (product of individual observation likelihoods)
fn estimate_likelihood(
&self,
subject: &Subject,
support_point: &Vec<f64>,
error_models: &ErrorModels,
cache: bool,
) -> Result<f64, PharmsolError>;

/// Estimate the log-likelihood of the subject given the support point and error model.
///
/// This function calculates the log of how likely the observed data is given the model
/// parameters and error model. It is numerically more stable than `estimate_likelihood`
/// for extreme values or many observations.
/// parameters and error model. It is numerically stable for extreme values or many observations.
///
/// # Parameters
/// - `subject`: The subject data
Expand Down Expand Up @@ -282,7 +260,8 @@ pub trait Equation: EquationPriv + 'static + Clone + Sync {
)?;
}
}
let ll = error_models.map(|_| likelihood.iter().product::<f64>());
// Return log-likelihood (sum of log-likelihoods) when error_models is provided
Copy link

Copilot AI Dec 30, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The inline comment correctly describes returning log-likelihood, which matches the implementation. However, the function's documentation comment at line 231 (outside the changed region) still says "optional likelihood" and should be updated in a follow-up to say "optional log-likelihood" for consistency.

Copilot uses AI. Check for mistakes.
let ll = error_models.map(|_| likelihood.iter().sum::<f64>());
Ok((output, ll))
}
}
Expand Down
29 changes: 2 additions & 27 deletions src/simulator/equation/ode/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,21 +69,6 @@ fn spphash(spp: &[f64]) -> u64 {

/// Hash a subject ID string to u64 for cache key generation.

fn _estimate_likelihood(
ode: &ODE,
subject: &Subject,
support_point: &Vec<f64>,
error_models: &ErrorModels,
cache: bool,
) -> Result<f64, PharmsolError> {
let ypred = if cache {
_subject_predictions(ode, subject, support_point)
} else {
_subject_predictions_no_cache(ode, subject, support_point)
}?;
ypred.likelihood(error_models)
}

#[inline(always)]
#[cached(
ty = "UnboundCache<(u64, u64), SubjectPredictions>",
Expand Down Expand Up @@ -174,16 +159,6 @@ impl EquationPriv for ODE {
}

impl Equation for ODE {
fn estimate_likelihood(
&self,
subject: &Subject,
support_point: &Vec<f64>,
error_models: &ErrorModels,
cache: bool,
) -> Result<f64, PharmsolError> {
_estimate_likelihood(self, subject, support_point, error_models, cache)
}

fn estimate_log_likelihood(
&self,
subject: &Subject,
Expand Down Expand Up @@ -324,7 +299,7 @@ impl Equation for ODE {
let pred =
observation.to_prediction(pred, solver.state().y.as_slice().to_vec());
if let Some(error_models) = error_models {
likelihood.push(pred.likelihood(error_models)?);
likelihood.push(pred.log_likelihood(error_models)?);
}
output.add_prediction(pred);
}
Expand Down Expand Up @@ -379,7 +354,7 @@ impl Equation for ODE {
}
}
}
let ll = error_models.map(|_| likelihood.iter().product::<f64>());
let ll = error_models.map(|_| likelihood.iter().sum::<f64>());
Ok((output, ll))
}
}
51 changes: 24 additions & 27 deletions src/simulator/equation/sde/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -306,23 +306,31 @@ impl EquationPriv for SDE {
//e = y[t] .- x[:,1]
// q = pdf.(Distributions.Normal(0, 0.5), e)
if let Some(em) = error_models {
let mut q: Vec<f64> = Vec::with_capacity(self.nparticles);
// Compute log-likelihoods for each particle
let mut log_q: Vec<f64> = Vec::with_capacity(self.nparticles);

pred.iter().for_each(|p| {
let lik = p.likelihood(em);
match lik {
Ok(l) => q.push(l),
Err(e) => panic!("Error in likelihood calculation: {:?}", e),
let log_lik = p.log_likelihood(em);
match log_lik {
Ok(l) => log_q.push(l),
Err(e) => panic!("Error in log-likelihood calculation: {:?}", e),
}
});
let sum_q: f64 = q.iter().sum();
let w: Vec<f64> = q.iter().map(|qi| qi / sum_q).collect();

// Use log-sum-exp trick for numerical stability
let max_log_q = log_q.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
let sum_exp: f64 = log_q.iter().map(|&lq| (lq - max_log_q).exp()).sum();
let log_sum_q = max_log_q + sum_exp.ln();

// Compute normalized weights from log-likelihoods
let w: Vec<f64> = log_q.iter().map(|&lq| (lq - log_sum_q).exp()).collect();
let i = sysresample(&w);
let a: Vec<DVector<f64>> = i.iter().map(|&i| x[i].clone()).collect();
*x = a;
likelihood.push(sum_q / self.nparticles as f64);
// let qq: Vec<f64> = i.iter().map(|&i| q[i]).collect();
// likelihood.push(qq.iter().sum::<f64>() / self.nparticles as f64);

// Push the average likelihood (in regular space) for final computation
// log(mean(likelihood)) = log(sum(exp(log_lik))) - log(n)
likelihood.push((log_sum_q - (self.nparticles as f64).ln()).exp());
Comment on lines +331 to +333
Copy link

Copilot AI Dec 30, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The SDE implementation is computing regular likelihoods (not log-likelihoods) in the particle filter, but then these values are being summed in the default simulate_subject method (line 264 in mod.rs). This is mathematically incorrect. Likelihoods should be multiplied (product) across observations, not summed. Since you're working in regular likelihood space here, you should either: (1) change line 333 to push log-likelihoods and keep the sum, or (2) override simulate_subject for SDE to compute the product instead of sum, then take the log in estimate_log_likelihood. The first approach is recommended for consistency with ODE and Analytical implementations.

Suggested change
// Push the average likelihood (in regular space) for final computation
// log(mean(likelihood)) = log(sum(exp(log_lik))) - log(n)
likelihood.push((log_sum_q - (self.nparticles as f64).ln()).exp());
// Push the average log-likelihood for final computation
// log(mean(likelihood)) = log(sum(exp(log_lik))) - log(n)
likelihood.push(log_sum_q - (self.nparticles as f64).ln());

Copilot uses AI. Check for mistakes.
}
Ok(())
}
Expand Down Expand Up @@ -351,7 +359,7 @@ impl EquationPriv for SDE {
}

impl Equation for SDE {
/// Estimates the likelihood of observed data given a model and parameters.
/// Estimates the log-likelihood of observed data given a model and parameters.
///
/// # Arguments
///
Expand All @@ -363,31 +371,20 @@ impl Equation for SDE {
/// # Returns
///
/// The log-likelihood of the observed data given the model and parameters.
fn estimate_likelihood(
fn estimate_log_likelihood(
&self,
subject: &Subject,
support_point: &Vec<f64>,
error_models: &ErrorModels,
cache: bool,
) -> Result<f64, PharmsolError> {
if cache {
// For SDE, the particle filter computes likelihood in regular space internally.
// We take the log of the final likelihood.
let lik = if cache {
_estimate_likelihood(self, subject, support_point, error_models)
} else {
_estimate_likelihood_no_cache(self, subject, support_point, error_models)
Copy link

Copilot AI Dec 30, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The function _estimate_likelihood_no_cache is called but never defined in this module or elsewhere in the codebase. This will cause a compilation error. You need to either define this function or change the implementation to use simulate_subject directly without caching, similar to how it's done for the cached version.

Copilot uses AI. Check for mistakes.
}
}

fn estimate_log_likelihood(
&self,
subject: &Subject,
support_point: &Vec<f64>,
error_models: &ErrorModels,
cache: bool,
) -> Result<f64, PharmsolError> {
// For SDE, the particle filter computes likelihood in regular space.
// We take the log of the cached/computed likelihood.
// Note: For extreme underflow cases, this may return -inf.
let lik = self.estimate_likelihood(subject, support_point, error_models, cache)?;
}?;
if lik > 0.0 {
Ok(lik.ln())
} else {
Expand Down
Loading
Loading