feat: Add support for log-likelihoods by Siel · Pull Request #231 · LAPKB/PMcore

Siel · 2025-11-27T02:04:18Z

No description provided.

github-actions · 2025-11-27T02:09:15Z

Bencher Report

Branch	log-likelihood
Testbed	rust-moan

Click to view all benchmark results

Benchmark	Latency	Benchmark Result seconds (s) (Result Δ%)	Upper Boundary seconds (s) (Limit %)
bimodal_ke_npag	📈 view plot 🚷 view threshold	11.93 s (-51.33%) Baseline: 24.50 s	32.42 s (36.79%)
bimodal_ke_npod	📈 view plot 🚷 view threshold	5.79 s (-43.29%) Baseline: 10.21 s	12.97 s (44.65%)
bimodal_ke_postprob	📈 view plot 🚷 view threshold	1.70 s (-49.66%) Baseline: 3.38 s	4.42 s (38.44%)

🐰 View full continuous benchmarking report in Bencher

Copilot

Pull request overview

This PR introduces log-space likelihood computations to improve numerical stability when dealing with many observations or extreme parameter values. The implementation adds a flag to track whether Psi matrices contain regular or log-space likelihoods, and provides unified dispatcher functions to automatically choose the appropriate calculation method.

Key Changes

Added is_log_space field to Psi struct with corresponding builder pattern and dispatcher functions for likelihood calculations
Implemented numerically stable logsumexp family of functions in new routines::math module
Extended IPM with log-space variant (burke_log) and unified dispatcher (burke_ipm)
Updated all algorithms (NPAG, NPOD, POSTPROB) and related code to handle both regular and log-space computations

Reviewed changes

Copilot reviewed 15 out of 15 changed files in this pull request and generated 8 comments.

Show a summary per file

File	Description
`src/structs/psi.rs`	Added `is_log_space` flag, `PsiBuilder`, `calculate_log_psi()`, and `calculate_psi_dispatch()` for unified likelihood calculations
`src/routines/settings.rs`	Added `log_space` setting (default: true) to enable log-space computations
`src/routines/output/posterior.rs`	Updated posterior calculation to handle log-space using logsumexp; removed duplicate function
`src/routines/mod.rs`	Added math module for numerical utilities
`src/routines/math.rs`	New module with `logsumexp`, `logsumexp_weighted`, and `logsumexp_rows` implementations
`src/routines/estimation/qr.rs`	Updated QR decomposition to apply softmax normalization for log-space matrices
`src/routines/estimation/ipm.rs`	Added `burke_log()` for log-space IPM and `burke_ipm()` dispatcher
`src/routines/condensation/mod.rs`	Updated to use `burke_ipm()` dispatcher
`src/bestdose/posterior.rs`	Updated to use `calculate_psi_dispatch()` and `burke_ipm()`
`src/algorithms/postprob.rs`	Updated to use `calculate_psi_dispatch()` and `burke_ipm()` with log_space setting
`src/algorithms/npod.rs`	Updated to use dispatchers and handle log-space in expansion phase
`src/algorithms/npag.rs`	Updated to use dispatchers and compute f1 correctly for log-space
`src/algorithms/mod.rs`	Updated `validate_psi()` to handle log-space (NEG_INFINITY is valid)
`Cargo.toml`	Bumped pharmsol dependency to 0.22.0 for `log_psi` function support

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copilot · 2025-11-27T02:09:48Z

src/routines/math.rs

+pub fn logsumexp_weighted(log_values: &[f64], log_weights: &[f64]) -> f64 {
+    assert_eq!(
+        log_values.len(),
+        log_weights.len(),
+        "log_values and log_weights must have the same length"
+    );
+
+    let combined: Vec<f64> = log_values
+        .iter()
+        .zip(log_weights.iter())
+        .map(|(&lv, &lw)| lv + lw)
+        .collect();
+
+    logsumexp(&combined)
+}


[nitpick] Consider adding a test case for logsumexp_weighted that verifies the behavior when weights have different magnitudes or when some log_values/log_weights are special values like NEG_INFINITY or INFINITY. This would improve test coverage and ensure edge cases are handled correctly.

Copilot · 2025-11-27T02:09:49Z

src/algorithms/npod.rs

