feat: BestDose (BD-Michael) - #315
Conversation
Add BestDose dose optimization feature ported from Pmetrics_rust bestdose branch: - R/PM_bestdose.R: PM_bestdose and PM_bestdose_problem R6 classes for Bayesian dose optimization with support for concentration and AUC targets - src/rust/src/bestdose_executor.rs: Rust backend for BestDose optimization using pmcore's BestDoseProblem with ODE model support - Updated lib.rs with bestdose, bestdose_prepare, bestdose_optimize exports - Updated extendr-wrappers.R with R-side wrapper functions - Updated NAMESPACE with PM_bestdose, PM_bestdose_problem, bestdose exports - Bumped pmcore dependency from 0.21.1 to 0.22.1 (required for bestdose) - Added libloading dependency for dynamic model loading - Added bestdose example data (past, prior, target CSVs) and test script - Fixed executor.rs mutability issue for pmcore 0.22.1 compatibility
Reworked bd and bd_post R6 classes: bd_post$optimize() now returns a bd object, and bd$new() supports one-shot usage. Removed dead Rust code (bestdose(), bestdose_ode(), bestdose_analytical()). Fixed bias_weight being NULL and print() breaking with variable dose counts. Updated PMcore dep to use relative time_offset gap semantics.
- Combine clean one-shot/two-stage API with inline future=list(...) support - Add top-5 weighted plot curves, triangle dose indicators, plot() method - Fix .sim_future() dose recycling when past doses are present
…sed most of the logs for bestdose, needs more work
…l some isues with the dose markers
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 83 out of 92 changed files in this pull request and generated 4 comments.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
|
|
||
| validate_paths(data_path, model_path); | ||
| let theta = parse_theta(theta); | ||
| let data = read_pmetrics(data_path).expect("Failed to parse data"); | ||
| let data = read_pmetrics_for_kind(data_path, kind).expect("Failed to parse data"); | ||
| let subjects = data.subjects(); |
There was a problem hiding this comment.
simulate_all() returns Result<...>, but it calls read_pmetrics_for_kind(...).expect(...), which will panic instead of returning an error to R. Please propagate the error (use ?) like simulate_one() does, so callers get a proper R error rather than a crash.
| id: id.to_string(), | ||
| time: pred.time(), | ||
| observed: pred.obs().unwrap_or(0.0), | ||
| pop_mean: pred.pop_mean(), | ||
| pop_median: pred.pop_median(), |
There was a problem hiding this comment.
Missing observations are currently coerced to 0.0 via unwrap_or(0.0), which is indistinguishable from a true zero and will bias downstream summaries/plots. Preserve missingness by emitting an NA-compatible value (e.g., store observed as Option<f64>/Rfloat and map None to NA_REAL on the R side).
| unsafe { pmcore::prelude::pharmsol::exa::load::load::<ODE>(model_path) }; | ||
|
|
||
| let settings = settings(params, meta.get_params(), "/tmp/bestdose") | ||
| .map_err(|e| format!("Failed to parse settings: {}", e))?; | ||
|
|
||
| let (population_theta, prior_weights) = |
There was a problem hiding this comment.
The settings initialization uses a hard-coded output path ("/tmp/bestdose"). This will fail or behave unexpectedly on non-Unix platforms (notably Windows) and can cause collisions between concurrent runs. Use a per-run temporary directory (e.g., std::env::temp_dir() + unique subdir) or accept an output path from the R side.
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 21 out of 22 changed files in this pull request and generated 5 comments.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| cond <- rust_args[[1]] | ||
| then_code <- expr_to_rust(args[[2]], params, covs, declared) | ||
| then_code <- expr_to_rust(args[[2]], params, covs, declared, index_mode = index_mode) | ||
| if (length(args) == 3) { | ||
| else_code <- expr_to_rust(args[[3]], params, covs, declared) | ||
| else_code <- expr_to_rust(args[[3]], params, covs, declared, index_mode = index_mode) |
There was a problem hiding this comment.
expr_to_rust() now calls itself with index_mode = index_mode in the if branch, but expr_to_rust() has no index_mode parameter and there is no index_mode variable in scope. This will error at runtime when transpiling models. Either add index_mode to the function signature and propagate it through all recursive calls, or remove these named arguments if they aren't needed.
| # Print target/future data for debugging | ||
| if (!quiet) { | ||
| cli::cli_h3("Target/Future data (for debugging)") | ||
| if (!is.null(x$future) && !is.null(x$future$standard_data)) { | ||
| print(x$future$standard_data) | ||
| } else { | ||
| cli::cli_alert_info("No future data available.") | ||
| } | ||
| } |
There was a problem hiding this comment.
plot.bd() prints the entire future dataset under a "for debugging" banner whenever quiet = FALSE. This will be very noisy in normal interactive use and in vignettes/reports. Consider removing this debug output or gating it behind an explicit debug flag / option so quiet = FALSE only controls informational messages.
| if (is.null(future$dose) || is.null(future$target)) { | ||
| cli::cli_abort(c( | ||
| "x" = "Future setup must include both {.code dose} and {.code target}.", | ||
| "i" = "Provide {.code future = list(dose = ..., target = ...)}." | ||
| )) | ||
| } | ||
|
|
||
| # Apply defaults | ||
| future$target_type <- future$target_type %||% "concentration" | ||
| future$frequency <- future$frequency %||% 24 | ||
| future$route <- future$route %||% 0 | ||
| future$number <- future$number %||% 1 | ||
|
|
||
| # Determine effective number from vector lengths of dose, frequency, route | ||
| n_from_vectors <- max( | ||
| length(future$dose), | ||
| length(future$frequency), | ||
| length(future$route) | ||
| ) | ||
| if (future$number < n_from_vectors) { | ||
| cli::cli_inform(c( | ||
| "!" = "{.arg number} ({future$number}) is less than the maximum length of | ||
| {.code dose}, {.code frequency}, or {.code route} ({n_from_vectors}).", | ||
| ">" = "Setting {.arg number} to {n_from_vectors}." | ||
| )) | ||
| future$number <- n_from_vectors | ||
| } | ||
| n <- as.integer(future$number) | ||
|
|
||
| # Recycle the last element of x to reach length n | ||
| recycle_last <- function(x, n) { | ||
| len <- length(x) | ||
| if (len >= n) { | ||
| return(x[seq_len(n)]) | ||
| } | ||
| c(x, rep(x[len], n - len)) | ||
| } |
There was a problem hiding this comment.
.build_future_data() only checks is.null(future$dose)/is.null(future$frequency)/is.null(future$route), but not zero-length inputs. If a user passes dose = numeric(0) (or empty frequency/route), recycle_last() will index x[len] with len = 0 and error. Add validation that these vectors are numeric and have length >= 1 after defaults are applied.
| fit <- actual |> | ||
| dplyr::left_join(pred, by = dplyr::join_by(id, time, outeq)) |> | ||
| dplyr::mutate( | ||
| obsSD = dplyr::coalesce(c0, 0) + dplyr::coalesce(c1, 0) * obs + dplyr::coalesce(c2, 0) * obs^2 + dplyr::coalesce(c3, 0) * obs^3, | ||
| obsSD = dplyr::if_else(!is.finite(obsSD) | obsSD <= 0, 1, obsSD), | ||
| d = pred - obs, | ||
| ds = d * d, | ||
| wd = d / obsSD, | ||
| wds = wd * wd, | ||
| kind = label | ||
| ) |> | ||
| dplyr::filter(!is.na(obs), !is.na(pred)) |
There was a problem hiding this comment.
bd_fit_data() defines d = pred - obs, but bd_fit_metrics() (and the package-wide PM_op summary) treat d as an absolute error when computing MAE/MBASE/etc. This will flip the meaning of MAE/bias in the BestDose report. Use d = abs(pred - obs) (and similarly define wd from that) to match the existing metric definitions used elsewhere in the package.
| future_predictions <- x$result$predictions | ||
| if (!is.null(future_predictions) && nrow(future_predictions) > 0) { | ||
| future_predictions <- future_predictions |> | ||
| dplyr::filter(outeq == !!(outeq - 1)) |> |
There was a problem hiding this comment.
bd_report_build() filters Rust future_predictions with outeq == outeq - 1, which will drop all predictions when outeq is the usual 1-based value used across Pmetrics. Unless Rust is intentionally emitting 0-based outeq (which would conflict with other Rust simulation exports), this should filter by outeq == outeq to keep the requested output equation.
| dplyr::filter(outeq == !!(outeq - 1)) |> | |
| dplyr::filter(outeq == !!outeq) |> |
mhovd
changed the title
Introduce a new simulate argument (default FALSE) to the bd R6 class constructor in R/PM_bestdose.R. The automatic calls to private$.sim_past() and private$.sim_future() are now executed only when simulate is TRUE, and a cli::cli_alert is shown when running simulations. This preserves previous behavior by default while allowing callers to opt out of immediate simulation (useful for performance or manual control).
|
@R-garreau Is there something in this branch that you still need? I can try to fix the conflicts if that is the case |
|
@Siel NB : On the bright side, I installed the desktop (tauri) app for the first time today (i hope i will have translated bestdose by the end of the year) |
Should we try to rebase that branch on main, and then complete the merge?
@R-garreau , when you say "better", in what way? Does that mean we should not trust the Shiny app dose recs? Is there an issue with Pmetrics? |
|
I can work on merging it. It would require updating rebasing it from main and fix all the issues that will arise, but it is doable. @R-garreau I would love to check the tauri app you're working on. I also started experimenting with bestdose and tauri but it is just on its infancy and nothing is really working yet |
|
I'm already working on the rebase with main and testing before merging! :) |
|
I'm creating a BD-Michael2 which is rebased on main, just so we don't mess anything up. I think we can close this PR and we'll eventually open a PR to merge BD-Michael2 into main, but only after we are satisfied that nothing in main will be overwritten that we don't want to be. |
|
Moved to #374 |
|
@mnneely No it just that this branch of Pmetrics was not patch with the proper posterior calculation (weight were all 0 and 1). @Siel you can pull BestDose_Rust. or i can send you the .exe but there is still so much to do ! |
Update BestDose objects and plots