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
This PR integrates the new BestDose workflow into the Pmetrics codebase, adding Rust-backed posterior preparation/optimization endpoints plus R-facing R6 objects (bd, bd_post), reporting, and supporting UX/options updates.
Changes:
- Add Rust extendr exports and a new
bestdose_executormodule to compute/posterior + optimize doses via PMcore’s BestDose API. - Add R BestDose R6 classes, plotting, and an HTML report (RMarkdown template + report builder utilities).
- Add user option support for platform-specific date parsing (
date_format) and helper functions to check/download the latest R release.
Reviewed changes
Copilot reviewed 36 out of 43 changed files in this pull request and generated 11 comments.
Show a summary per file
| File | Description |
|---|---|
| src/rust/src/lib.rs | Registers new extendr exports for BestDose and adjusts default log level. |
| src/rust/src/executor.rs | Minor fit flow tweak (mutable result). |
| src/rust/src/bestdose_executor.rs | New Rust BestDose executor: posterior handle + optimization + R conversions. |
| src/rust/Cargo.toml | Updates Rust deps (PMcore source + libloading). |
| R/zzz.R | Adds interactive attach checks for latest Pmetrics and platform-specific latest R. |
| R/PMutilities.R | Adds latestR() and downloadR() helpers. |
| R/PMoptions.R | Adds date_format option initialization + Shiny UI control. |
| R/PM_sim.R | Adds validation of unrecognized ... args; sim logic tweaks. |
| R/PM_post.R | Adjusts post-processing of pred.csv columns. |
| R/PM_model.R | Adds quiet support to compile/run/sim; changes compile invocation behavior. |
| R/PM_data.R | Date parsing improvements (PMoptions-aware) and plot pred handling tweaks. |
| R/PM_bestdose.R | New BestDose R6 API (bd, bd_post) + plotting method. |
| R/PM_bestdose_report.R | New BestDose report builder + metrics/plots helpers. |
| R/extendr-wrappers.R | Adds R .Call() wrappers for new Rust BestDose exports. |
| NAMESPACE | Exports bd, bd_post, latestR, downloadR; registers plot.bd. |
| man/summary.PM_sim.Rd | Updates docs to use getPMoptions("digits") + group arg naming. |
| man/PMcheck.Rd | Updates PMcheck() signature default for path. |
| man/PM_model.Rd | Documents new quiet params and updated method signatures. |
| man/plot.bd.Rd | New documentation for plot.bd. |
| man/latestR.Rd | New documentation for latestR(). |
| man/downloadR.Rd | New documentation for downloadR(). |
| man/bd.Rd | New documentation for bd. |
| man/bd_post.Rd | New documentation for bd_post. |
| inst/report/templates/bestdose.Rmd | New RMarkdown template for BestDose HTML reports. |
| inst/options/PMoptions.json | Adds date_format default. |
| inst/Examples/src/bestdose_target.csv | New BestDose example target dataset. |
| inst/Examples/src/bestdose_target_pmcore.csv | New BestDose example target dataset (PMcore variant). |
| inst/Examples/src/bestdose_prior.csv | New BestDose example prior file. |
| inst/Examples/src/bestdose_prior_pmcore.csv | New BestDose example prior file (PMcore variant). |
| inst/Examples/src/bestdose_past.csv | New BestDose example past data file. |
| inst/Examples/Rscript/test_status_column.R | New example/test script for report dose status handling. |
| inst/Examples/Rscript/test_future_vectors.R | New example/test script for vectorized future specs. |
| inst/Examples/Rscript/run_bestdose_comparison.R | New example script comparing bias/prior weights. |
| inst/Examples/Rscript/bestdose_simple_test.R | New example script for basic BestDose usage. |
| DESCRIPTION | Updates package metadata (version/URLs). |
| .gitignore | Ensures inst/Examples is tracked while keeping other example outputs ignored. |
| _pkgdown.yml | Adds BestDose + latestR/downloadR to pkgdown reference. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| let (library, (eq, meta)) = | ||
| 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))?; | ||
|
|
There was a problem hiding this comment.
settings(..., "/tmp/bestdose") hard-codes a Unix-only output path. This will fail on Windows and can also violate sandboxed environments. Use a platform-independent temp directory (e.g., std::env::temp_dir() / tempfile) or pass an output path from R/settings.
| Self { | ||
| id: id.to_string(), | ||
| time: pred.time(), | ||
| observed: pred.obs().unwrap_or(0.0), | ||
| pop_mean: pred.pop_mean(), | ||
| pop_median: pred.pop_median(), | ||
| post_mean: pred.post_mean(), | ||
| post_median: pred.post_median(), | ||
| outeq: pred.outeq(), |
There was a problem hiding this comment.
Missing observations are converted to 0.0 (unwrap_or(0.0)), which will be indistinguishable from a real zero observation and can bias plots/metrics. Preserve missingness by emitting an R NA (e.g., use Rfloat/Option<f64> mapping) or a dedicated missing sentinel handled on the R side.
There was a problem hiding this comment.
This is a good catch - should be Option<f64>
| // Predictions as data frame | ||
| let pred_rows: Vec<BestDosePredictionRow> = result | ||
| .predictions() | ||
| .predictions() | ||
| .iter() | ||
| .map(|p| BestDosePredictionRow::from_np_prediction(p, "subject_1")) | ||
| .collect(); |
There was a problem hiding this comment.
Prediction rows are tagged with a hard-coded id ("subject_1"). If input data uses a different id (or multiple subjects), downstream joins/plots will be wrong. Use the actual subject id from the parsed target/past data when constructing BestDosePredictionRow.
| bestdose_prepare <- function(model_path, prior_path, past_data_path, params, kind) .Call(wrap__bestdose_prepare, model_path, prior_path, past_data_path, params, kind) | ||
|
|
||
| bestdose_optimize <- function(handle, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type) .Call(wrap__bestdose_optimize, handle, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type) | ||
|
|
||
|
|
There was a problem hiding this comment.
This file is marked "Generated by extendr: Do not edit by hand", but new wrappers were added manually. To avoid future regeneration overwriting changes, consider regenerating wrappers via rextendr::document() after adding the Rust #[extendr] exports, and keep this file fully generated.
| bestdose_prepare <- function(model_path, prior_path, past_data_path, params, kind) .Call(wrap__bestdose_prepare, model_path, prior_path, past_data_path, params, kind) | |
| bestdose_optimize <- function(handle, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type) .Call(wrap__bestdose_optimize, handle, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type) |
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 20 out of 25 changed files in this pull request and generated 8 comments.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| pub fn from_np_prediction( | ||
| pred: &pmcore::routines::output::predictions::NPPredictionRow, | ||
| id: &str, | ||
| ) -> Self { | ||
| Self { | ||
| id: id.to_string(), | ||
| time: pred.time(), | ||
| observed: pred.obs().unwrap_or(0.0), | ||
| pop_mean: pred.pop_mean(), | ||
| pop_median: pred.pop_median(), | ||
| post_mean: pred.post_mean(), | ||
| post_median: pred.post_median(), | ||
| outeq: pred.outeq(), | ||
| } |
There was a problem hiding this comment.
BestDosePredictionRow::from_np_prediction() converts missing observations to 0.0 via unwrap_or(0.0), which makes missing indistinguishable from a true 0 and will bias plots/metrics. Prefer preserving missingness (e.g., store observed as Option<f64> and map to R NA_real_).
| ) -> std::result::Result<Self, String> { | ||
| let (library, (eq, meta)) = | ||
| 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))?; | ||
|
|
There was a problem hiding this comment.
BestDosePosteriorHandle::new() hard-codes the settings output path to /tmp/bestdose. This is non-portable (e.g., Windows) and can fail in restricted environments. Prefer using a temp dir derived from the R session (e.g., R tempdir() passed down) or std::env::temp_dir() and ensure the directory exists before writing settings/logs.
| if (has_past) { | ||
| max(x$past$standard_data$time, na.rm = TRUE) | ||
| } else { | ||
| x$start_offset %||% 0 |
There was a problem hiding this comment.
This file uses the %||% infix operator (e.g., x$start_offset %||% 0), but %||% is not defined/imported in the package namespace. This will error at runtime unless %||% is provided. Consider using rlang::%||% (and import it) or defining a package-level %||% helper.
| x$start_offset %||% 0 | |
| if (is.null(x$start_offset)) 0 else x$start_offset |
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).
Update BestDose objects and plots