Conversation
|
This PR is a draft, as it requires the new changes to pharmsol |
|
Bencher Report| Branch | mmopt |
| Testbed | pmcore-runner |
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 | 14.76 s(-0.85%)Baseline: 14.88 s | 25.35 s (58.23%) |
| bimodal_ke_npod | 📈 view plot 🚷 view threshold | 12.08 s(+22.21%)Baseline: 9.88 s | 13.59 s (88.84%) |
| bimodal_ke_postprob | 📈 view plot 🚷 view threshold | 3.70 s(+9.81%)Baseline: 3.37 s | 3.96 s (93.37%) |
|
| Branch | mmopt |
| Testbed | mhovd-pgx |
Click to view all benchmark results
| Benchmark | Latency | Benchmark Result milliseconds (ms) (Result Δ%) | Upper Boundary milliseconds (ms) (Limit %) |
|---|---|---|---|
| bimodal_ke_npag | 📈 view plot 🚷 view threshold | 4,964.30 ms(+0.48%)Baseline: 4,940.55 ms | 5,168.07 ms (96.06%) |
| bimodal_ke_npod | 📈 view plot 🚷 view threshold | 1,406.10 ms(0.00%)Baseline: 1,406.10 ms | 3,865.83 ms (36.37%) |
| bimodal_ke_postprob | 📈 view plot 🚷 view threshold | 147.54 ms(-65.52%)Baseline: 427.85 ms | 1,862.47 ms (7.92%) |
There was a problem hiding this comment.
Pull request overview
Adds a new mmopt module to compute Bayes-risk–based optimal sampling times (MMopt), plus examples and tests to validate core behavior and input validation.
Changes:
- Introduces
mmoptoptimization andbayes_riskscoring APIs withMmoptResult. - Adds integration tests covering common scenarios and validation failures.
- Adds runnable examples replicating (parts of) Bayard & Neely (2017) experiments.
Reviewed changes
Copilot reviewed 5 out of 5 changed files in this pull request and generated 6 comments.
Show a summary per file
| File | Description |
|---|---|
src/mmopt/mod.rs |
Implements MMopt search, Bayes risk calculation, and internal helpers/tests. |
src/lib.rs |
Exposes the new mmopt module and re-exports it in the prelude. |
tests/mmopt_tests.rs |
Adds integration tests for mmopt behavior and validation errors. |
examples/mmopt.rs |
Example driver for paper replication (Section 4/6). |
examples/bayes_risk.rs |
Example showing how to score a fixed design via bayes_risk. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| return Err(anyhow::anyhow!( | ||
| "C({}, {}) = {} exceeds the maximum allowed combinations ({}). \ | ||
| Reduce the number of candidate times or increase nsamp.", | ||
| times.len(), | ||
| nsamp, | ||
| n_combinations, | ||
| MAX_COMBINATIONS | ||
| )); |
There was a problem hiding this comment.
The guidance in this error message is mathematically incorrect: increasing nsamp does not necessarily reduce C(m, nsamp) (it increases until nsamp reaches ~m/2). Consider rewording to avoid suggesting “increase nsamp” as a general remedy (e.g., only recommend reducing candidate times, or mention choosing nsamp farther from m/2).
| .par_iter() | ||
| .map(|combo| { | ||
| let risk = calculate_risk(combo, &pred_matrix, &errormodel, &weights_vec); | ||
| (combo.clone(), risk) |
There was a problem hiding this comment.
This eagerly allocates all C(m, nsamp) combinations into memory and then clones each combo again inside the parallel map. Even with MAX_COMBINATIONS = 1_000_000, this can be very memory-heavy and slow. Consider iterating combinations lazily (streaming generator) and avoiding per-iteration cloning (e.g., keep a reference during reduction and only copy the best combo once at the end).
| (combo.clone(), risk) | |
| (combo, risk) |
| } | ||
|
|
||
| let weights_vec = weights.to_vec(); | ||
|
|
There was a problem hiding this comment.
The docstrings for mmopt/bayes_risk state that weights “must sum to ~1.0”, but the implementation only checks length. To keep the API contract consistent (and avoid returning Bayes risks on an invalid probability distribution), consider validating that weights are non-negative and sum to 1 within a tolerance, or relax/update the docs if normalization is intentionally left to callers.
| // Validate that weights form a proper probability distribution | |
| if weights_vec.iter().any(|w| *w < 0.0) { | |
| return Err(anyhow::anyhow!( | |
| "Weights must be non-negative, but found a negative value" | |
| )); | |
| } | |
| let weight_sum: f64 = weights_vec.iter().sum(); | |
| let tol: f64 = 1e-6; | |
| if (weight_sum - 1.0).abs() > tol { | |
| return Err(anyhow::anyhow!( | |
| "Weights must sum to 1.0 within a tolerance of {}, but sum to {}", | |
| tol, | |
| weight_sum | |
| )); | |
| } |
| // Only proceed if the subject actually has multiple occasions | ||
| if subject.occasions().len() > 1 { | ||
| let result = mmopt( | ||
| &theta, | ||
| &subject, | ||
| eq, | ||
| additive_error_model(), | ||
| 0, | ||
| 1, | ||
| &Weights::from_vec(vec![0.5, 0.5]), | ||
| ); | ||
| assert!(result.is_err()); | ||
| } |
There was a problem hiding this comment.
This if subject.occasions().len() > 1 guard can cause the test to silently pass without asserting anything if the builder behavior changes (or if occasions are unexpectedly collapsed). It’s safer to assert the precondition (e.g., assert!(... > 1)) and then always check that mmopt returns an error.
| // Only proceed if the subject actually has multiple occasions | |
| if subject.occasions().len() > 1 { | |
| let result = mmopt( | |
| &theta, | |
| &subject, | |
| eq, | |
| additive_error_model(), | |
| 0, | |
| 1, | |
| &Weights::from_vec(vec![0.5, 0.5]), | |
| ); | |
| assert!(result.is_err()); | |
| } | |
| // Assert the precondition that the subject has multiple occasions | |
| assert!(subject.occasions().len() > 1); | |
| let result = mmopt( | |
| &theta, | |
| &subject, | |
| eq, | |
| additive_error_model(), | |
| 0, | |
| 1, | |
| &Weights::from_vec(vec![0.5, 0.5]), | |
| ); | |
| assert!(result.is_err()); |
| /// Paper results (Table 6.2): | ||
| /// n=1: t* = {4.25}, Bayes Risk = 0.5474 | ||
| /// n=2: t* = {1.0, 9.5}, Bayes Risk = 0.2947 | ||
| /// n=3: t* = {1.0, 1.0, 10.5}, Bayes Risk = 0.2325 |
There was a problem hiding this comment.
The paper’s reported 3-sample design includes a repeated time (two samples at 1.0). The current mmopt implementation searches combinations without repetition, so it can never return duplicate times; this example text may therefore set expectations that can’t be met. Consider clarifying this limitation here (and/or updating mmopt to support repeated sampling times if that’s required for replication).
| /// Paper results (Table 6.2): | |
| /// n=1: t* = {4.25}, Bayes Risk = 0.5474 | |
| /// n=2: t* = {1.0, 9.5}, Bayes Risk = 0.2947 | |
| /// n=3: t* = {1.0, 1.0, 10.5}, Bayes Risk = 0.2325 | |
| /// Paper results (Table 6.2, allowing repeated sampling times): | |
| /// n=1: t* = {4.25}, Bayes Risk = 0.5474 | |
| /// n=2: t* = {1.0, 9.5}, Bayes Risk = 0.2947 | |
| /// n=3: t* = {1.0, 1.0, 10.5}, Bayes Risk = 0.2325 | |
| /// | |
| /// Note: The current `mmopt` implementation searches over combinations of | |
| /// distinct sampling times (no repetition). As a result, designs with | |
| /// repeated times like the paper's n=3 optimum {1.0, 1.0, 10.5} cannot be | |
| /// reproduced exactly here; instead, we obtain the best 3-point design with | |
| /// all sampling times distinct. |
| @@ -0,0 +1,502 @@ | |||
| use anyhow::{Ok, Result}; | |||
There was a problem hiding this comment.
anyhow does not export Ok; importing it here will fail to compile. Import only anyhow::Result (or just use the Ok(...) variant without importing it).
| use anyhow::{Ok, Result}; | |
| use anyhow::Result; |
|
We need to discuss if this should be allowed to run on multiple subjects, and if the cost should be weighted |
Implements the MMopt algorithm for determining the optimal sample times