+        let pyl = if self.psi.is_log_space() {
+            // pyl[i] = sum_j(exp(log_psi[i,j]) * w[j]) = sum_j(exp(log_psi[i,j] + log(w[j])))
+            // Using logsumexp for stability, then exp to get regular values
+            let log_w: Array1<f64> = w.iter().map(|&x| x.ln()).collect();


If any weight in w is zero or negative, x.ln() will produce NEG_INFINITY or NaN, which may lead to unexpected behavior. While weights should typically be positive (from the IPM), consider adding a validation check or documenting this assumption to prevent potential issues.

Copilot · 2025-11-27T02:09:49Z

src/routines/output/posterior.rs

+        // Calculate py[i] = sum_j(psi[i,j] * w[j]) for each subject i
+        // In log-space: py[i] = logsumexp_j(log_psi[i,j] + log(w[j]))
+        let py: Vec<f64> = if is_log_space {
+            let log_w: Vec<f64> = (0..w.len()).map(|j| w.weights().get(j).ln()).collect();


If any weight in w is zero or negative, w.weights().get(j).ln() will produce NEG_INFINITY or NaN, potentially causing numerical issues in the posterior calculation. While weights from the IPM should be positive, consider adding validation or documenting this assumption.

Copilot · 2025-11-27T02:09:49Z

src/routines/output/posterior.rs

        let posterior = Mat::from_fn(psi_matrix.nrows(), psi_matrix.ncols(), |i, j| {
-            psi_matrix.get(i, j) * w.weights().get(j) / py.get(i)
+            if is_log_space {
+                let log_w_j = w.weights().get(j).ln();


If the weight at index j is zero or negative, w.weights().get(j).ln() will produce NEG_INFINITY or NaN, potentially causing numerical issues. While weights from the IPM should be positive, consider adding validation or documenting this assumption.

Copilot · 2025-11-27T02:09:50Z

src/algorithms/npag.rs

+                    let log_w: Vec<f64> = w.weights().iter().map(|&x| x.ln()).collect();
+                    (0..psi.nrows())
+                        .map(|i| {


If any weight in w is zero or negative, .ln() will produce NEG_INFINITY or NaN, potentially causing numerical issues. While weights from the IPM should be positive, consider adding validation or documenting this assumption.

Suggested change

let log_w: Vec<f64> = w.weights().iter().map(|&x| x.ln()).collect();

(0..psi.nrows())

.map(|i| {

if w.weights().iter().any(|&x| x <= 0.0) {

bail!("All weights must be positive before taking logarithm, found zero or negative weight.");

}

let log_w: Vec<f64> = w.weights().iter().map(|&x| x.ln()).collect();

(0..psi.nrows())

Copilot · 2025-11-27T02:09:50Z

src/routines/estimation/qr.rs


+/// Compute log-sum-exp of a slice for numerical stability
+#[inline]
+fn logsumexp(values: &[f64]) -> f64 {
+    let max_val = values.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
+    if max_val.is_infinite() {
+        return max_val;
+    }
+    max_val
+        + values
+            .iter()
+            .map(|&x| (x - max_val).exp())
+            .sum::<f64>()
+            .ln()
+}
+


This logsumexp function duplicates the implementation in src/routines/math.rs. Consider importing and using crate::routines::math::logsumexp instead to maintain a single source of truth and avoid potential inconsistencies.

Suggested change

/// Compute log-sum-exp of a slice for numerical stability

#[inline]

fn logsumexp(values: &[f64]) -> f64 {

let max_val = values.iter().cloned().fold(f64::NEG_INFINITY, f64::max);

if max_val.is_infinite() {

return max_val;

}

max_val

+ values

.iter()

.map(|&x| (x - max_val).exp())

.sum::<f64>()

.ln()

}

use crate::routines::math::logsumexp;

Copilot · 2025-11-27T02:09:50Z

src/routines/estimation/ipm.rs

+    // Validate that all entries are finite
+    for row in log_psi_mat.row_iter() {
+        for &x in row.iter() {
+            if !x.is_finite() {
+                bail!("Input log-psi matrix must have finite entries");


The validation rejects NEG_INFINITY values in log-space, but NEG_INFINITY is a valid representation of zero probability in log-space (log(0) = -∞). Consider allowing NEG_INFINITY but still rejecting NaN and positive infinity. This may cause the function to fail on valid log-space inputs.

Suggested change

// Validate that all entries are finite

for row in log_psi_mat.row_iter() {

for &x in row.iter() {

if !x.is_finite() {

bail!("Input log-psi matrix must have finite entries");

// Validate that all entries are valid: allow NEG_INFINITY, but reject NaN and positive infinity

for row in log_psi_mat.row_iter() {

for &x in row.iter() {

if x.is_nan() || x == f64::INFINITY {

bail!("Input log-psi matrix must not contain NaN or positive infinity");

* wip * wip * 2 tests are failing * Fix tests * Add to_space for converting

log-likelihood

3115445

Copilot AI review requested due to automatic review settings November 27, 2025 02:04

Copilot started reviewing on behalf of Siel November 27, 2025 02:04 View session

Copilot finished reviewing on behalf of Siel November 27, 2025 02:08

Copilot AI reviewed Nov 27, 2025

View reviewed changes

mhovd and others added 2 commits December 2, 2025 15:23

Merge branch 'main' into log-likelihood

f373ee7

Formatting

bcd365b

mhovd changed the title ~~log-likelihood~~ feat: Add support for log-likelihoods Dec 11, 2025

chore: Log likelihood suggestions (#235)

7648624

* wip * wip * 2 tests are failing * Fix tests * Add to_space for converting

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

feat: Add support for log-likelihoods#231

feat: Add support for log-likelihoods#231
Siel wants to merge 4 commits intomainfrom
log-likelihood

Siel commented Nov 27, 2025

Uh oh!

github-actions bot commented Nov 27, 2025 •

edited

Loading

Uh oh!

Copilot AI left a comment

Uh oh!

Copilot AI Nov 27, 2025

Uh oh!

Copilot AI Nov 27, 2025

Uh oh!

Copilot AI Nov 27, 2025

Uh oh!

Copilot AI Nov 27, 2025

Uh oh!

Copilot AI Nov 27, 2025

Uh oh!

Copilot AI Nov 27, 2025

Uh oh!

Copilot AI Nov 27, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

Siel commented Nov 27, 2025

Uh oh!

github-actions bot commented Nov 27, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Bencher Report

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Key Changes

Reviewed changes

Uh oh!

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

Uh oh!

Copilot AI Nov 27, 2025

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

github-actions bot commented Nov 27, 2025 •

edited

Loading