diff --git a/.gitignore b/.gitignore
index d76efe5c..28eb76c2 100755
--- a/.gitignore
+++ b/.gitignore
@@ -63,7 +63,7 @@ inst/Learn/_freeze/
data-raw/Runs/
tests/testthat/Runs/
Examples
-!inst/Examples/
+inst/Examples/
Experimental
other/
1
@@ -85,3 +85,4 @@ src/Makevars.win
# Old files
Archived
+
diff --git a/Cargo.lock b/Cargo.lock
index 80e95efd..070e26b2 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -1380,9 +1380,8 @@ dependencies = [
[[package]]
name = "pmcore"
-version = "0.25.2"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "56359b47b12acc12b49b61e98c2405f8750e586a3241e8417bdfa5c59d6aa71d"
+version = "0.25.1"
+source = "git+https://github.com/LAPKB/PMcore.git?branch=fix%2Fpmetrics-analytical#8e99683e1d33fafb7c6eb56227aa353b530c48a7"
dependencies = [
"anyhow",
"argmin",
@@ -2085,9 +2084,9 @@ dependencies = [
[[package]]
name = "typenum"
-version = "1.19.0"
+version = "1.20.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "562d481066bde0658276a35467c4af00bdc6ee726305698a55b86e61d7ad82bb"
+checksum = "40ce102ab67701b8526c123c1bab5cbe42d7040ccfd0f64af1a385808d2f43de"
[[package]]
name = "ucd-trie"
@@ -2148,11 +2147,11 @@ checksum = "ccf3ec651a847eb01de73ccad15eb7d99f80485de043efb2f370cd654f4ea44b"
[[package]]
name = "wasip2"
-version = "1.0.2+wasi-0.2.9"
+version = "1.0.3+wasi-0.2.9"
source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "9517f9239f02c069db75e65f174b3da828fe5f5b945c4dd26bd25d89c03ebcf5"
+checksum = "20064672db26d7cdc89c7798c48a0fdfac8213434a1186e5ef29fd560ae223d6"
dependencies = [
- "wit-bindgen",
+ "wit-bindgen 0.57.1",
]
[[package]]
@@ -2161,7 +2160,7 @@ version = "0.4.0+wasi-0.3.0-rc-2026-01-06"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "5428f8bf88ea5ddc08faddef2ac4a67e390b88186c703ce6dbd955e1c145aca5"
dependencies = [
- "wit-bindgen",
+ "wit-bindgen 0.51.0",
]
[[package]]
@@ -2362,6 +2361,12 @@ dependencies = [
"wit-bindgen-rust-macro",
]
+[[package]]
+name = "wit-bindgen"
+version = "0.57.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1ebf944e87a7c253233ad6766e082e3cd714b5d03812acc24c318f549614536e"
+
[[package]]
name = "wit-bindgen-core"
version = "0.51.0"
diff --git a/DESCRIPTION b/DESCRIPTION
index 0f28a973..4121d777 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -77,7 +77,7 @@ Encoding: UTF-8
VignetteBuilder: knitr
LazyData: true
LazyDataCompression: xz
-Config/rextendr/version: 0.4.2
+Config/rextendr/version: 0.5.0
Config/testthat/edition: 3
Config/Needs/website:
r-lib/pkgdown,
diff --git a/NAMESPACE b/NAMESPACE
index 26f7df4d..df0c318d 100755
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -17,6 +17,7 @@ S3method(plot,PM_pta)
S3method(plot,PM_sim)
S3method(plot,PM_valid)
S3method(plot,PMvalid)
+S3method(plot,bd)
S3method(print,PM_compare)
S3method(print,PMerr)
S3method(print,PMnpc)
@@ -98,6 +99,8 @@ export(add_smooth)
export(additive)
export(all_is_numeric)
export(apps)
+export(bd)
+export(bd_post)
export(check_updates)
export(clear_build)
export(cli_ask)
@@ -183,6 +186,7 @@ importFrom(dplyr,group_by)
importFrom(dplyr,group_map)
importFrom(dplyr,if_else)
importFrom(dplyr,inner_join)
+importFrom(dplyr,left_join)
importFrom(dplyr,mutate)
importFrom(dplyr,n)
importFrom(dplyr,nest_by)
@@ -196,6 +200,7 @@ importFrom(dplyr,rowwise)
importFrom(dplyr,select)
importFrom(dplyr,slice)
importFrom(dplyr,slice_head)
+importFrom(dplyr,slice_max)
importFrom(dplyr,slice_tail)
importFrom(dplyr,starts_with)
importFrom(dplyr,summarize)
diff --git a/R/PM_bestdose.R b/R/PM_bestdose.R
new file mode 100644
index 00000000..b11b978a
--- /dev/null
+++ b/R/PM_bestdose.R
@@ -0,0 +1,1247 @@
+#' @title
+#' Object to contain BestDose optimization results
+#'
+#' @description
+#' `r lifecycle::badge("experimental")`
+#'
+#' This object is created after a successful BestDose optimization run.
+#' BestDose finds optimal dosing regimens to achieve target drug concentrations
+#' or AUC values using Bayesian optimization.
+#'
+#' @importFrom dplyr filter left_join group_by slice_max select ungroup bind_rows arrange pull group_map mutate
+#' @importFrom plotly plot_ly add_lines add_markers layout
+#' @importFrom purrr list_rbind
+#' @export
+bd <- R6::R6Class(
+ "bd",
+ public = list(
+ #' @field past PM_data object containing past patient data used in the optimization (if applicable)
+ past = NULL,
+ #' @field past_pred PM_sim object containing model predictions for the past data
+ past_pred = NULL,
+ #' @field future PM_data object containing future data used in the optimization
+ future = NULL,
+ #' @field future_pred PM_sim object containing model predictions for the future data
+ future_pred = NULL,
+ #' @field result List containing optimization results, including optimal doses, predictions, and objective function value
+ result = NULL,
+ #' @field posterior The `bd_post` object used to compute the posterior distribution (if applicable)
+ posterior = NULL,
+ #' @field prior_weight The prior weight (lambda) used in the optimization (if applicable)
+ prior_weight = NULL,
+ #' @field start Start specification provided by user (numeric hours or datetime string)
+ start = NULL,
+ #' @field start_offset Resolved numeric offset (hours) used internally; NULL means use target times as-is
+ start_offset = 0,
+ #' @field future_requested_doses Numeric vector of originally requested future doses (0 = optimize, non-zero = fixed)
+ future_requested_doses = NULL,
+ #' @field future_target_info List with target_type, target, and target_time from the original future specification
+ future_target_info = NULL,
+ #' @description
+ #' Initialize a `bd` object by running a one-shot BestDose optimization.
+ #' Creates the posterior internally, then optimizes. For reusing the posterior
+ #' across multiple optimizations, use `bd_post$new()` followed by `$optimize()`.
+ #' @param model PM_model object or path to compiled model
+ #' @param prior Prior parameter value distributions for the model, can be a [PM_result], [PM_final], or path to `theta.csv` file.
+ #' @param past_data PM_data object or path to CSV file with past patient data (optional). Omitting this will plan a new regimen for the patient
+ #' based on the population prior without any individualization other than inclusion of any covariates which are in the model.
+ #' @param future The future regimen to optimize. It can be one of three things:
+ #' * A list describing future dosing/target setup with elements:
+ #' - `dose` The amount to administer (required). Use `0` to optimize that dose to achieve the target; any other value is treated as fixed.
+ #' Can be a scalar or a numeric vector. If a vector is provided, each element corresponds to each dose event in order.
+ #' - `frequency` How often to administer dose(s) (default 24 hours). Can be a scalar or numeric vector; when vectorized,
+ #' element `i` is the interval after dose `i` and before dose `i+1`.
+ #' - `route` Either 0 for oral/bolus, or a value >0 indicating infusion duration. Can be a scalar or numeric vector;
+ #' when vectorized, element `i` applies to dose `i`.
+ #' - `number` The total number of doses to administer in the future (default 1).
+ #' If `number` is smaller than `max(length(dose), length(frequency), length(route))`, it is increased to that maximum and an informational message is shown.
+ #' If any of `dose`, `frequency`, or `route` are shorter than `number`, their last value is recycled to length `number`.
+ #' - `target` The target value to achieve after every dose (required).
+ #' - `target_type` The target type (default `"concentration"`). Must be one of `"concentration"`, `"auc"`, or `"time"`.
+ #' `"concentration"` targets a concentration at `target_time`; `"auc"` targets an AUC from dose time to `target_time`;
+ #' `"time"` uses `target_time` as a proportion between 0 and 1 of each dose interval.
+ #' - `target_time` For `target_type = "concentration"` or `"auc"`, the time after each dose at which to evaluate the target,
+ #' default 24 hours if unspecified. For `target_type = "time"`, this must be a single numeric proportion between 0 and 1,
+ #' and each target observation time is set to `round(frequency[i] * target_time)` hours after dose `i`.
+ #' - `covariates` Named list of covariates in the model with values to use. If there are no covariates in the model, this can be omitted.
+ #' * A `PM_data` object containing the future regimen and target observations. Dose events should be included with `evid=1` and target observations should be included with `evid=0` and `out` equal to the target value.
+ #' When `future` is supplied this way, concentration targets are assumed.
+ #' * A path to a CSV file containing the future regimen and target observations in the same format as described for the `PM_data` object above.
+ #' When `future` is supplied this way, concentration targets are assumed.
+ #' @param dose_range List with 'min' and 'max' elements defining the dose search range (default: 0 to 1000)
+ #' @param prior_weight Numeric between 0 and 1 indicating the balance between the prior and the posterior in the optimization (default: 0.5).
+ #' If 1, the model prior parameter
+ #' value distribution will be used as the sole basis for optimization. If 0, only the Bayesian posterior parameter distribution will be used,
+ #' which incorporates the patient history. Values between 0 and 1 will weight the contribution of the prior vs posterior parameter distributions in the optimization.
+ #' Choose values closer to 1 when the patient history is sparse or not believed to be very informative, and values closer to 0 when the patient history is rich and believed to be highly informative.
+ #' @param start Start time for the future regimen. Can be either:
+ #' * `NULL`: use the supplied future/target times as-is, without appending them after past data.
+ #' * Numeric hours (default `0`): relative to the last event date/time in past data, or if no past data,
+ #' relative to the next local clock hour.
+ #' * Date-time character string: interpreted as an absolute start date/time and converted to an internal hour offset.
+ #' Ambiguous formats (e.g. `01/02/26`) are resolved using the `date_format` Pmetrics option
+ #' (set via [setPMoptions]), which defaults to `%m/%d/%y` for US locales and `%d/%m/%y` otherwise.
+ #' Accepted formats include:
+ #' - `YYYY-mm-dd HH:MM[:SS]`
+ #' - `mm/dd/[YY]YY HH:MM[:SS]` or `mm-dd-[YY]YY HH:MM[:SS]`
+ #' - `dd/mm/[YY]YY HH:MM[:SS]` or `dd-mm-[YY]YY HH:MM[:SS]`
+ #' @param max_cycles Maximum number of optimization cycles for computing the posterior (default: 500)
+ #' @param settings List of additional settings for posterior computation (optional)
+ #' @param posterior `bd_post` object to use instead of computing a new one (optional).
+ #' When provided, `prior`, `model`, `past_data`, `max_cycles`, and `settings` are ignored, as they are already included in the `posterior` object.
+ #' @param quiet Logical indicating whether to suppress verbose simulation output (default: FALSE)
+ #' @return A `bd` object containing the optimization results and associated information
+
+ initialize = function(model = NULL,
+ prior = NULL,
+ past_data = NULL,
+ future = NULL,
+ dose_range = list(min = 0, max = 1000),
+ prior_weight = 0.5,
+ start = 0,
+ max_cycles = 500,
+ settings = NULL,
+ posterior = NULL,
+ quiet = FALSE,
+ simulate = FALSE) {
+ if (!is.null(future)) {
+ if (is.list(future)) {
+ future_spec <- private$.build_future_data(future)
+ future_data <- future_spec$data
+ target_type <- future_spec$target_type
+ self$future_target_info <- list(
+ target_type = tolower(as.character(future$target_type %||% "concentration")),
+ target = future$target,
+ target_time = future$target_time
+ )
+ } else if (inherits(future, "PM_data") || is.character(future)) {
+ future_data <- PM_data$new(future, quiet = TRUE)
+ target_type <- "concentration"
+ self$future_target_info <- list(target_type = "concentration", target = NA, target_time = NA)
+ } else {
+ cli::cli_abort("{.arg future} must be either a list describing the future setup, a PM_data object, or a path to a CSV file.")
+ }
+ } else {
+ cli::cli_abort("A {.arg future} plan must be supplied.")
+ }
+
+ # Capture originally requested doses before optimization may mutate the PM_data object
+ if (!is.null(future_data$standard_data)) {
+ orig_evid1 <- future_data$standard_data |> dplyr::filter(evid == 1)
+ self$future_requested_doses <- if (nrow(orig_evid1) > 0) orig_evid1$dose else numeric(0)
+ }
+
+ if (is.null(posterior)) {
+ if (is.null(prior) || is.null(model)) {
+ cli::cli_abort("Both {.arg prior} and {.arg model} are required when {.arg posterior} is not provided.")
+ }
+ posterior <- bd_post$new(
+ prior = prior,
+ model = model,
+ past_data = past_data,
+ max_cycles = max_cycles,
+ settings = settings,
+ quiet = quiet
+ )
+ }
+
+ start_info <- private$.resolve_start(start, posterior$past)
+
+ raw <- private$.optimize(posterior, future_data, dose_range, prior_weight, target_type, start_info$start_offset, quiet)
+
+ private$.set_result(future_data, raw, posterior, prior_weight)
+ self$past <- posterior$past
+ self$start <- start_info$start
+ self$start_offset <- start_info$start_offset
+
+ if (!is.null(self$past) && !is.null(self$start_offset) && self$start_offset > 0) {
+ private$.extend_past_to_start(self$start_offset)
+ }
+
+ if (simulate) {
+ cli::cli_alert("Simulating predictions for past and future data...")
+ private$.sim_past()
+ private$.sim_future()
+ }
+
+ invisible(self)
+ },
+
+ #' @description
+ #' Print summary of BestDose results
+ print = function() {
+ cat("BestDose Optimization Results\n")
+ cat("==============================\n\n")
+ dose_str <- paste(sprintf("%.2f", self$doses), collapse = ", ")
+ cat(sprintf("Optimal doses: [%s] mg\n", dose_str))
+ cat(sprintf("Objective function: %.10f\n", self$objf))
+ cat(sprintf("ln(Objective): %.4f\n", log(self$objf)))
+ cat(sprintf("Method: %s\n", self$method))
+ cat(sprintf("Status: %s\n", self$status))
+ if (!is.null(self$prior_weight)) {
+ cat(sprintf("Prior weight: %.2f\n", self$prior_weight))
+ }
+ if (!is.null(self$start)) {
+ if (is.character(self$start)) {
+ cat(sprintf("Start: %s\n", self$start))
+ } else {
+ cat(sprintf("Start offset (h): %.2f\n", self$start_offset))
+ }
+ }
+ cat(sprintf("\nNumber of predictions: %d\n", nrow(self$result$predictions)))
+ if (!is.null(self$result$auc_predictions)) {
+ cat(sprintf("Number of AUC predictions: %d\n", nrow(self$result$auc_predictions)))
+ }
+ invisible(self)
+ },
+ #' @description
+ #' Save results to RDS file
+ #' @param filename Path to save file. Default: "bestdose_result.rds"
+ save = function(filename = "bestdose_result.rds") {
+ saveRDS(self, filename)
+ cli::cli_alert_success("Results saved to {filename}")
+ invisible(self)
+ },
+ #' @description
+ #' Plot observed and predicted concentrations over time for both past and future data
+ #' @param ... Additional arguments passed to plot.bd function
+ plot = function(...) {
+ plot.bd(self, ...)
+ },
+ #' @description
+ #' Generate an HTML report for a BestDose result
+ #' @param ... Parameters passed to the internal BestDose report generator.
+ report = function(...) {
+ bd_report(self, ...)
+ }
+ ),
+ private = list(
+ .build_future_data = function(future) {
+ 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))
+ }
+ doses <- recycle_last(future$dose, n)
+ freqs <- recycle_last(future$frequency, n)
+ routes <- recycle_last(future$route, n)
+ target_type <- tolower(as.character(future$target_type[[1]]))
+
+ if (!target_type %in% c("concentration", "auc", "time")) {
+ cli::cli_abort("{.arg future$target_type} must be one of: concentration, auc, time")
+ }
+
+ if (target_type == "time") {
+ if (is.null(future$target_time) ||
+ length(future$target_time) != 1 ||
+ !is.numeric(future$target_time) ||
+ is.na(future$target_time) ||
+ future$target_time < 0 ||
+ future$target_time > 1) {
+ cli::cli_abort("For {.code future$target_type = 'time'}, {.arg future$target_time} must be a single numeric value between 0 and 1.")
+ }
+ target_offsets <- round(freqs * future$target_time)
+ rust_target_type <- "concentration"
+ } else {
+ future$target_time <- future$target_time %||% 24
+ if (length(future$target_time) != 1 ||
+ !is.numeric(future$target_time) ||
+ is.na(future$target_time) ||
+ future$target_time < 0) {
+ cli::cli_abort("{.arg future$target_time} must be a single non-negative numeric value.")
+ }
+ target_offsets <- rep(as.numeric(future$target_time), n)
+ rust_target_type <- if (target_type == "auc") "auc_from_zero" else "concentration"
+ }
+
+ # Cumulative dose times: dose 1 at t=0, dose i at sum of preceding frequencies
+ dose_times <- c(0, cumsum(head(freqs, n - 1)))
+
+ future_data <- PM_data$new()
+
+ # Add each dose event individually
+ for (i in seq_len(n)) {
+ dose_args <- list(
+ id = 1, time = dose_times[i], evid = 1,
+ dose = doses[i], out = NA,
+ dur = routes[i],
+ validate = FALSE
+ )
+ do.call(future_data$addEvent, dose_args)
+ }
+
+ # Add covariates if specified
+ if (!is.null(future$covariates) && length(future$covariates) > 0) {
+ cov_args <- c(list(id = 1, validate = FALSE), future$covariates)
+ do.call(future_data$addEvent, cov_args)
+ }
+
+ # Add observation events at target_time after each dose
+ for (i in seq_len(n)) {
+ obs_time <- dose_times[i] + target_offsets[i]
+ obs_args <- list(
+ id = 1, time = obs_time, evid = 0,
+ dose = NA, out = future$target,
+ validate = (i == n)
+ )
+ do.call(future_data$addEvent, obs_args)
+ }
+ list(data = future_data, target_type = rust_target_type)
+ },
+ .optimize = function(posterior, future_data, dose_range, prior_weight, target_type, start_offset, quiet = FALSE) {
+ if (is.null(posterior$handle)) {
+ cli::cli_abort(c(
+ "x" = "bd_post object is not properly initialized.",
+ "i" = "Create a new {.cls bd_post} with {.code bd_post$new()}."
+ ))
+ }
+
+ if (!target_type %in% c("concentration", "auc_from_zero")) {
+ cli::cli_abort("Internal {.arg target_type} must be one of: concentration, auc_from_zero")
+ }
+
+ if (prior_weight < 0 || prior_weight > 1) {
+ cli::cli_abort("{.arg prior_weight} must be between 0 and 1")
+ }
+
+ if (is.null(dose_range$min) || is.null(dose_range$max)) {
+ cli::cli_abort("{.arg dose_range} must have both 'min' and 'max' elements")
+ }
+
+ if (dose_range$min >= dose_range$max) {
+ cli::cli_abort("{.arg dose_range$min} must be less than {.arg dose_range$max}")
+ }
+
+ future_data_path <- bestdose_parse_data(future_data)
+
+ res <- bestdose_optimize(
+ posterior$handle,
+ future_data_path,
+ start_offset,
+ dose_range$min,
+ dose_range$max,
+ bias_weight = prior_weight,
+ target_type
+ )
+
+ if (is.character(res)) {
+ cli::cli_abort(res)
+ }
+
+ res
+ },
+ .set_result = function(future, result, posterior, prior_weight) {
+ if (!is.null(future)) {
+ self$future <- if (inherits(future, "PM_data")) future else PM_data$new(future, quiet = TRUE)
+ }
+ self$result <- result
+ self$posterior <- posterior
+ self$prior_weight <- prior_weight
+ invisible(self)
+ },
+ .resolve_start = function(start, past) {
+ if (is.null(start)) {
+ start_offset <- NULL
+ } else if (is.numeric(start) && length(start) == 1 && !is.na(start)) {
+ start_offset <- as.numeric(start)
+ } else if (is.character(start) && length(start) == 1 && nzchar(start)) {
+ start_datetime <- private$.parse_datetime(start)
+ if (is.na(start_datetime)) {
+ cli::cli_abort(c(
+ "x" = "Unable to parse {.arg start} datetime string.",
+ "i" = "Use an ISO-like value such as {.code '2026-02-01 14:00:00'}."
+ ))
+ }
+
+ has_past <- !is.null(past) && !is.null(past$data) && nrow(past$data) > 0
+ if (has_past) {
+ ref_datetime <- private$.reference_datetime(past)
+ start_offset <- as.numeric(difftime(start_datetime, ref_datetime, units = "hours"))
+ } else {
+ # Without past data, a datetime string defines the plotting origin itself.
+ # Keep optimization on the future timeline by using a zero internal offset.
+ start_offset <- 0
+ }
+ } else {
+ cli::cli_abort("{.arg start} must be NULL, a single numeric value, or a date-time character string.")
+ }
+
+ if (!is.null(start_offset) && (!is.finite(start_offset) || start_offset < 0)) {
+ cli::cli_abort("{.arg start} resolves to a negative or invalid offset relative to the reference time.")
+ }
+
+ list(start = start, start_offset = start_offset)
+ },
+ .reference_datetime = function(past) {
+ if (!is.null(past) && !is.null(past$data) && nrow(past$data) > 0 &&
+ all(c("date", "time") %in% names(past$data))) {
+ dt_raw <- paste(as.character(past$data$date), as.character(past$data$time))
+ dt_vals <- vapply(
+ dt_raw,
+ function(x) private$.parse_datetime(x),
+ FUN.VALUE = as.POSIXct(NA_real_, origin = "1970-01-01", tz = Sys.timezone())
+ )
+
+ if (!all(is.na(dt_vals))) {
+ if (!is.null(past$standard_data) && "time" %in% names(past$standard_data) &&
+ nrow(past$standard_data) == length(dt_vals)) {
+ idx <- which.max(past$standard_data$time)
+ if (length(idx) == 1 && !is.na(dt_vals[idx])) {
+ return(dt_vals[idx])
+ }
+ }
+ return(max(dt_vals, na.rm = TRUE))
+ }
+ }
+
+ now <- Sys.time()
+ ceiling_hour <- as.POSIXct(
+ format(now, "%Y-%m-%d %H:00:00"),
+ format = "%Y-%m-%d %H:%M:%S",
+ tz = Sys.timezone()
+ )
+ if (ceiling_hour <= now) ceiling_hour <- ceiling_hour + 3600
+ ceiling_hour
+ },
+ .parse_datetime = function(x) {
+ # Get user-preferred date format from Pmetrics options; fall back to locale
+ user_date_fmt <- getPMoptions("date_format", warn = FALSE, quiet = TRUE)
+ if (!is.character(user_date_fmt) || !nzchar(user_date_fmt)) {
+ user_date_fmt <- if (grepl("en_US", Sys.getlocale("LC_TIME"), fixed = TRUE)) "%m/%d/%y" else "%d/%m/%y"
+ }
+ user_4yr <- sub("%y", "%Y", user_date_fmt, fixed = TRUE)
+ formats <- unique(c(
+ paste(user_4yr, "%H:%M:%S"),
+ paste(user_4yr, "%H:%M"),
+ paste(user_date_fmt, "%H:%M:%S"),
+ paste(user_date_fmt, "%H:%M"),
+ "%Y-%m-%d %H:%M:%S",
+ "%Y-%m-%d %H:%M",
+ "%m/%d/%Y %H:%M:%S",
+ "%m/%d/%Y %H:%M",
+ "%m/%d/%y %H:%M:%S",
+ "%m/%d/%y %H:%M",
+ "%m-%d-%Y %H:%M",
+ "%m-%d-%Y %H:%M:%S",
+ "%m-%d-%y %H:%M",
+ "%m-%d-%y %H:%M:%S",
+ "%d/%m/%Y %H:%M:%S",
+ "%d/%m/%Y %H:%M",
+ "%d/%m/%y %H:%M:%S",
+ "%d/%m/%y %H:%M",
+ "%d-%m-%Y %H:%M:%S",
+ "%d-%m-%Y %H:%M",
+ "%d-%m-%y %H:%M:%S",
+ "%d-%m-%y %H:%M"
+ ))
+ for (fmt in formats) {
+ dt <- suppressWarnings(as.POSIXct(x, format = fmt, tz = Sys.timezone()))
+ if (!is.na(dt)) {
+ return(dt)
+ }
+ }
+ suppressWarnings(as.POSIXct(x, tz = Sys.timezone()))
+ },
+ .extend_past_to_start = function(start_offset) {
+ if (is.null(self$past) || is.null(self$past$standard_data) || nrow(self$past$standard_data) == 0) {
+ return(invisible(NULL))
+ }
+
+ past_std <- self$past$standard_data
+ last_rows <- past_std |>
+ dplyr::group_by(id) |>
+ dplyr::slice_max(time, n = 1, with_ties = FALSE) |>
+ dplyr::ungroup()
+
+ new_rows <- last_rows
+ new_rows$time <- new_rows$time + start_offset
+ if ("evid" %in% names(new_rows)) new_rows$evid <- 0
+ if ("dose" %in% names(new_rows)) new_rows$dose <- NA_real_
+ if ("dur" %in% names(new_rows)) new_rows$dur <- NA_real_
+ if ("addl" %in% names(new_rows)) new_rows$addl <- NA_real_
+ if ("ii" %in% names(new_rows)) new_rows$ii <- NA_real_
+ if ("out" %in% names(new_rows)) new_rows$out <- NA_real_
+
+ self$past$standard_data <- dplyr::bind_rows(past_std, new_rows) |>
+ dplyr::arrange(id, time)
+
+ invisible(self)
+ },
+ .sim_past = function() {
+ if (is.null(self$past)) {
+ return(invisible(NULL))
+ }
+ sim1 <- PM_sim$new(
+ poppar = as.data.frame(self$posterior$theta),
+ model = self$posterior$model_info$model,
+ data = self$past,
+ predInt = 1,
+ limits = NA,
+ quiet = TRUE
+ )
+ self$past_pred <- sim1
+ invisible(self)
+ },
+ .sim_future = function() {
+ if (is.null(self$future)) {
+ return(invisible(NULL))
+ }
+ # Replace placeholder doses with optimized values (use tail to skip past fixed doses)
+ future_sim_data <- self$future$standard_data
+ n_future_doses <- sum(future_sim_data$evid != 0)
+ future_sim_data$dose[future_sim_data$evid != 0] <- tail(self$result$doses, n_future_doses)
+ self$future$standard_data <- future_sim_data
+
+ sim2 <- PM_sim$new(
+ poppar = as.data.frame(self$posterior$theta),
+ model = self$posterior$model_info$model,
+ data = self$future,
+ predInt = 1,
+ limits = NA,
+ quiet = TRUE
+ )
+ self$future_pred <- sim2
+ invisible(self)
+ }
+ ),
+ active = list(
+ #' @field doses Optimal doses found by BestDose optimization
+ doses = function() {
+ self$result$doses
+ },
+ #' @field objf Objective function value at the optimal doses
+ objf = function() {
+ self$result$objf
+ },
+ #' @field method Optimization method used
+ method = function() {
+ self$result$method
+ },
+ #' @field status Status of the optimization run
+ status = function() {
+ self$result$status
+ }
+ )
+)
+
+#' @title
+#' Compute a reusable BestDose posterior
+#'
+#' @description
+#' `r lifecycle::badge("experimental")`
+#'
+#' Use `bd_post` to compute the Bayesian posterior once from
+#' prior population data and patient history, then call `$optimize()` multiple
+#' times with different targets, dose ranges, or bias weights.
+#'
+#' @export
+bd_post <- R6::R6Class(
+ "bd_post",
+ public = list(
+ #' @field handle Memory pointer to the computed posterior (opaque to users)
+ handle = NULL,
+ #' @field theta Matrix of support points in the posterior distribution
+ theta = NULL,
+ #' @field theta_dim Dimensions of the theta matrix
+ theta_dim = NULL,
+ #' @field param_names Names of the parameters in the posterior
+ param_names = NULL,
+ #' @field posterior_weights Weights of the posterior support points
+ posterior_weights = NULL,
+ #' @field population_weights Weights of the population support points
+ population_weights = NULL,
+ #' @field past PM_data object containing past patient data used in the posterior computation (if applicable)
+ past = NULL,
+ #' @field model_info Information about the model used
+ model_info = NULL,
+ #' @field settings Settings used for the posterior computation
+ settings = NULL,
+ #' @description
+ #' Initialize the `bd_post` object by computing the posterior distribution from the given prior, model, and past data
+ #' @param prior Prior information for the model, can be a PM_result, PM_final, or path to theta.csv
+ #' @param model PM_model object or path to compiled model
+ #' @param past_data PM_data object or path to CSV file with past patient data (optional)
+ #' @param max_cycles Maximum number of optimization cycles for computing the posterior (default: 500)
+ #' @param settings List of additional settings for posterior computation (optional)
+ #' @param quiet Logical indicating whether to suppress verbose simulation output (default: FALSE)
+ #' @return A `bd_post` object containing the computed posterior distribution and associated information
+
+ initialize = function(prior,
+ model,
+ past_data = NULL,
+ max_cycles = 500,
+ settings = NULL,
+ quiet = FALSE) {
+ prior_path <- bestdose_parse_prior(prior)
+ model_info <- bestdose_parse_model(model)
+ past_data_path <- if (!is.null(past_data)) bestdose_parse_data(past_data) else NULL
+
+ if (is.null(settings)) {
+ model_for_settings <- if (!is.null(model_info$model)) model_info$model else model
+ settings <- bestdose_default_settings(prior, model_for_settings, max_cycles = max_cycles)
+ }
+
+ prep <- bestdose_prepare(
+ model_path = model_info$path,
+ prior_path = prior_path,
+ past_data_path = past_data_path,
+ params = settings,
+ kind = model_info$kind
+ )
+
+ if (is.character(prep)) {
+ cli::cli_abort(prep)
+ }
+
+ dim <- as.integer(prep$theta_dim)
+ theta_matrix <- matrix(prep$theta_values, nrow = dim[1], ncol = dim[2])
+ colnames(theta_matrix) <- prep$param_names
+
+ self$handle <- prep$handle
+ self$theta <- theta_matrix
+ self$theta_dim <- dim
+ self$param_names <- prep$param_names
+ self$posterior_weights <- prep$posterior_weights
+ self$population_weights <- prep$population_weights
+ self$past <- if (!is.null(past_data)) PM_data$new(past_data, quiet = TRUE) else NULL
+ self$model_info <- model_info
+ self$settings <- settings
+ private$.quiet <- quiet
+
+ cli::cli_alert_success("BestDose posterior computed with {dim[1]} support points")
+ },
+ #' @description
+ #' Run optimization and return a `bd` result object
+ #' @param target Future regimen/target specification. Accepts the same forms as the `future` argument to `bd$new()`:
+ #' a list with elements such as `dose`, `frequency`, `route`, `number`, `target`, `target_type`, and `target_time`;
+ #' a `PM_data` object; or a path to a CSV file.
+ #' @param dose_range List with 'min' and 'max' elements defining the dose search range (default: 0 to 1000)
+ #' @param prior_weight Numeric between 0 and 1 indicating the weight of bias in the optimization (default: 0.5)
+ #' @param start Start time for future regimen. Can be either `NULL`, numeric hours (default: `0`), or a date-time
+ #' character string. Ambiguous formats are resolved using the `date_format` Pmetrics option
+ #' (see [setPMoptions]). Accepted string formats:
+ #' - `YYYY-mm-dd HH:MM[:SS]`
+ #' - `mm/dd/[YY]YY HH:MM[:SS]` or `mm-dd-[YY]YY HH:MM[:SS]`
+ #' - `dd/mm/[YY]YY HH:MM[:SS]` or `dd-mm-[YY]YY HH:MM[:SS]`
+ #' @param quiet Logical indicating whether to suppress verbose simulation output. If NULL, uses the quiet setting from posterior computation (default: NULL)
+ #' @return A `bd` object containing the optimization results
+ optimize = function(target,
+ dose_range = list(min = 0, max = 1000),
+ prior_weight = 0.5,
+ start = 0,
+ quiet = NULL) {
+ if (is.null(quiet)) quiet <- private$.quiet
+ bd$new(
+ future = target,
+ dose_range = dose_range,
+ prior_weight = prior_weight,
+ start = start,
+ posterior = self,
+ quiet = quiet
+ )
+ }
+ ),
+ private = list(
+ .quiet = FALSE,
+ finalize = function() {
+ self$handle <- NULL
+ }
+ )
+)
+
+#' @export
+bd$load <- function(filename = "bestdose_result.rds") {
+ if (!file.exists(filename)) {
+ cli::cli_abort("File not found: {filename}")
+ }
+ readRDS(filename)
+}
+
+
+############ HELPER FUNCTIONS ############
+
+bestdose_parse_prior <- function(prior) {
+ if (inherits(prior, "PM_result")) {
+ theta_path <- file.path(prior$rundir, "outputs", "theta.csv")
+ if (!file.exists(theta_path)) {
+ cli::cli_abort("theta.csv not found in PM_result outputs")
+ }
+ theta_path
+ } else if (inherits(prior, "PM_final")) {
+ temp_path <- tempfile(fileext = ".csv")
+ bestdose_write_prior_csv(prior, temp_path)
+ temp_path
+ } else if (is.character(prior)) {
+ if (!file.exists(prior)) {
+ cli::cli_abort("Prior file not found: {prior}")
+ }
+ prior
+ } else {
+ cli::cli_abort("prior must be PM_result, PM_final, or path to theta.csv")
+ }
+}
+
+bestdose_write_prior_csv <- function(prior, path) {
+ df <- as.data.frame(prior$popPoints)
+ df$prob <- prior$popProb
+ write.csv(df, path, row.names = FALSE, quote = FALSE)
+}
+
+bestdose_parse_model <- function(model) {
+ if (inherits(model, "PM_model")) {
+ compiled_path <- model$binary_path
+ if (is.null(compiled_path) || !file.exists(compiled_path)) {
+ cli::cli_abort("Model must be compiled first. Use model$compile()")
+ }
+
+ kind <- if (!is.null(model$model_list$analytical) && model$model_list$analytical) {
+ "analytical"
+ } else {
+ "ode"
+ }
+
+ list(path = compiled_path, kind = kind, model = model)
+ } else if (is.character(model)) {
+ if (!file.exists(model)) {
+ cli::cli_abort("Model file not found: {model}")
+ }
+ kind <- if (grepl("analytical", model, ignore.case = TRUE)) {
+ "analytical"
+ } else {
+ "ode"
+ }
+ list(path = model, kind = kind, model = NULL)
+ } else {
+ cli::cli_abort("model must be PM_model or path to compiled model")
+ }
+}
+
+bestdose_parse_data <- function(data) {
+ if (inherits(data, "PM_data")) {
+ temp_path <- tempfile(fileext = ".csv")
+ write.csv(data$standard_data, temp_path, row.names = FALSE, quote = FALSE)
+ temp_path
+ } else if (is.character(data)) {
+ if (!file.exists(data)) {
+ cli::cli_abort("Data file not found: {data}")
+ }
+ data
+ } else {
+ cli::cli_abort("data must be PM_data or path to CSV file")
+ }
+}
+
+bestdose_default_settings <- function(prior, model, max_cycles = 500) {
+ param_ranges <- lapply(model$model_list$pri, function(x) {
+ c(x$min, x$max)
+ })
+ names(param_ranges) <- tolower(names(param_ranges))
+
+ list(
+ algorithm = "NPAG",
+ ranges = param_ranges,
+ error_models = lapply(model$model_list$err, function(x) x$flatten()),
+ max_cycles = max_cycles,
+ points = 2028,
+ seed = 22,
+ prior = "prior.csv",
+ idelta = 0.25,
+ tad = 0.0
+ )
+}
+
+
+#' Plot BestDose predictions
+#'
+#' @description
+#' Plot observed and predicted concentrations over time for both past and future data,
+#' with options to apply a multiplier to the concentrations, and customize the plot appearance. The top 5 posterior support
+#' points are highlighted with distinct colors.
+#'
+#' @param x A `bd` object containing the best dose predictions.
+#' @param mult Numeric multiplier to apply to the concentrations (default: 1, no scaling).
+#' @param outeq Numeric value of outeq to filter observations for plotting (default: 1).
+#' @param quiet Logical indicating whether to suppress messages (default: FALSE).
+#' @param legend Logical indicating whether to display a legend (default: TRUE).
+#' @param log Logical indicating whether to use a logarithmic scale for the y-axis (default: FALSE).
+#' @param grid Logical indicating whether to display a grid (default: FALSE).
+#' @param future_region Logical indicating whether to display the divider line and shaded region for future data (default: FALSE).
+#' @param xlab Label for the x-axis.
+#' @param ylab Label for the y-axis.
+#' @param title Title of the plot.
+#' @param xlim Limits for the x-axis.
+#' @param ylim Limits for the y-axis.
+#' @param print Logical indicating whether to print the plot (default: TRUE).
+#' @param ... Additional arguments passed to the plotting function.
+#'
+#' @return Invisibly returns a plotly object.
+#'
+#' @export
+#' @method plot bd
+plot.bd <- function(x, mult = 1, outeq = 1,
+ quiet = FALSE, legend = TRUE, log = FALSE,
+ grid = FALSE, future_region = FALSE,
+ xlab = "Time", ylab = "Concentration",
+ title = NULL, xlim = NULL, ylim = NULL, print = TRUE, ...) {
+ if (!inherits(x, "bd")) {
+ stop("Object must be of class 'bd'")
+ }
+
+ # 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.")
+ }
+ }
+
+ # Check if we have future predictions (past is optional)
+ if (is.null(x$future_pred)) {
+ cli::cli_abort(c(
+ "x" = "Future predictions are required for plotting.",
+ "i" = "BestDose object must have {.code future_pred}."
+ ))
+ }
+
+ # Get starting date/time for x-axis datetime conversion
+ start_datetime <- bd_plot_start_datetime(x)
+
+ # Compute a single, consistent future_time_offset (in hours) for plotting.
+ # This mirrors the Rust optimization logic:
+ # with past: future starts at max past time (past is optionally extended to start)
+ # without past: future starts at start_offset hours from plot origin
+ future_time_offset <- bd_future_time_offset(x)
+ has_past <- !is.null(x$past) && !is.null(x$past$standard_data) && nrow(x$past$standard_data) > 0
+
+ if (!quiet) {
+ cli::cli_alert_info("start = {deparse(x$start)}, start_offset = {x$start_offset %||% 0}, future_time_offset = {future_time_offset} h")
+ }
+
+ # Extract observation data from both simulations
+ past_obs <- NULL
+ if (!is.null(x$past_pred)) {
+ past_obs <- x$past_pred$data$obs |>
+ dplyr::filter(outeq == !!outeq)
+ }
+
+ future_obs <- x$future_pred$data$obs |>
+ dplyr::filter(outeq == !!outeq)
+
+
+ # Adjust future simulation prediction times using the unified offset
+ future_obs_adjusted <- future_obs |>
+ dplyr::mutate(time = time + future_time_offset)
+
+ # Combine past and future observations
+ if (!is.null(past_obs) && nrow(past_obs) > 0) {
+ combined_obs <- dplyr::bind_rows(
+ past_obs |> dplyr::mutate(source = "past"),
+ future_obs_adjusted |> dplyr::mutate(source = "future")
+ ) |>
+ dplyr::arrange(id, nsim, time) |>
+ dplyr::mutate(datetime = start_datetime + time * 3600)
+ } else {
+ combined_obs <- future_obs_adjusted |>
+ dplyr::mutate(source = "future") |>
+ dplyr::arrange(id, nsim, time) |>
+ dplyr::mutate(datetime = start_datetime + time * 3600)
+ }
+
+ # Apply multiplier
+ combined_obs$out <- combined_obs$out * mult
+
+ # Handle log scale
+ if (log) {
+ if (any(combined_obs$out <= 0, na.rm = TRUE)) {
+ if (!quiet) {
+ cat("Values <= 0 omitted from log plot.\n")
+ }
+ combined_obs$out[combined_obs$out <= 0] <- NA
+ }
+ }
+
+ # Remove NA values (but keep intentional break rows for gap rendering)
+ combined_obs <- combined_obs |> dplyr::filter(!is.na(out) | source == "break")
+
+ # Extract actual observations from past and future data
+ past_data_obs <- NULL
+ future_data_obs <- NULL
+ past_doses <- NULL
+ future_doses <- NULL
+
+ if (!is.null(x$past)) {
+ past_data_obs <- x$past$standard_data |>
+ dplyr::filter(evid == 0, outeq == !!outeq) |>
+ dplyr::select(id, time, out)
+
+ past_data_obs$out <- past_data_obs$out * mult
+
+ if (log && any(past_data_obs$out <= 0, na.rm = TRUE)) {
+ past_data_obs$out[past_data_obs$out <= 0] <- NA
+ }
+
+ past_data_obs <- past_data_obs |> dplyr::filter(!is.na(out))
+
+ past_doses <- x$past$standard_data |>
+ dplyr::filter(evid == 1) |>
+ dplyr::select(id, time, dur, dose)
+ }
+
+ if (!is.null(x$future)) {
+ future_data_obs <- x$future$standard_data |>
+ dplyr::filter(evid == 0, outeq == !!outeq) |>
+ dplyr::select(id, time, out)
+
+ # Use the same unified future_time_offset for consistency
+ future_data_obs <- future_data_obs |>
+ dplyr::mutate(time = time + future_time_offset)
+
+ future_data_obs$out <- future_data_obs$out * mult
+
+ if (log && any(future_data_obs$out <= 0, na.rm = TRUE)) {
+ future_data_obs$out[future_data_obs$out <= 0] <- NA
+ }
+
+ future_data_obs <- future_data_obs |> dplyr::filter(!is.na(out))
+
+ future_doses <- bd_future_dose_events(x)
+ if (!is.null(future_doses) && nrow(future_doses) > 0) {
+ future_doses <- future_doses |>
+ dplyr::select(id, time, dur, dose)
+ }
+ }
+
+ # Initialize plot
+ p <- plotly::plot_ly()
+
+ # Get posterior weights and robustly map them to nsim values
+ posterior_weights <- as.numeric(x$posterior$posterior_weights)
+ top5_colors <- c("#FF0000", "#0000FF", "#008000", "#800080", "#FFA500")
+
+ unique_nsim <- sort(unique(combined_obs$nsim))
+ nsim_keys <- as.character(unique_nsim)
+
+ sim_index_by_nsim <- setNames(seq_along(unique_nsim), nsim_keys)
+
+ weight_by_nsim <- setNames(rep(NA_real_, length(nsim_keys)), nsim_keys)
+ for (k in nsim_keys) {
+ idx <- sim_index_by_nsim[[k]]
+ if (!is.na(idx) && idx >= 1 && idx <= length(posterior_weights)) {
+ weight_by_nsim[[k]] <- posterior_weights[idx]
+ }
+ }
+
+ # Identify top 5 nsim keys by weight (largest to smallest)
+ ranked_keys <- names(sort(weight_by_nsim, decreasing = TRUE, na.last = NA))
+ top_keys <- ranked_keys[seq_len(min(5, length(ranked_keys)))]
+
+ color_by_nsim <- setNames(rep("#cccccc", length(nsim_keys)), nsim_keys)
+ legend_name_by_nsim <- setNames(rep(NA_character_, length(nsim_keys)), nsim_keys)
+ for (i in seq_along(top_keys)) {
+ k <- top_keys[i]
+ color_by_nsim[[k]] <- top5_colors[i]
+ legend_name_by_nsim[[k]] <- sprintf("%.4e", weight_by_nsim[[k]])
+ }
+
+ legend_added <- setNames(rep(FALSE, length(nsim_keys)), nsim_keys)
+ id_values <- unique(combined_obs$id)
+ non_top_keys <- setdiff(nsim_keys, top_keys)
+
+ # 1) Plot non-top lines first (light grey, no legend)
+ for (nsim_key in non_top_keys) {
+ weight_val <- weight_by_nsim[[nsim_key]]
+ for (id_val in id_values) {
+ this_sim <- combined_obs |>
+ dplyr::filter(id == id_val, as.character(nsim) == nsim_key) |>
+ dplyr::arrange(datetime)
+
+ if (nrow(this_sim) > 0) {
+ p <- p |>
+ plotly::add_lines(
+ x = ~datetime, y = ~out, data = this_sim,
+ line = list(color = "#cccccc", width = 0.8),
+ opacity = 0.6,
+ name = NULL,
+ legendgroup = NULL,
+ hovertemplate = sprintf(
+ "Time: %%{x}
Out: %%{y}
Weight: %s",
+ ifelse(is.na(weight_val), "NA", sprintf("%.4e", weight_val))
+ ),
+ showlegend = FALSE
+ )
+ }
+ }
+ }
+
+ # 2) Plot top lines second (colored, legend in descending weight order)
+ for (nsim_key in top_keys) {
+ line_color <- color_by_nsim[[nsim_key]]
+ legend_name <- legend_name_by_nsim[[nsim_key]]
+ weight_val <- weight_by_nsim[[nsim_key]]
+
+ for (id_val in id_values) {
+ this_sim <- combined_obs |>
+ dplyr::filter(id == id_val, as.character(nsim) == nsim_key) |>
+ dplyr::arrange(datetime)
+
+ if (nrow(this_sim) > 0) {
+ show_in_legend <- !legend_added[[nsim_key]] && legend
+ if (show_in_legend) {
+ legend_added[[nsim_key]] <- TRUE
+ }
+ legend_rank <- match(nsim_key, top_keys)
+
+ p <- p |>
+ plotly::add_lines(
+ x = ~datetime, y = ~out, data = this_sim,
+ opacity = 0.6,
+ line = list(color = line_color, width = 2),
+ name = legend_name,
+ legendgroup = legend_name,
+ legendrank = legend_rank,
+ hovertemplate = sprintf(
+ "Time: %%{x}
Out: %%{y}
Weight: %s",
+ ifelse(is.na(weight_val), "NA", sprintf("%.4e", weight_val))
+ ),
+ showlegend = show_in_legend
+ )
+ }
+ }
+ }
+
+ # Add past data observations as markers
+ if (!is.null(past_data_obs) && nrow(past_data_obs) > 0) {
+ past_data_obs$datetime <- start_datetime + past_data_obs$time * 3600
+ p <- p |>
+ plotly::add_markers(
+ x = ~datetime, y = ~out, data = past_data_obs,
+ marker = list(color = "black", size = 8, symbol = "circle"),
+ name = "Past Observations",
+ hovertemplate = "Time: %{x}
Out: %{y}",
+ showlegend = FALSE
+ )
+ }
+
+ # Add future data observations as markers
+ if (!is.null(future_data_obs) && nrow(future_data_obs) > 0) {
+ future_data_obs$datetime <- start_datetime + future_data_obs$time * 3600
+ p <- p |>
+ plotly::add_markers(
+ x = ~datetime, y = ~out, data = future_data_obs,
+ marker = list(color = "red", size = 8, symbol = "circle"),
+ name = "Future Observations",
+ hovertemplate = "Time: %{x}
Out: %{y}",
+ showlegend = FALSE
+ )
+ }
+
+ # Add dose indicators as markers/segments directly on the plot
+ # Past doses (black triangles)
+ if (!is.null(past_doses) && nrow(past_doses) > 0) {
+ y_min <- min(combined_obs$out, na.rm = TRUE)
+
+ for (i in 1:nrow(past_doses)) {
+ dose_time <- past_doses$time[i]
+ dose_datetime <- start_datetime + dose_time * 3600
+ dose_dur <- past_doses$dur[i]
+ dose_amt <- past_doses$dose[i]
+
+ dur_text <- if (is.na(dose_dur) || dose_dur == 0) "Bolus" else sprintf("Duration: %.1f", dose_dur)
+ p <- p |>
+ plotly::add_markers(
+ x = dose_datetime,
+ y = y_min,
+ marker = list(
+ symbol = "triangle-up-open",
+ size = 10,
+ color = "black",
+ line = list(color = "black", width = 1)
+ ),
+ customdata = list(list(dose_amt, dur_text)),
+ hovertemplate = "Time: %{x}
Dose: %{customdata[0]}
%{customdata[1]}",
+ showlegend = FALSE
+ )
+
+ # If infusion, add marker at end time too
+ if (!is.na(dose_dur) && dose_dur > 0) {
+ p <- p |>
+ plotly::add_markers(
+ x = dose_datetime + dose_dur * 3600,
+ y = y_min,
+ marker = list(
+ symbol = "triangle-down-open",
+ size = 10,
+ color = "black",
+ line = list(color = "black", width = 1)
+ ),
+ text = "End infusion",
+ hovertemplate = "Time: %{x}
%{text}",
+ showlegend = FALSE
+ )
+ }
+ }
+ }
+
+ # Future doses (red triangles)
+ if (!is.null(future_doses) && nrow(future_doses) > 0) {
+ y_min <- min(combined_obs$out, na.rm = TRUE)
+
+ for (i in 1:nrow(future_doses)) {
+ dose_time <- future_doses$time[i]
+ dose_datetime <- start_datetime + dose_time * 3600
+ dose_dur <- future_doses$dur[i]
+ dose_amt <- future_doses$dose[i]
+
+ dur_text <- if (is.na(dose_dur) || dose_dur == 0) "Bolus" else sprintf("Duration: %.1f", dose_dur)
+ p <- p |>
+ plotly::add_markers(
+ x = dose_datetime,
+ y = y_min,
+ marker = list(
+ symbol = "triangle-up-open",
+ size = 10,
+ color = "red",
+ line = list(color = "red", width = 1)
+ ),
+ customdata = list(list(dose_amt, dur_text)),
+ hovertemplate = "Time: %{x}
Dose: %{customdata[0]}
%{customdata[1]}",
+ showlegend = FALSE
+ )
+
+ if (!is.na(dose_dur) && dose_dur > 0) {
+ p <- p |>
+ plotly::add_markers(
+ x = dose_datetime + dose_dur * 3600,
+ y = y_min,
+ marker = list(
+ symbol = "triangle-down-open",
+ size = 10,
+ color = "red",
+ line = list(color = "red", width = 1)
+ ),
+ text = "End infusion",
+ hovertemplate = "Time: %{x}
%{text}",
+ showlegend = FALSE
+ )
+ }
+ }
+ }
+
+ # Add vertical dashed line and background shading for the future region
+ future_boundary_datetime <- NULL
+ if (future_region && has_past && any(combined_obs$source == "future")) {
+ # Use future_time_offset directly for exact alignment with doses/obs
+ future_boundary_datetime <- start_datetime + future_time_offset * 3600
+ x_max <- max(combined_obs$datetime)
+ y_range <- range(combined_obs$out, na.rm = TRUE)
+
+ # Add light background rectangle for the future region
+ p <- p |>
+ plotly::layout(
+ shapes = list(
+ list(
+ type = "rect",
+ x0 = future_boundary_datetime,
+ x1 = x_max,
+ y0 = 0,
+ y1 = 1,
+ xref = "x",
+ yref = "paper",
+ fillcolor = "rgba(255, 200, 200, 0.15)",
+ line = list(width = 0),
+ layer = "below"
+ ),
+ list(
+ type = "line",
+ x0 = future_boundary_datetime,
+ x1 = future_boundary_datetime,
+ y0 = 0,
+ y1 = 1,
+ xref = "x",
+ yref = "paper",
+ line = list(
+ color = "rgba(180, 0, 0, 0.6)",
+ width = 2,
+ dash = "dash"
+ )
+ )
+ ),
+ annotations = list(
+ list(
+ x = future_boundary_datetime,
+ y = 1.02,
+ xref = "x",
+ yref = "paper",
+ text = "Future starts",
+ showarrow = FALSE,
+ font = list(color = "rgba(180, 0, 0, 0.8)", size = 11)
+ )
+ )
+ )
+ }
+
+ # Set axis labels and layout
+ xlab_val <- if (missing(xlab)) "Time" else xlab
+ ylab_val <- if (missing(ylab)) "Concentration" else ylab
+ title_val <- if (missing(title)) "BestDose Predictions: Past and Future" else title
+
+ p$x$layout$title <- title_val
+ p$x$layout$xaxis <- list(
+ title = xlab_val,
+ zeroline = FALSE,
+ showgrid = grid
+ )
+ p$x$layout$yaxis <- list(
+ title = ylab_val,
+ type = if (log) "log" else "linear",
+ zeroline = FALSE,
+ showgrid = grid
+ )
+ p$x$layout$hovermode <- "closest"
+ p$x$layout$showlegend <- legend
+ p$x$layout$legend <- list(traceorder = "normal")
+
+ if (!missing(xlim)) {
+ p$x$layout$xaxis$range <- xlim
+ }
+ if (!missing(ylim)) {
+ p$x$layout$yaxis$range <- ylim
+ }
+
+ if (print) print(p)
+ return(invisible(p))
+}
diff --git a/R/PM_bestdose_report.R b/R/PM_bestdose_report.R
new file mode 100644
index 00000000..ff319ff2
--- /dev/null
+++ b/R/PM_bestdose_report.R
@@ -0,0 +1,1286 @@
+bd_parse_datetime <- function(x) {
+ user_date_fmt <- getPMoptions("date_format", warn = FALSE, quiet = TRUE)
+ if (!is.character(user_date_fmt) || length(user_date_fmt) != 1 || !nzchar(user_date_fmt)) {
+ user_date_fmt <- if (grepl("en_US", Sys.getlocale("LC_TIME"), fixed = TRUE)) "%m/%d/%y" else "%d/%m/%y"
+ }
+
+ user_4yr <- sub("%y", "%Y", user_date_fmt, fixed = TRUE)
+ formats <- unique(c(
+ paste(user_4yr, "%H:%M:%S"),
+ paste(user_4yr, "%H:%M"),
+ paste(user_date_fmt, "%H:%M:%S"),
+ paste(user_date_fmt, "%H:%M"),
+ "%Y-%m-%d %H:%M:%S",
+ "%Y-%m-%d %H:%M",
+ "%m/%d/%Y %H:%M:%S",
+ "%m/%d/%Y %H:%M",
+ "%m/%d/%y %H:%M:%S",
+ "%m/%d/%y %H:%M",
+ "%m-%d-%Y %H:%M:%S",
+ "%m-%d-%Y %H:%M",
+ "%m-%d-%y %H:%M:%S",
+ "%m-%d-%y %H:%M",
+ "%d/%m/%Y %H:%M:%S",
+ "%d/%m/%Y %H:%M",
+ "%d/%m/%y %H:%M:%S",
+ "%d/%m/%y %H:%M",
+ "%d-%m-%Y %H:%M:%S",
+ "%d-%m-%Y %H:%M",
+ "%d-%m-%y %H:%M:%S",
+ "%d-%m-%y %H:%M"
+ ))
+
+ for (fmt in formats) {
+ dt <- suppressWarnings(as.POSIXct(x, format = fmt, tz = Sys.timezone()))
+ if (!is.na(dt)) {
+ return(dt)
+ }
+ }
+
+ suppressWarnings(as.POSIXct(x, tz = Sys.timezone()))
+}
+
+bd_plot_start_datetime <- function(x) {
+ if (!is.null(x$past) && !is.null(x$past$data) && nrow(x$past$data) > 0 &&
+ all(c("date", "time") %in% names(x$past$data))) {
+ dt <- bd_parse_datetime(paste(as.character(x$past$data$date[1]), as.character(x$past$data$time[1])))
+ if (!is.na(dt)) {
+ return(dt)
+ }
+ }
+
+ if (is.character(x$start) && length(x$start) == 1 && nzchar(x$start)) {
+ dt <- bd_parse_datetime(x$start)
+ if (!is.na(dt)) {
+ return(dt)
+ }
+ }
+
+ now <- Sys.time()
+ ceiling_hour <- as.POSIXct(
+ format(now, "%Y-%m-%d %H:00:00"),
+ format = "%Y-%m-%d %H:%M:%S",
+ tz = Sys.timezone()
+ )
+ if (ceiling_hour <= now) ceiling_hour <- ceiling_hour + 3600
+ ceiling_hour
+}
+
+bd_future_time_offset <- function(x) {
+ has_past <- !is.null(x$past) && !is.null(x$past$standard_data) && nrow(x$past$standard_data) > 0
+ if (has_past) {
+ if (is.null(x$start_offset)) {
+ 0
+ } else {
+ max(x$past$standard_data$time, na.rm = TRUE)
+ }
+ } else {
+ x$start_offset %||% 0
+ }
+}
+
+bd_expand_dose_events <- function(data) {
+ if (is.null(data) || nrow(data) == 0) {
+ return(data.frame())
+ }
+
+ rows <- lapply(seq_len(nrow(data)), function(i) {
+ row <- data[i, , drop = FALSE]
+ addl <- if ("addl" %in% names(row) && !is.na(row$addl[1])) as.integer(row$addl[1]) else 0L
+ n_doses <- addl + 1L
+ ii <- if ("ii" %in% names(row) && !is.na(row$ii[1])) row$ii[1] else 0
+ times <- row$time[1] + seq.int(0, n_doses - 1L) * ii
+ expanded <- row[rep(1, n_doses), , drop = FALSE]
+ expanded$time <- times
+ expanded$addl <- NA
+ expanded$ii <- NA
+ expanded
+ })
+
+ do.call(rbind, rows)
+}
+
+bd_future_dose_events <- function(x) {
+ if (is.null(x$future) || is.null(x$future$standard_data)) {
+ return(NULL)
+ }
+
+ dose_rows <- x$future$standard_data |>
+ dplyr::filter(evid == 1)
+
+ if (nrow(dose_rows) == 0) {
+ return(NULL)
+ }
+
+ dose_rows <- bd_expand_dose_events(dose_rows)
+ if (nrow(dose_rows) == 0) {
+ return(NULL)
+ }
+
+ result_doses <- x$result$doses %||% numeric()
+ if (length(result_doses) == nrow(dose_rows)) {
+ dose_rows$dose <- as.numeric(result_doses)
+ } else if (length(result_doses) > 0 && all(is.na(dose_rows$dose) | dose_rows$dose == 0)) {
+ dose_rows$dose[seq_len(min(length(result_doses), nrow(dose_rows)))] <- as.numeric(result_doses[seq_len(min(length(result_doses), nrow(dose_rows)))])
+ }
+
+ dose_rows |>
+ dplyr::mutate(time = time + bd_future_time_offset(x))
+}
+
+bd_weight_table <- function(x, nsim_values) {
+ nsim_values <- sort(unique(nsim_values))
+ tibble::tibble(
+ nsim = nsim_values,
+ weight = as.numeric(x$posterior$posterior_weights)[seq_along(nsim_values)]
+ )
+}
+
+bd_weighted_pred_summary <- function(x, sim_data) {
+ if (is.null(sim_data) || nrow(sim_data) == 0) {
+ return(NULL)
+ }
+
+ sim_data <- sim_data |>
+ dplyr::mutate(id = as.character(id))
+
+ weights <- bd_weight_table(x, sim_data$nsim)
+
+ sim_data |>
+ dplyr::left_join(weights, by = dplyr::join_by(nsim)) |>
+ dplyr::group_by(id, time, outeq) |>
+ dplyr::summarise(
+ pred = stats::weighted.mean(out, w = weight, na.rm = TRUE),
+ .groups = "drop"
+ )
+}
+
+bd_fit_data <- function(x, source = c("past", "future"), outeq = 1) {
+ source <- match.arg(source)
+
+ if (source == "past") {
+ if (is.null(x$past) || is.null(x$past$standard_data) || is.null(x$past_pred)) {
+ return(NULL)
+ }
+
+ actual <- x$past$standard_data |>
+ dplyr::filter(evid == 0, outeq == !!outeq) |>
+ dplyr::transmute(id = as.character(id), time, outeq, obs = out, c0, c1, c2, c3)
+
+ sim <- x$past_pred$data$obs |>
+ dplyr::filter(outeq == !!outeq)
+
+ label <- "Measured"
+ } else {
+ if (is.null(x$future) || is.null(x$future$standard_data) || is.null(x$future_pred)) {
+ return(NULL)
+ }
+
+ offset <- bd_future_time_offset(x)
+ actual <- x$future$standard_data |>
+ dplyr::filter(evid == 0, outeq == !!outeq) |>
+ dplyr::transmute(id = as.character(id), time = time + offset, outeq, obs = out, c0, c1, c2, c3)
+
+ sim <- x$future_pred$data$obs |>
+ dplyr::filter(outeq == !!outeq) |>
+ dplyr::mutate(time = time + offset)
+
+ label <- "Target"
+ }
+
+ if (nrow(actual) == 0 || nrow(sim) == 0) {
+ return(NULL)
+ }
+
+ pred <- bd_weighted_pred_summary(x, sim)
+ if (is.null(pred) || nrow(pred) == 0) {
+ return(NULL)
+ }
+
+ 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))
+
+ if (nrow(fit) == 0) {
+ return(NULL)
+ }
+
+ fit
+}
+
+bd_fit_metrics <- function(data) {
+ if (is.null(data) || nrow(data) == 0) {
+ pe <- data.frame(type = rep(NA_character_, 8), absolute = rep(NA_real_, 8), percent = rep(NA_real_, 8))
+ return(list(pe = pe, metric_info = get_metric_info(pe)))
+ }
+
+ N <- nrow(data)
+ mae <- sum(data$d, na.rm = TRUE) / N
+ percent_mae <- mean(data$d / data$obs, na.rm = TRUE) * 100
+
+ mwe <- sum(data$wd, na.rm = TRUE) / N
+ percent_mwe <- sum(data$wd, na.rm = TRUE) / sum(data$obs / data$obsSD, na.rm = TRUE) * 100
+
+ mean_obs <- mean(data$obs, na.rm = TRUE)
+ wmean_obs <- sum(data$obs / data$obsSD, na.rm = TRUE) / N
+
+ mse <- sum(data$ds, na.rm = TRUE) / N
+ percent_mse <- mean(data$ds, na.rm = TRUE) / (mean_obs^2) * 100
+
+ mwse <- sum(data$wds, na.rm = TRUE) / N
+ percent_mwse <- mwse / (wmean_obs^2) * 100
+
+ rmse <- sqrt(mse)
+ percent_rmse <- rmse / mean_obs * 100
+
+ mbase <- mse - mae^2
+ percent_mbase <- mbase / (mean_obs^2) * 100
+
+ mbawse <- mwse - mwe^2
+ percent_mbawse <- mbawse / (wmean_obs^2) * 100
+
+ rmbawse <- sqrt(mbawse)
+ percent_rmbawse <- rmbawse * 100 / wmean_obs
+
+ pe <- data.frame(
+ type = c("mae", "mwe", "mse", "mwse", "rmse", "mbase", "mbawse", "rmbawse"),
+ absolute = c(mae, mwe, mse, mwse, rmse, mbase, mbawse, rmbawse),
+ percent = c(percent_mae, percent_mwe, percent_mse, percent_mwse, percent_rmse, percent_mbase, percent_mbawse, percent_rmbawse)
+ )
+
+ list(pe = pe, metric_info = get_metric_info(pe))
+}
+
+bd_fit_plot <- function(data, title) {
+ if (is.null(data) || nrow(data) == 0) {
+ return(NULL)
+ }
+
+ lims <- range(c(data$pred, data$obs), na.rm = TRUE)
+ fit <- if (nrow(data) >= 2) stats::lm(obs ~ pred, data = data) else NULL
+ ann_text <- NULL
+ if (!is.null(fit) && !anyNA(stats::coef(fit))) {
+ ann_text <- sprintf(
+ "R² = %.3f
Intercept = %.3f
Slope = %.3f",
+ summary(fit)$r.squared,
+ stats::coef(fit)[1],
+ stats::coef(fit)[2]
+ )
+ }
+
+ p <- plotly::plot_ly(
+ data = data, x = ~pred, y = ~obs, type = "scatter", mode = "markers",
+ marker = list(size = 9, color = "rgba(214, 86, 69, 0.75)", line = list(color = "rgba(120, 40, 30, 0.9)", width = 1)),
+ hovertemplate = "Predicted: %{x}
Observed: %{y}"
+ ) |>
+ plotly::layout(
+ title = title,
+ xaxis = list(title = "Predicted", range = lims),
+ yaxis = list(title = "Observed", range = lims),
+ shapes = list(list(
+ type = "line",
+ x0 = lims[1], x1 = lims[2],
+ y0 = lims[1], y1 = lims[2],
+ line = list(color = "rgba(80, 80, 80, 0.8)", dash = "dash")
+ )),
+ hovermode = "closest"
+ )
+
+ if (!is.null(fit) && !anyNA(stats::coef(fit))) {
+ reg_line <- data.frame(x = lims, y = stats::coef(fit)[1] + stats::coef(fit)[2] * lims)
+ p <- p |>
+ plotly::add_lines(
+ data = reg_line, x = ~x, y = ~y, inherit = FALSE,
+ line = list(color = "#2c3e50", width = 2),
+ hoverinfo = "skip", showlegend = FALSE
+ ) |>
+ plotly::layout(annotations = list(list(
+ x = 0.02, y = 0.98, xref = "paper", yref = "paper",
+ text = ann_text, align = "left",
+ showarrow = FALSE,
+ bgcolor = "rgba(255,255,255,0.85)",
+ bordercolor = "rgba(44,62,80,0.4)",
+ font = list(size = 11)
+ )))
+ }
+
+ p
+}
+
+bd_percent_metric <- function(pe, method) {
+ if (is.null(pe) || !all(c("type", "percent") %in% names(pe))) {
+ return(NA_real_)
+ }
+
+ method0 <- gsub("^percent_", "", method)
+ idx <- which(pe$type == method0)
+ if (length(idx) == 0) {
+ return(NA_real_)
+ }
+
+ as.numeric(pe$percent[idx[1]])
+}
+
+bd_fit_rainbow_plot <- function(past_metrics, future_metrics) {
+ bias_method <- getPMoptions("bias_method")
+ imp_method <- getPMoptions("imp_method")
+
+ past_bias <- bd_percent_metric(past_metrics$pe, bias_method)
+ past_imp <- bd_percent_metric(past_metrics$pe, imp_method)
+ future_bias <- bd_percent_metric(future_metrics$pe, bias_method)
+ future_imp <- bd_percent_metric(future_metrics$pe, imp_method)
+
+ x_vals <- c(-100, 100, past_bias, future_bias)
+ y_vals <- c(0, 100, past_imp, future_imp)
+ x_vals <- x_vals[is.finite(x_vals)]
+ y_vals <- y_vals[is.finite(y_vals)]
+
+ if (length(x_vals) == 0 || length(y_vals) == 0) {
+ return(NULL)
+ }
+
+ xlim <- c(min(x_vals), max(x_vals))
+ ylim <- c(0, max(y_vals))
+
+ if (diff(xlim) < 40) xlim <- c(xlim[1] - 20, xlim[2] + 20)
+ if (ylim[2] < 20) ylim[2] <- 20
+
+ b <- seq(100, 0, -20)
+ a <- seq(50, 0, -10)
+ cols <- c(
+ "rgba(255,0,0,0.7)",
+ "rgba(255,165,0,0.7)",
+ "rgba(255,255,0,0.7)",
+ "rgba(154,205,50,0.7)",
+ "rgba(0,128,0,0.7)"
+ )
+
+ p <- plotly::plot_ly()
+
+ for (i in 1:5) {
+ a1 <- a[i]
+ a2 <- a[i + 1]
+ b1 <- b[i]
+ b2 <- b[i + 1]
+
+ x1 <- seq(-a1, a1, by = 1)
+ x2 <- seq(-a2, a2, by = 1)
+ y1 <- sqrt(pmax(0, b1^2 * (1 - x1^2 / a1^2)))
+ y2 <- sqrt(pmax(0, b2^2 * (1 - x2^2 / a2^2)))
+
+ poly_x <- c(x1, rev(x2))
+ poly_y <- c(y1, rev(y2))
+
+ p <- p |>
+ plotly::add_trace(
+ x = poly_x,
+ y = poly_y,
+ type = "scatter",
+ mode = "lines",
+ fill = "toself",
+ fillcolor = cols[i],
+ line = list(color = "rgba(0,0,0,0)", width = 0),
+ hoverinfo = "skip",
+ showlegend = FALSE
+ )
+ }
+
+ if (is.finite(past_bias) && is.finite(past_imp)) {
+ p <- p |>
+ plotly::add_markers(
+ x = past_bias,
+ y = past_imp,
+ marker = list(color = "black", size = 11, line = list(color = "white", width = 2)),
+ name = "Past (% bias, % imprecision)",
+ hoverlabel = list(bgcolor = bd_rainbow_region_color(past_bias, past_imp, alpha = 1)),
+ hovertemplate = "Past
%Bias: %{x:.2f}
%Imprecision: %{y:.2f}",
+ showlegend = TRUE
+ )
+ }
+
+ if (is.finite(future_bias) && is.finite(future_imp)) {
+ p <- p |>
+ plotly::add_markers(
+ x = future_bias,
+ y = future_imp,
+ marker = list(color = "red", size = 11, line = list(color = "white", width = 2)),
+ name = "Future (% bias, % imprecision)",
+ hoverlabel = list(bgcolor = bd_rainbow_region_color(future_bias, future_imp, alpha = 1)),
+ hovertemplate = "Future
%Bias: %{x:.2f}
%Imprecision: %{y:.2f}",
+ showlegend = TRUE
+ )
+ }
+
+ p |>
+ plotly::layout(
+ title = "Fit quality rainbow",
+ xaxis = list(title = "% Bias", range = xlim, zeroline = FALSE),
+ yaxis = list(title = "% Imprecision", range = ylim, zeroline = FALSE),
+ shapes = list(
+ list(
+ type = "rect",
+ x0 = xlim[1],
+ x1 = xlim[2],
+ y0 = 0,
+ y1 = ylim[2],
+ xref = "x",
+ yref = "y",
+ fillcolor = "#808080",
+ line = list(width = 0),
+ layer = "below"
+ ),
+ list(
+ type = "line",
+ x0 = 0,
+ x1 = 0,
+ y0 = 0,
+ y1 = ylim[2],
+ line = list(color = "rgba(128,128,128,0.9)", dash = "dash"),
+ layer = "below"
+ )
+ ),
+ legend = list(orientation = "h", y = -0.2)
+ )
+}
+
+bd_rainbow_region_color <- function(bias, imp, alpha = 1) {
+ if (!is.finite(bias) || !is.finite(imp)) {
+ return(if (alpha >= 1) "#808080" else "rgba(128,128,128,0.7)")
+ }
+
+ a <- seq(50, 0, -10)
+ b <- seq(100, 0, -20)
+
+ inside_ellipse <- function(x, y, ax, by) {
+ if (!is.finite(ax) || !is.finite(by) || ax <= 0 || by <= 0) {
+ return(FALSE)
+ }
+ ((x^2) / (ax^2) + (y^2) / (by^2)) <= 1
+ }
+
+ base_cols <- c("#FF0000", "#FFA500", "#FFFF00", "#9ACD32", "#008000")
+ rgba_cols <- c(
+ "rgba(255,0,0,0.7)",
+ "rgba(255,165,0,0.7)",
+ "rgba(255,255,0,0.7)",
+ "rgba(154,205,50,0.7)",
+ "rgba(0,128,0,0.7)"
+ )
+
+ # smallest enclosing ellipse determines visible region color
+ for (i in 5:1) {
+ if (inside_ellipse(bias, imp, a[i], b[i])) {
+ return(if (alpha >= 1) base_cols[i] else rgba_cols[i])
+ }
+ }
+
+ if (alpha >= 1) "#808080" else "rgba(128,128,128,0.7)"
+}
+
+bd_fit_rainbow_plot_single <- function(metrics, color = "black", name = "Point", title = "Fit quality") {
+ bias_method <- getPMoptions("bias_method")
+ imp_method <- getPMoptions("imp_method")
+
+ bias_val <- bd_percent_metric(metrics$pe, bias_method)
+ imp_val <- bd_percent_metric(metrics$pe, imp_method)
+
+ x_vals <- c(-100, 100)
+ y_vals <- c(0, 100)
+ if (is.finite(bias_val)) x_vals <- c(x_vals, bias_val)
+ if (is.finite(imp_val)) y_vals <- c(y_vals, imp_val)
+
+ xlim <- c(min(x_vals), max(x_vals))
+ ylim <- c(0, max(y_vals))
+
+ if (diff(xlim) < 40) xlim <- c(xlim[1] - 20, xlim[2] + 20)
+ if (ylim[2] < 20) ylim[2] <- 20
+
+ b <- seq(100, 0, -20)
+ a <- seq(50, 0, -10)
+ cols <- c(
+ "rgba(255,0,0,0.7)",
+ "rgba(255,165,0,0.7)",
+ "rgba(255,255,0,0.7)",
+ "rgba(154,205,50,0.7)",
+ "rgba(0,128,0,0.7)"
+ )
+
+ p <- plotly::plot_ly()
+
+ for (i in 1:5) {
+ a1 <- a[i]
+ a2 <- a[i + 1]
+ b1 <- b[i]
+ b2 <- b[i + 1]
+ x1 <- seq(-a1, a1, by = 1)
+ x2 <- seq(-a2, a2, by = 1)
+ y1 <- sqrt(pmax(0, b1^2 * (1 - x1^2 / a1^2)))
+ y2 <- sqrt(pmax(0, b2^2 * (1 - x2^2 / a2^2)))
+ poly_x <- c(x1, rev(x2))
+ poly_y <- c(y1, rev(y2))
+ p <- p |> plotly::add_trace(
+ x = poly_x, y = poly_y,
+ type = "scatter", mode = "lines",
+ fill = "toself", fillcolor = cols[i],
+ line = list(color = "rgba(0,0,0,0)", width = 0),
+ hoverinfo = "skip", showlegend = FALSE
+ )
+ }
+
+ if (is.finite(bias_val) && is.finite(imp_val)) {
+ p <- p |> plotly::add_markers(
+ x = bias_val, y = imp_val,
+ marker = list(color = color, size = 11, line = list(color = "white", width = 2)),
+ name = name,
+ hoverlabel = list(bgcolor = bd_rainbow_region_color(bias_val, imp_val, alpha = 1)),
+ hovertemplate = paste0(name, "
%Bias: %{x:.2f}
%Imprecision: %{y:.2f}"),
+ showlegend = TRUE
+ )
+ }
+
+ p |> plotly::layout(
+ title = list(text = ""),
+ xaxis = list(title = "% Bias", range = xlim, zeroline = FALSE),
+ yaxis = list(title = "% Imprecision", range = ylim, zeroline = FALSE),
+ shapes = list(
+ list(
+ type = "rect",
+ x0 = xlim[1], x1 = xlim[2],
+ y0 = 0, y1 = ylim[2],
+ xref = "x", yref = "y",
+ fillcolor = "#808080",
+ line = list(width = 0),
+ layer = "below"
+ ),
+ list(
+ type = "line",
+ x0 = 0, x1 = 0,
+ y0 = 0, y1 = ylim[2],
+ line = list(color = "rgba(128,128,128,0.9)", dash = "dash"),
+ layer = "below"
+ )
+ ),
+ legend = list(orientation = "h", y = -0.2)
+ )
+}
+
+bd_sim_weighted_conc <- function(x, sim_obj, outeq = 1, time_offset = 0) {
+ if (is.null(sim_obj) || is.null(sim_obj$data) || is.null(sim_obj$data$obs)) {
+ return(NULL)
+ }
+
+ sim <- sim_obj$data$obs |>
+ dplyr::filter(outeq == !!outeq) |>
+ dplyr::mutate(time = time + time_offset)
+
+ pred <- bd_weighted_pred_summary(x, sim)
+ if (is.null(pred) || nrow(pred) == 0) {
+ return(NULL)
+ }
+
+ pred |>
+ dplyr::group_by(time) |>
+ dplyr::summarise(conc = mean(pred, na.rm = TRUE), .groups = "drop") |>
+ dplyr::arrange(time)
+}
+
+bd_past_dose_events <- function(x) {
+ if (is.null(x$past) || is.null(x$past$standard_data)) {
+ return(NULL)
+ }
+
+ doses <- x$past$standard_data |>
+ dplyr::filter(evid == 1)
+
+ if (nrow(doses) == 0) {
+ return(NULL)
+ }
+
+ doses <- bd_expand_dose_events(doses)
+ if (nrow(doses) == 0) {
+ return(NULL)
+ }
+
+ doses |>
+ dplyr::mutate(source = "Past")
+}
+
+bd_auc_interval <- function(conc_data, t0, t1) {
+ seg <- conc_data |>
+ dplyr::filter(time >= t0, time <= t1) |>
+ dplyr::arrange(time)
+
+ if (nrow(seg) < 2) {
+ return(0)
+ }
+
+ dt <- diff(seg$time)
+ y0 <- seg$conc[-nrow(seg)]
+ y1 <- seg$conc[-1]
+ sum((y0 + y1) / 2 * dt, na.rm = TRUE)
+}
+
+bd_auc_table <- function(x, outeq = 1, start_datetime = NULL) {
+ offset <- bd_future_time_offset(x)
+
+ past_conc <- bd_sim_weighted_conc(x, x$past_pred, outeq = outeq, time_offset = 0)
+ future_conc <- bd_sim_weighted_conc(x, x$future_pred, outeq = outeq, time_offset = offset)
+ conc <- dplyr::bind_rows(past_conc, future_conc)
+
+ if (is.null(conc) || nrow(conc) == 0) {
+ return(NULL)
+ }
+
+ conc <- conc |>
+ dplyr::group_by(time) |>
+ dplyr::summarise(conc = mean(conc, na.rm = TRUE), .groups = "drop") |>
+ dplyr::arrange(time)
+
+ past_doses <- bd_past_dose_events(x)
+ future_doses <- bd_future_dose_events(x)
+ if (!is.null(future_doses) && nrow(future_doses) > 0) {
+ future_doses <- future_doses |>
+ dplyr::mutate(source = "Future")
+ }
+
+ doses <- dplyr::bind_rows(past_doses, future_doses)
+ if (is.null(doses) || nrow(doses) == 0) {
+ return(NULL)
+ }
+
+ doses <- doses |>
+ dplyr::arrange(time)
+
+ n <- nrow(doses)
+ postdose <- numeric(n)
+ final_time <- max(conc$time, na.rm = TRUE)
+
+ for (i in seq_len(n)) {
+ t0 <- doses$time[i]
+ t1 <- if (i < n) doses$time[i + 1] else final_time
+ if (is.finite(t0) && is.finite(t1) && t1 > t0) {
+ postdose[i] <- bd_auc_interval(conc, t0, t1)
+ } else {
+ postdose[i] <- 0
+ }
+ }
+
+ out <- doses |>
+ dplyr::transmute(
+ time = time,
+ dose = dose,
+ source = source,
+ postdose_auc = postdose,
+ cumulative_auc = cumsum(postdose)
+ )
+
+ bd_add_datetime_column(out, start_datetime)
+}
+
+bd_auc_plot <- function(auc_tbl) {
+ if (is.null(auc_tbl) || nrow(auc_tbl) == 0) {
+ return(NULL)
+ }
+
+ date_fmt <- getPMoptions("date_format", warn = FALSE, quiet = TRUE)
+ if (!is.character(date_fmt) || length(date_fmt) != 1 || !nzchar(date_fmt)) {
+ date_fmt <- if (grepl("en_US", Sys.getlocale("LC_TIME"), fixed = TRUE)) "%m/%d/%y" else "%d/%m/%y"
+ }
+
+ df <- auc_tbl |> dplyr::mutate(
+ datetime = suppressWarnings(as.POSIXct(datetime, tz = Sys.timezone())),
+ dose_idx = seq_len(dplyr::n()),
+ dose_amt = ifelse(is.na(dose), NA_character_, formatC(as.numeric(dose), format = "f", digits = 2)),
+ date_label = format(datetime, date_fmt),
+ time_label = format(datetime, "%H:%M:%S")
+ )
+
+ # colors per source
+ cols <- ifelse(df$source == "Past", "black", "red")
+
+ p <- plotly::plot_ly()
+
+ # bars for post-dose AUC
+ p <- p |> plotly::add_bars(
+ x = ~ df$datetime, y = ~ df$postdose_auc,
+ marker = list(color = cols),
+ opacity = 0.5,
+ name = "Post-dose AUC",
+ showlegend = FALSE,
+ text = paste0(
+ "Dose: ", df$dose_idx,
+ "
Amount: ", df$dose_amt,
+ "
Date: ", df$date_label,
+ "
Time: ", df$time_label,
+ "
Post-dose AUC: ", sprintf("%.2f", df$postdose_auc)
+ ),
+ hoverinfo = "text",
+ textposition = "none"
+ )
+
+ # cumulative line split into past / future to change color
+ past_idx <- which(df$source == "Past")
+ future_idx <- which(df$source != "Past")
+
+ if (length(past_idx) > 0) {
+ p <- p |> plotly::add_trace(
+ x = df$datetime[past_idx], y = df$cumulative_auc[past_idx],
+ type = "scatter", mode = "lines+markers", name = "Cumulative (Past)",
+ line = list(color = "black"), marker = list(color = "black"), yaxis = "y2", showlegend = FALSE,
+ text = paste0(
+ "Dose: ", df$dose_idx[past_idx],
+ "
Amount: ", df$dose_amt[past_idx],
+ "
Date: ", df$date_label[past_idx],
+ "
Time: ", df$time_label[past_idx],
+ "
Cumulative AUC: ", sprintf("%.2f", df$cumulative_auc[past_idx])
+ ),
+ hoverinfo = "text"
+ )
+ }
+
+ if (length(future_idx) > 0) {
+ p <- p |> plotly::add_trace(
+ x = df$datetime[future_idx], y = df$cumulative_auc[future_idx],
+ type = "scatter", mode = "lines+markers", name = "Cumulative (Future)",
+ line = list(color = "red"), marker = list(color = "red"), yaxis = "y2", showlegend = FALSE,
+ text = paste0(
+ "Dose: ", df$dose_idx[future_idx],
+ "
Amount: ", df$dose_amt[future_idx],
+ "
Date: ", df$date_label[future_idx],
+ "
Time: ", df$time_label[future_idx],
+ "
Cumulative AUC: ", sprintf("%.2f", df$cumulative_auc[future_idx])
+ ),
+ hoverinfo = "text"
+ )
+ }
+
+ # bridge segment from last past to first future cumulative point
+ if (length(past_idx) > 0 && length(future_idx) > 0) {
+ i_last_past <- max(past_idx)
+ i_first_future <- min(future_idx)
+ p <- p |> plotly::add_trace(
+ x = c(df$datetime[i_last_past], df$datetime[i_first_future]),
+ y = c(df$cumulative_auc[i_last_past], df$cumulative_auc[i_first_future]),
+ type = "scatter", mode = "lines", yaxis = "y2",
+ line = list(color = "black", width = 2),
+ showlegend = FALSE,
+ hoverinfo = "skip"
+ )
+ }
+
+ # x-axis tick text: dose number only
+ ticktext <- paste0("Dose ", df$dose_idx)
+
+ p |> plotly::layout(
+ xaxis = list(title = "", tickvals = df$datetime, ticktext = ticktext),
+ yaxis = list(title = "Post-dose AUC", zeroline = FALSE),
+ yaxis2 = list(
+ overlaying = "y", side = "right",
+ title = list(text = "Cumulative AUC", standoff = 40),
+ automargin = TRUE,
+ zeroline = FALSE
+ ),
+ margin = list(t = 90),
+ bargap = 0.2,
+ showlegend = FALSE
+ )
+}
+
+bd_pd_table <- function(x, outeq = 1, start_datetime = NULL) {
+ if (is.null(x$future) || is.null(x$future$standard_data)) {
+ return(NULL)
+ }
+
+ target_rows <- x$future$standard_data |>
+ dplyr::filter(evid == 0, outeq == !!outeq)
+ if (nrow(target_rows) == 0) {
+ return(NULL)
+ }
+
+ target_time <- as.numeric(target_rows$time[1])
+ target_value <- as.numeric(target_rows$out[1])
+ if (!is.finite(target_time) || !is.finite(target_value)) {
+ return(NULL)
+ }
+
+ past_doses <- bd_past_dose_events(x)
+ future_doses <- bd_future_dose_events(x)
+ if (!is.null(future_doses) && nrow(future_doses) > 0) {
+ future_doses <- future_doses |>
+ dplyr::mutate(source = "Future")
+ }
+
+ doses <- dplyr::bind_rows(past_doses, future_doses)
+ if (is.null(doses) || nrow(doses) == 0) {
+ return(NULL)
+ }
+
+ doses <- doses |>
+ dplyr::arrange(time) |>
+ dplyr::mutate(eval_time = time + target_time)
+
+ offset <- bd_future_time_offset(x)
+ past_sim <- NULL
+ future_sim <- NULL
+
+ if (!is.null(x$past_pred) && !is.null(x$past_pred$data) && !is.null(x$past_pred$data$obs)) {
+ past_sim <- x$past_pred$data$obs |>
+ dplyr::filter(outeq == !!outeq)
+ }
+ if (!is.null(x$future_pred) && !is.null(x$future_pred$data) && !is.null(x$future_pred$data$obs)) {
+ future_sim <- x$future_pred$data$obs |>
+ dplyr::filter(outeq == !!outeq) |>
+ dplyr::mutate(time = time + offset)
+ }
+
+ sim <- dplyr::bind_rows(past_sim, future_sim)
+ if (is.null(sim) || nrow(sim) == 0) {
+ return(NULL)
+ }
+
+ sim <- sim |>
+ dplyr::mutate(nsim = as.character(nsim))
+
+ weights <- bd_weight_table(x, sim$nsim) |>
+ dplyr::mutate(nsim = as.character(nsim))
+
+ sim_split <- split(sim, sim$nsim)
+ nsim_keys <- intersect(names(sim_split), weights$nsim)
+ if (length(nsim_keys) == 0) {
+ return(NULL)
+ }
+
+ eval_times <- doses$eval_time
+ success_mat <- matrix(FALSE, nrow = length(nsim_keys), ncol = length(eval_times))
+ w <- numeric(length(nsim_keys))
+
+ for (j in seq_along(nsim_keys)) {
+ key <- nsim_keys[j]
+ d <- sim_split[[key]] |>
+ dplyr::arrange(time)
+ d <- d[!is.na(d$time) & !is.na(d$out), , drop = FALSE]
+ if (nrow(d) < 2) next
+
+ pred_at <- stats::approx(
+ x = d$time,
+ y = d$out,
+ xout = eval_times,
+ rule = 2,
+ ties = mean
+ )$y
+
+ success <- pred_at >= target_value
+ success_mat[j, ] <- as.logical(success)
+
+ w_row <- weights$weight[weights$nsim == key]
+ w[j] <- if (length(w_row) == 0 || !is.finite(w_row[1])) 0 else as.numeric(w_row[1])
+ }
+
+ if (sum(w, na.rm = TRUE) <= 0) {
+ return(NULL)
+ }
+
+ cum_prob <- vapply(seq_along(eval_times), function(i) {
+ stats::weighted.mean(as.numeric(success_mat[, i]), w = w, na.rm = TRUE)
+ }, numeric(1))
+
+ out <- doses |>
+ dplyr::transmute(
+ dose_number = dplyr::row_number(),
+ time = time,
+ source = source,
+ cumulative_probability = cum_prob
+ )
+
+ bd_add_datetime_column(out, start_datetime)
+}
+
+bd_add_datetime_column <- function(data, start_datetime) {
+ if (is.null(data) || nrow(data) == 0 || is.null(start_datetime) || is.na(start_datetime) || !"time" %in% names(data)) {
+ return(data)
+ }
+
+ data |>
+ dplyr::mutate(datetime = format(start_datetime + time * 3600, "%Y-%m-%d %H:%M:%S")) |>
+ dplyr::relocate(datetime, .after = time)
+}
+
+bd_normalize_weights <- function(w) {
+ w <- as.numeric(w)
+ w[!is.finite(w)] <- 0
+ s <- sum(w, na.rm = TRUE)
+ if (!is.finite(s) || s <= 0) {
+ return(rep(0, length(w)))
+ }
+ w / s
+}
+
+bd_parameter_shift_data <- function(x, bins = 25) {
+ if (is.null(x$posterior) || is.null(x$posterior$theta)) {
+ return(NULL)
+ }
+
+ theta <- as.data.frame(x$posterior$theta)
+ if (nrow(theta) == 0 || ncol(theta) == 0) {
+ return(NULL)
+ }
+
+ param_names <- colnames(theta)
+ if (is.null(param_names) || !length(param_names)) {
+ param_names <- x$posterior$param_names
+ colnames(theta) <- param_names
+ }
+
+ prior_w <- bd_normalize_weights(x$posterior$population_weights)
+ post_w <- bd_normalize_weights(x$posterior$posterior_weights)
+
+ n_points <- nrow(theta)
+ if (length(prior_w) != n_points) prior_w <- rep(1 / n_points, n_points)
+ if (length(post_w) != n_points) post_w <- rep(1 / n_points, n_points)
+
+ bins <- max(5L, as.integer(bins))
+
+ long <- lapply(seq_along(param_names), function(i) {
+ p <- param_names[i]
+ vals <- as.numeric(theta[[p]])
+
+ ok <- is.finite(vals)
+ if (!any(ok)) {
+ return(NULL)
+ }
+
+ vals_ok <- vals[ok]
+ prior_ok <- bd_normalize_weights(prior_w[ok])
+ post_ok <- bd_normalize_weights(post_w[ok])
+
+ if (length(vals_ok) == 1 || diff(range(vals_ok, na.rm = TRUE)) == 0) {
+ bin_idx <- rep(1L, length(vals_ok))
+ breaks <- range(vals_ok, na.rm = TRUE)
+ if (length(breaks) < 2 || !is.finite(diff(breaks)) || diff(breaks) == 0) {
+ breaks <- c(vals_ok[1] - 0.5, vals_ok[1] + 0.5)
+ }
+ bins_use <- 1L
+ } else {
+ breaks <- seq(min(vals_ok, na.rm = TRUE), max(vals_ok, na.rm = TRUE), length.out = bins + 1L)
+ bin_idx <- cut(vals_ok, breaks = breaks, include.lowest = TRUE, labels = FALSE)
+ bins_use <- bins
+ }
+
+ bin_tbl <- tibble::tibble(
+ parameter = p,
+ bin = as.integer(bin_idx),
+ value = vals_ok,
+ prior_weight = prior_ok,
+ post_weight = post_ok
+ ) |>
+ dplyr::group_by(parameter, bin) |>
+ dplyr::summarise(
+ value_low = min(value, na.rm = TRUE),
+ value_high = max(value, na.rm = TRUE),
+ value_mid = stats::weighted.mean(value, w = post_weight, na.rm = TRUE),
+ prior_prob = sum(prior_weight, na.rm = TRUE),
+ post_prob = sum(post_weight, na.rm = TRUE),
+ .groups = "drop"
+ )
+
+ all_bins <- tibble::tibble(bin = seq_len(bins_use))
+ out <- all_bins |>
+ dplyr::left_join(bin_tbl, by = dplyr::join_by(bin)) |>
+ dplyr::mutate(
+ parameter = p,
+ prior_prob = dplyr::coalesce(prior_prob, 0),
+ post_prob = dplyr::coalesce(post_prob, 0),
+ delta_prob = post_prob - prior_prob,
+ bin_pct = (bin - 0.5) / bins_use * 100,
+ q_label = sprintf("Q%02d", bin)
+ )
+
+ out
+ }) |>
+ dplyr::bind_rows()
+
+ if (is.null(long) || nrow(long) == 0) {
+ return(NULL)
+ }
+
+ summary_tbl <- long |>
+ dplyr::group_by(parameter) |>
+ dplyr::summarise(
+ total_variation = 0.5 * sum(abs(delta_prob), na.rm = TRUE),
+ max_shift = max(abs(delta_prob), na.rm = TRUE),
+ .groups = "drop"
+ ) |>
+ dplyr::arrange(dplyr::desc(total_variation), dplyr::desc(max_shift))
+
+ means_tbl <- tibble::tibble(
+ parameter = param_names,
+ prior_mean = vapply(seq_along(param_names), function(i) {
+ stats::weighted.mean(as.numeric(theta[[i]]), w = prior_w, na.rm = TRUE)
+ }, numeric(1)),
+ posterior_mean = vapply(seq_along(param_names), function(i) {
+ stats::weighted.mean(as.numeric(theta[[i]]), w = post_w, na.rm = TRUE)
+ }, numeric(1))
+ ) |>
+ dplyr::mutate(mean_shift = posterior_mean - prior_mean)
+
+ summary_tbl <- summary_tbl |>
+ dplyr::left_join(means_tbl, by = dplyr::join_by(parameter))
+
+ long <- long |>
+ dplyr::mutate(
+ parameter = factor(parameter, levels = summary_tbl$parameter)
+ )
+
+ list(long = long, summary = summary_tbl, bins = bins)
+}
+
+bd_parameter_shift_heatmap <- function(shift_data) {
+ if (is.null(shift_data) || is.null(shift_data$long) || nrow(shift_data$long) == 0) {
+ return(NULL)
+ }
+
+ df <- shift_data$long |>
+ dplyr::arrange(parameter, bin)
+
+ df <- df |>
+ dplyr::mutate(
+ x_value = dplyr::coalesce(value_mid, (value_low + value_high) / 2)
+ )
+
+ df_prior <- df |>
+ dplyr::filter(is.finite(x_value), is.finite(prior_prob), prior_prob > 0)
+
+ df_post <- df |>
+ dplyr::filter(is.finite(x_value), is.finite(post_prob), post_prob > 0)
+
+ gg <- ggplot2::ggplot() +
+ ggplot2::geom_point(
+ data = df_prior,
+ ggplot2::aes(x = x_value, y = prior_prob),
+ color = "grey55",
+ size = 2.3,
+ shape = 16
+ ) +
+ ggplot2::geom_point(
+ data = df_post,
+ ggplot2::aes(x = x_value, y = post_prob),
+ color = "red",
+ alpha = 0.5,
+ size = 2.3,
+ shape = 16
+ ) +
+ ggplot2::facet_wrap(~parameter, scales = "free") +
+ ggplot2::scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) +
+ ggplot2::scale_y_continuous(breaks = scales::pretty_breaks(n = 5)) +
+ ggplot2::labs(
+ title = "Marginal parameter distributions: prior (gray) vs posterior (red)",
+ x = "Parameter value",
+ y = "Probability"
+ ) +
+ ggplot2::theme_minimal(base_size = 12) +
+ ggplot2::theme(
+ legend.position = "none",
+ panel.grid.minor = ggplot2::element_blank()
+ )
+
+ plotly::ggplotly(gg)
+}
+
+bd_report_build <- function(x, outeq = 1) {
+ start_datetime <- bd_plot_start_datetime(x)
+ future_offset <- bd_future_time_offset(x)
+
+ overview <- tibble::tibble(
+ item = c("Status", "Method", "Objective", "log(Objective)", "Prior weight", "Start", "Start offset (h)", "Future starts at (h)"),
+ value = c(
+ x$status %||% NA_character_,
+ x$method %||% NA_character_,
+ signif(x$objf %||% NA_real_, 6),
+ if (is.null(x$objf) || is.na(x$objf) || x$objf <= 0) NA_real_ else signif(log(x$objf), 6),
+ x$prior_weight %||% NA_real_,
+ paste(x$start %||% NA_character_, collapse = ", "),
+ x$start_offset %||% NA_real_,
+ future_offset
+ )
+ )
+
+ past_last_dose <- NULL
+ past_last_obs <- NULL
+ if (!is.null(x$past) && !is.null(x$past$standard_data) && nrow(x$past$standard_data) > 0) {
+ past_last_dose <- x$past$standard_data |>
+ dplyr::filter(evid == 1) |>
+ dplyr::arrange(dplyr::desc(time)) |>
+ dplyr::slice(1) |>
+ dplyr::select(time, dose, dur, input)
+
+ past_last_obs <- x$past$standard_data |>
+ dplyr::filter(evid == 0, outeq == !!outeq) |>
+ dplyr::arrange(dplyr::desc(time)) |>
+ dplyr::slice(1) |>
+ dplyr::select(time, out, outeq)
+ }
+
+ # Original requested doses (0 = Optimized, non-zero = Fixed), captured before optimization
+ original_requested_doses <- x$future_requested_doses
+
+ future_doses <- bd_future_dose_events(x)
+ if (!is.null(future_doses) && nrow(future_doses) > 0) {
+ n_fd <- nrow(future_doses)
+ status_vec <- rep("Fixed", n_fd)
+ if (!is.null(original_requested_doses) && length(original_requested_doses) >= n_fd) {
+ status_vec <- dplyr::if_else(
+ as.numeric(original_requested_doses[seq_len(n_fd)]) == 0,
+ "Optimized", "Fixed"
+ )
+ }
+ future_doses <- future_doses |>
+ dplyr::select(time, dose, dur, input) |>
+ dplyr::mutate(status = status_vec) |>
+ bd_add_datetime_column(start_datetime)
+ }
+
+ future_targets <- NULL
+ if (!is.null(x$future) && !is.null(x$future$standard_data)) {
+ future_targets <- x$future$standard_data |>
+ dplyr::filter(evid == 0, outeq == !!outeq) |>
+ dplyr::mutate(time = time + future_offset) |>
+ dplyr::select(time, out, outeq) |>
+ bd_add_datetime_column(start_datetime)
+ }
+
+ future_predictions <- x$result$predictions
+ if (!is.null(future_predictions) && nrow(future_predictions) > 0) {
+ future_predictions <- future_predictions |>
+ dplyr::filter(outeq == !!(outeq - 1)) |>
+ bd_add_datetime_column(start_datetime)
+ }
+
+ auc_table <- bd_auc_table(x, outeq = outeq, start_datetime = start_datetime)
+ auc_plot <- bd_auc_plot(auc_table)
+ pd_table <- bd_pd_table(x, outeq = outeq, start_datetime = start_datetime)
+
+ pd_target_summary <- NULL
+ if (!is.null(x$future_target_info)) {
+ ti <- x$future_target_info
+ tt <- ti$target_type %||% "concentration"
+ tgt <- ti$target
+ ttm <- ti$target_time
+ if (!is.null(tt) && tt == "time") {
+ pct <- if (!is.null(ttm) && is.numeric(ttm) && !is.na(ttm)) round(ttm * 100) else NA
+ pd_target_summary <- paste0(
+ "Therapeutic target: Time above ", tgt,
+ " for ", pct, "% of each dosing interval."
+ )
+ } else {
+ type_label <- if (!is.null(tt) && tt == "auc") "AUC" else "Concentration"
+ ttm_display <- if (!is.null(ttm) && length(ttm) == 1 && is.numeric(ttm) && !is.na(ttm)) ttm else 24
+ pd_target_summary <- paste0(
+ "Therapeutic target: ", type_label, " of ", tgt,
+ ", ", ttm_display, " hours after each dose."
+ )
+ }
+ }
+
+ past_fit <- bd_fit_data(x, source = "past", outeq = outeq)
+ future_fit <- bd_fit_data(x, source = "future", outeq = outeq)
+
+ past_fit <- bd_add_datetime_column(past_fit, start_datetime)
+ future_fit <- bd_add_datetime_column(future_fit, start_datetime)
+
+ past_fit_metrics <- bd_fit_metrics(past_fit)
+ future_fit_metrics <- bd_fit_metrics(future_fit)
+
+ posterior_table <- as.data.frame(x$posterior$theta)
+ posterior_table$weight <- as.numeric(x$posterior$posterior_weights)
+ posterior_table$population_weight <- as.numeric(x$posterior$population_weights)
+
+ posterior_summary <- tibble::tibble(
+ parameter = x$posterior$param_names,
+ weighted_mean = vapply(seq_along(x$posterior$param_names), function(i) {
+ stats::weighted.mean(x$posterior$theta[, i], x$posterior$posterior_weights, na.rm = TRUE)
+ }, numeric(1)),
+ weighted_median = vapply(seq_along(x$posterior$param_names), function(i) {
+ vals <- x$posterior$theta[, i]
+ wts <- x$posterior$posterior_weights
+ ord <- order(vals)
+ vals <- vals[ord]
+ wts <- wts[ord] / sum(wts)
+ vals[which.max(cumsum(wts) >= 0.5)]
+ }, numeric(1)),
+ weighted_sd = vapply(seq_along(x$posterior$param_names), function(i) {
+ vals <- x$posterior$theta[, i]
+ mu <- stats::weighted.mean(vals, x$posterior$posterior_weights, na.rm = TRUE)
+ sqrt(stats::weighted.mean((vals - mu)^2, x$posterior$posterior_weights, na.rm = TRUE))
+ }, numeric(1))
+ )
+
+ shift_data <- bd_parameter_shift_data(x, bins = 25)
+
+ list(
+ overview = overview,
+ plot = plot.bd(x, quiet = TRUE, print = FALSE, future_region = FALSE),
+ past_rainbow_plot = bd_fit_rainbow_plot_single(past_fit_metrics, color = "black", name = "Past"),
+ future_rainbow_plot = bd_fit_rainbow_plot_single(future_fit_metrics, color = "red", name = "Future"),
+ past_last_dose = bd_add_datetime_column(past_last_dose, start_datetime),
+ past_last_obs = bd_add_datetime_column(past_last_obs, start_datetime),
+ future_doses = future_doses,
+ future_targets = future_targets,
+ future_predictions = future_predictions,
+ auc_plot = auc_plot,
+ auc_table = auc_table,
+ pd_table = pd_table,
+ pd_target_summary = pd_target_summary,
+ past_fit = past_fit,
+ future_fit = future_fit,
+ past_fit_metrics = past_fit_metrics,
+ future_fit_metrics = future_fit_metrics,
+ past_fit_plot = bd_fit_plot(past_fit, "Past observations vs predictions"),
+ future_fit_plot = bd_fit_plot(future_fit, "Future targets vs predictions"),
+ posterior_summary = posterior_summary,
+ posterior_table = posterior_table,
+ parameter_shift_plot = bd_parameter_shift_heatmap(shift_data),
+ parameter_shift_summary = if (!is.null(shift_data)) shift_data$summary else NULL,
+ model_summary = paste(capture.output(print(x$posterior$model_info$model)), collapse = "\n")
+ )
+}
+
+bd_report <- function(x, path, show = TRUE, quiet = TRUE, title = "BestDose Report") {
+ if (!inherits(x, "bd")) {
+ cli::cli_abort(c("x" = "This function expects a {.cls bd} object."))
+ }
+
+ template_file <- system.file("report/templates/bestdose.Rmd", package = "Pmetrics")
+ if (!nzchar(template_file) || !file.exists(template_file)) {
+ cli::cli_abort(c("x" = "BestDose report template was not found in the package."))
+ }
+
+ out_path <- if (missing(path)) tempdir() else normalizePath(path, winslash = "/", mustWork = FALSE)
+ fs::dir_create(out_path)
+
+ rmarkdown::render(
+ input = template_file,
+ output_file = file.path(out_path, "bestdose_report.html"),
+ params = list(bd = x, title = title),
+ clean = TRUE,
+ quiet = quiet
+ )
+
+ out_file <- file.path(out_path, "bestdose_report.html")
+ if (file.exists(out_file)) {
+ if (show) {
+ pander::openFileInOS(out_file)
+ }
+ return(invisible(1))
+ }
+
+ invisible(-1)
+}
diff --git a/R/Pmetrics-package.R b/R/Pmetrics-package.R
index ad023234..14d94a7d 100755
--- a/R/Pmetrics-package.R
+++ b/R/Pmetrics-package.R
@@ -10,6 +10,7 @@
#' @aliases Pmetrics
#' @author Michael Neely MD
#' \url{http://www.lapk.org}
+#' @useDynLib Pmetrics, .registration = TRUE
#' @keywords internal
#' "_PACKAGE"
#'
diff --git a/R/extendr-wrappers.R b/R/extendr-wrappers.R
index 20be6617..b38258d3 100755
--- a/R/extendr-wrappers.R
+++ b/R/extendr-wrappers.R
@@ -1,22 +1,12 @@
-# Generated by extendr: Do not edit by hand
+#' Simulates the first subject in the data set using the model at the given path.
# nolint start
-#
-# This file was created with the following call:
-# .Call("wrap__make_Pmetrics_wrappers", use_symbols = TRUE, package_name = "Pmetrics")
-
-#' @usage NULL
-#' @useDynLib Pmetrics, .registration = TRUE
-NULL
-
-#' Simulates the first subject in the data set using the model at the given path.
#' @param data_path Path to the data file.
#' @param model_path Path to the compiled model file.
#' @param spp One support point as a numeric vector with probabiltity.
#' @param kind Kind of model, which can either be "ODE" or "Analytical".
#' @return Simulation results.
-#' @keywords internal
#'@export
simulate_one <- function(data_path, model_path, spp, kind) .Call(wrap__simulate_one, data_path, model_path, spp, kind)
@@ -26,7 +16,6 @@ simulate_one <- function(data_path, model_path, spp, kind) .Call(wrap__simulate_
#' @param theta Data frame of support points.
#' @param kind Kind of model, which can either be "ODE" or "Analytical".
#' @return Simulation results.
-#' @keywords internal
#' @export
simulate_all <- function(data_path, model_path, theta, kind) .Call(wrap__simulate_all, data_path, model_path, theta, kind)
@@ -37,20 +26,17 @@ simulate_all <- function(data_path, model_path, theta, kind) .Call(wrap__simulat
#' @param template_path Path to the template directory.
#' @param kind Kind of model, which can either be "ODE" or "Analytical".
#' @return Result of the compilation process.
-#' @keywords internal
#' @export
compile_model <- function(model_path, output_path, params, template_path, kind) .Call(wrap__compile_model, model_path, output_path, params, template_path, kind)
#' Dummy function to cache compilation artifacts.
#' @param template_path Path to the template directory.
#' @return Path to the build directory.
-#' @keywords internal
#' @export
dummy_compile <- function(template_path) .Call(wrap__dummy_compile, template_path)
#' Checks if Cargo is installed on the system.
#' @return TRUE if Cargo is installed, FALSE otherwise.
-#' @keywords internal
#' @export
is_cargo_installed <- function() .Call(wrap__is_cargo_installed)
@@ -61,7 +47,6 @@ is_cargo_installed <- function() .Call(wrap__is_cargo_installed)
#' @param output_path Path to save the fitting results.
#' @param kind Kind of model, which can either be "ODE" or "Analytical".
#' @return Result of the fitting process.
-#' @keywords internal
#' @export
fit <- function(model_path, data, params, output_path, kind) .Call(wrap__fit, model_path, data, params, output_path, kind)
@@ -69,13 +54,11 @@ fit <- function(model_path, data, params, output_path, kind) .Call(wrap__fit, mo
#' @param model_path Path to the compiled model file.
#' @param kind Kind of model, which can either be "ODE" or "Analytical".
#' @return List of model parameters.
-#' @keywords internal
#' @export
model_parameters <- function(model_path, kind) .Call(wrap__model_parameters, model_path, kind)
#' Retrieves the temporary path used for building models.
#' @return Temporary build path.
-#' @keywords internal
#' @export
temporary_path <- function() .Call(wrap__temporary_path)
@@ -84,5 +67,9 @@ temporary_path <- function() .Call(wrap__temporary_path)
#' @export
setup_logs <- function() .Call(wrap__setup_logs)
+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)
+
# nolint end
diff --git a/_pkgdown.yml b/_pkgdown.yml
index b3a9865d..847f6714 100755
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -51,6 +51,8 @@ reference:
- PM_load
- PM_result
- PM_valid
+ - bd
+ - bd_post
- title: Project management
desc: Functions for project management and workflow
- contents:
diff --git a/inst/Examples/Rscript/bestdose_simple_test.R b/inst/Examples/Rscript/bestdose_simple_test.R
new file mode 100644
index 00000000..58cd4cce
--- /dev/null
+++ b/inst/Examples/Rscript/bestdose_simple_test.R
@@ -0,0 +1,148 @@
+devtools::load_all()
+
+mod_onecomp <- PM_model$new(
+ pri = list(
+ ke = ab(0.001, 3.0),
+ v = ab(25.0, 250.0)
+ ),
+ eqn = function() {
+ dx[1] <- -ke * X[1] + B[1]
+ },
+ out = function() {
+ Y[1] <- X[1] / v
+ },
+ err = list(
+ additive(0, c(0, 0.20, 0, 0))
+ )
+)
+
+past_file <- "inst/Examples/src/bestdose_past.csv"
+target_file <- "inst/Examples/src/bestdose_target.csv"
+prior_file <- "inst/Examples/src/bestdose_prior.csv"
+
+cat("Using PMcore prior file\n")
+
+posterior <- bd_post$new(
+ prior = prior_file,
+ model = mod_onecomp,
+ past_data = PM_data$new(past_file, quiet = TRUE),
+ max_cycles = 500,
+ quiet = TRUE
+)
+
+cat("\nPosterior theta (first 5 rows):\n")
+print(head(posterior$theta, 5))
+cat("\nPosterior weights (first 10):\n")
+print(head(posterior$posterior_weights, 10))
+
+prior_weights <- seq(0, 1, by = 0.1)
+for (lambda in prior_weights) {
+ r <- posterior$optimize(
+ target = target_file,
+ dose_range = list(min = 0, max = 300),
+ prior_weight = lambda, # ought to be called "prior_weight"
+ start = NULL
+ )
+ cat(sprintf(
+ "Prior weight: %.2f\t\tOptimal dose: [%s]\t\tCost: %.6f\t\tln Cost: %.4f\t\tMethod: %s\n",
+ lambda,
+ paste(r$doses, collapse = ", "),
+ r$objf,
+ log(r$objf),
+ r$method
+ ))
+}
+
+# # Last result predictions
+# r <- posterior$optimize(
+# target = target_file,
+# dose_range = list(min = 0, max = 300),
+# prior_weight = 0.0
+# )
+# cat("\nConcentration-time predictions for prior_weight=0.0:\n")
+# preds <- r$result$predictions
+# for (j in seq_len(nrow(preds))) {
+# p <- preds[j, ]
+# cat(sprintf(
+# "Time: %.2f h, Observed: %.2f, (Pop Mean: %.4f, Pop Median: %.4f, Post Mean: %.4f, Post Median: %.4f)\n",
+# p$time, p$obs, p$pop_mean, p$pop_median, p$post_mean, p$post_median
+# ))
+# }
+
+
+# # ===== One-shot API =====
+# bd1 <- bd$new(
+# prior = prior_file,
+# model = mod_onecomp,
+# past_data = PM_data$new(past_file, quiet = TRUE),
+# max_cycles = 500,
+# future = PM_data$new(target_file, quiet = TRUE),
+# dose_range = list(min = 0, max = 5000),
+# prior_weight = 0,
+# start = 0
+# )
+
+# bd1
+
+# bd1$plot()
+
+# bd2 <- bd$new(
+# prior = prior_file,
+# model = mod_onecomp,
+# max_cycles = 500,
+# future = target_file,
+# dose_range = list(min = 0, max = 300),
+# prior_weight = 0.0
+# )
+
+# bd2
+
+
+
+# bd2$plot()
+
+# plot(bd1)
+
+
+# ## adding ability to specify future as argument to bd$new() instead of target file.
+# ## This will allow for more flexible future specifications and avoid the need for a separate target file.
+
+# future_list <- list(dose = c(120,0),
+# frequency = 12,
+# route = 0,
+# number = 3,
+# target_time = 11.5,
+# target = 0.3,
+# target_type = "concentration")
+
+# future_list <- list(dose = 320,
+# frequency = 12,
+# route = 0,
+# number = 3,
+# target_time = 0.6,
+# target = 0.3,
+# target_type = "time")
+
+# #AUC target - crashing as of 3/13/26
+# future_list <- list(dose = 320,
+# frequency = 12,
+# route = 0,
+# number = 3,
+# target_time = 11.5,
+# target = 10,
+# target_type = "auc")
+
+
+# bd_new <- bd$new(
+# prior = prior_file,
+# model = mod_onecomp,
+# past_data = PM_data$new(past_file, quiet = TRUE),
+# max_cycles = 50,
+# future = future_list,
+# dose_range = list(min = 0, max = 5000),
+# prior_weight = 0,
+# start = "02/02/26 00:00",
+# )
+
+# # bd_new$plot()
+# bd_new$report()
diff --git a/inst/Examples/Rscript/run_bestdose_comparison.R b/inst/Examples/Rscript/run_bestdose_comparison.R
new file mode 100644
index 00000000..f6888324
--- /dev/null
+++ b/inst/Examples/Rscript/run_bestdose_comparison.R
@@ -0,0 +1,71 @@
+devtools::load_all()
+
+mod_onecomp <- PM_model$new(
+ pri = list(
+ ke = ab(0.001, 3.0),
+ v = ab(25.0, 250.0)
+ ),
+ eqn = function() {
+ dx[1] <- -ke * X[1] + B[1]
+ },
+ out = function() {
+ Y[1] <- X[1] / v
+ },
+ err = list(
+ additive(0, c(0, 0.20, 0, 0))
+ )
+)
+
+past_file <- "inst/Examples/src/bestdose_past.csv"
+target_file <- "inst/Examples/src/bestdose_target.csv"
+prior_file <- "inst/Examples/src/bestdose_prior.csv"
+
+cat("Using PMcore prior file\n")
+
+posterior <- bd_post$new(
+ prior = prior_file,
+ model = mod_onecomp,
+ past_data = PM_data$new(past_file, quiet = TRUE),
+ max_cycles = 500,
+ quiet = TRUE
+)
+
+cat("\nPosterior theta (first 5 rows):\n")
+print(head(posterior$theta, 5))
+cat("\nPosterior weights (first 10):\n")
+print(head(posterior$posterior_weights, 10))
+
+prior_weights <- seq(0, 1, by = 0.1)
+for (lambda in prior_weights) {
+ r <- posterior$optimize(
+ target = target_file,
+ dose_range = list(min = 0, max = 300),
+ prior_weight = lambda,
+ start = NULL
+ )
+ cat(sprintf(
+ "Prior weight: %.2f\t\tOptimal dose: [%.4f, %.4f]\t\tCost: %.6f\t\tln Cost: %.4f\t\tMethod: %s\n",
+ lambda,
+ r$doses[1], r$doses[2],
+ r$objf,
+ log(r$objf),
+ r$method
+ ))
+}
+
+# Last result predictions
+r <- posterior$optimize(
+ target = target_file,
+ dose_range = list(min = 0, max = 300),
+ prior_weight = 1.0,
+ start = NULL
+)
+cat("\nConcentration-time predictions for prior_weight=1.0:\n")
+preds <- r$result$predictions
+for (j in seq_len(nrow(preds))) {
+ p <- preds[j, ]
+ cat(sprintf(
+ "Time: %.2f h, Observed: %.2f, (Pop Mean: %.4f, Pop Median: %.4f, Post Mean: %.4f, Post Median: %.4f)\n",
+ p$time, p$obs, p$pop_mean, p$pop_median, p$post_mean, p$post_median
+ ))
+}
diff --git a/inst/Examples/Rscript/test_future_vectors.R b/inst/Examples/Rscript/test_future_vectors.R
new file mode 100644
index 00000000..ecb2cdd8
--- /dev/null
+++ b/inst/Examples/Rscript/test_future_vectors.R
@@ -0,0 +1,91 @@
+devtools::load_all(quiet = TRUE)
+
+mod_onecomp <- PM_model$new(
+ pri = list(ke = ab(0.001, 3.0), v = ab(25.0, 250.0)),
+ eqn = function() { dx[1] <- -ke * X[1] + B[1] },
+ out = function() { Y[1] <- X[1] / v },
+ err = list(additive(0, c(0, 0.20, 0, 0)))
+)
+
+# Test 1: scalar (legacy behavior)
+cat("--- Test 1: scalar dose/freq/route, number=3 ---\n")
+bd1 <- bd$new(
+ prior = "inst/Examples/src/bestdose_prior.csv",
+ model = mod_onecomp,
+ past_data = PM_data$new("inst/Examples/src/bestdose_past.csv", quiet = TRUE),
+ max_cycles = 5,
+ future = list(dose = 1000, frequency = 12, route = 0, number = 3, target_time = 11.5, target = 0.3, target_type = "concentration"),
+ dose_range = list(min = 0, max = 5000), prior_weight = 0, start = "02/02/26 00:00",
+ quiet = TRUE
+)
+cat("doses:", paste(bd1$future$data$dose[!is.na(bd1$future$data$dose)], collapse = ","), "\n")
+cat("times:", paste(bd1$future$data$time[!is.na(bd1$future$data$dose)], collapse = ","), "\n\n")
+
+# Test 2: vector dose (3 elements, number=3)
+cat("--- Test 2: vector dose c(500,750,1000), number=3 ---\n")
+bd2 <- bd$new(
+ prior = "inst/Examples/src/bestdose_prior.csv",
+ model = mod_onecomp,
+ past_data = PM_data$new("inst/Examples/src/bestdose_past.csv", quiet = TRUE),
+ max_cycles = 5,
+ future = list(dose = c(500, 750, 1000), frequency = 12, route = 0, number = 3, target_time = 11.5, target = 0.3, target_type = "concentration"),
+ dose_range = list(min = 0, max = 5000), prior_weight = 0, start = "02/02/26 00:00",
+ quiet = TRUE
+)
+cat("doses:", paste(bd2$future$data$dose[!is.na(bd2$future$data$dose)], collapse = ","), "\n")
+cat("times:", paste(bd2$future$data$time[!is.na(bd2$future$data$dose)], collapse = ","), "\n\n")
+
+# Test 3: number < vector length - should inform and set number=3
+cat("--- Test 3: number=2 < dose length 3 (should inform, expand to 3) ---\n")
+bd3 <- bd$new(
+ prior = "inst/Examples/src/bestdose_prior.csv",
+ model = mod_onecomp,
+ past_data = PM_data$new("inst/Examples/src/bestdose_past.csv", quiet = TRUE),
+ max_cycles = 5,
+ future = list(dose = c(500, 750, 1000), frequency = 12, route = 0, number = 2, target_time = 11.5, target = 0.3, target_type = "concentration"),
+ dose_range = list(min = 0, max = 5000), prior_weight = 0, start = "02/02/26 00:00",
+ quiet = TRUE
+)
+cat("doses:", paste(bd3$future$data$dose[!is.na(bd3$future$data$dose)], collapse = ","), "\n\n")
+
+# Test 4: short dose vector recycled (dose length 2, number=4 -> last value 1000 recycled)
+cat("--- Test 4: dose c(500,1000), number=4 (recycle last -> 500,1000,1000,1000) ---\n")
+bd4 <- bd$new(
+ prior = "inst/Examples/src/bestdose_prior.csv",
+ model = mod_onecomp,
+ past_data = PM_data$new("inst/Examples/src/bestdose_past.csv", quiet = TRUE),
+ max_cycles = 5,
+ future = list(dose = c(500, 1000), frequency = 12, route = 0, number = 4, target_time = 11.5, target = 0.3, target_type = "concentration"),
+ dose_range = list(min = 0, max = 5000), prior_weight = 0, start = "02/02/26 00:00",
+ quiet = TRUE
+)
+cat("doses:", paste(bd4$future$data$dose[!is.na(bd4$future$data$dose)], collapse = ","), "\n\n")
+
+# Test 5: vector frequency (variable intervals)
+cat("--- Test 5: dose=1000, freq=c(8,12,24), number=3 ---\n")
+bd5 <- bd$new(
+ prior = "inst/Examples/src/bestdose_prior.csv",
+ model = mod_onecomp,
+ past_data = PM_data$new("inst/Examples/src/bestdose_past.csv", quiet = TRUE),
+ max_cycles = 5,
+ future = list(dose = 1000, frequency = c(8, 12, 24), route = 0, number = 3, target_time = 7.5, target = 0.3, target_type = "concentration"),
+ dose_range = list(min = 0, max = 5000), prior_weight = 0, start = "02/02/26 00:00",
+ quiet = TRUE
+)
+cat("doses:", paste(bd5$future$data$dose[!is.na(bd5$future$data$dose)], collapse = ","), "\n")
+cat("times:", paste(bd5$future$data$time[!is.na(bd5$future$data$dose)], collapse = ","), "\n\n")
+
+# Test 6: target_type = time (proportional target_time per interval)
+cat("--- Test 6: target_type='time', freq=c(8,12,24), target_time=0.5 ---\n")
+bd6 <- bd$new(
+ prior = "inst/Examples/src/bestdose_prior.csv",
+ model = mod_onecomp,
+ past_data = PM_data$new("inst/Examples/src/bestdose_past.csv", quiet = TRUE),
+ max_cycles = 5,
+ future = list(dose = 1000, frequency = c(8, 12, 24), route = 0, number = 3, target_time = 0.5, target = 0.3, target_type = "time"),
+ dose_range = list(min = 0, max = 5000), prior_weight = 0, start = "02/02/26 00:00",
+ quiet = TRUE
+)
+cat("obs times:", paste(bd6$future$data$time[is.na(bd6$future$data$dose)], collapse = ","), "\n\n")
+
+cat("All tests completed.\n")
diff --git a/inst/Examples/Rscript/test_status_column.R b/inst/Examples/Rscript/test_status_column.R
new file mode 100644
index 00000000..079c8c70
--- /dev/null
+++ b/inst/Examples/Rscript/test_status_column.R
@@ -0,0 +1,35 @@
+devtools::load_all(quiet = TRUE)
+
+mod_onecomp <- PM_model$new(
+ pri = list(ke = ab(0.001, 3.0), v = ab(25.0, 250.0)),
+ eqn = function() { dx[1] <- -ke * X[1] + B[1] },
+ out = function() { Y[1] <- X[1] / v },
+ err = list(additive(0, c(0, 0.20, 0, 0)))
+)
+
+# Mixed: dose 1 = fixed (1000), dose 2 = optimize (0), dose 3 = fixed (1500)
+bd_obj <- bd$new(
+ prior = "inst/Examples/src/bestdose_prior.csv",
+ model = mod_onecomp,
+ past_data = PM_data$new("inst/Examples/src/bestdose_past.csv", quiet = TRUE),
+ max_cycles = 5,
+ future = list(dose = c(1000, 0, 1500), frequency = 12, route = 0, number = 3,
+ target_time = 11.5, target = 0.3, target_type = "concentration"),
+ dose_range = list(min = 0, max = 5000), prior_weight = 0, start = "02/02/26 00:00",
+ quiet = TRUE
+)
+
+# Debug: check what standard_data looks like for future evid==1 rows
+cat("x$future$standard_data evid==1 rows:\n")
+sd <- bd_obj$future$standard_data
+print(sd[sd$evid == 1, ])
+
+report <- Pmetrics:::bd_report_build(bd_obj)
+cat("future_doses status column:\n")
+print(report$future_doses[, c("dose", "status")])
+
+out <- file.path(tempdir(), "bd-report-status-test")
+dir.create(out, showWarnings = FALSE, recursive = TRUE)
+status <- bd_obj$report(path = out, show = FALSE, quiet = TRUE)
+cat("\nReport status:", status, "\n")
+cat("HTML exists:", file.exists(file.path(out, "bestdose_report.html")), "\n")
diff --git a/inst/Examples/src/bestdose_past.csv b/inst/Examples/src/bestdose_past.csv
new file mode 100644
index 00000000..a367c173
--- /dev/null
+++ b/inst/Examples/src/bestdose_past.csv
@@ -0,0 +1,9 @@
+id,evid,date,time,dur,dose,addl,ii,input,out,outeq,c0,c1,c2,c3
+1,1,2/1/26,00:00:00,0,150,0,0,1,.,.,.,.,.,.
+1,0,2/1/26,02:00:00,.,.,.,.,.,0.759050698,1,.,.,.,.
+1,0,2/1/26,04:00:00,.,.,.,.,.,0.38408517,1,.,.,.,.
+1,0,2/1/26,06:00:00,.,.,.,.,.,0.194349887,1,.,.,.,.
+1,1,2/1/26,12:00:00,0,75,0,0,1,.,.,.,.,.,.
+1,0,2/1/26,14:00:00,.,.,.,.,.,0.392266578,1,.,.,.,.
+1,0,2/1/26,16:00:00,.,.,.,.,.,0.198489739,1,.,.,.,.
+1,0,2/1/26,18:00:00,.,.,.,.,.,0.100437251,1,.,.,.,.
\ No newline at end of file
diff --git a/inst/Examples/src/bestdose_prior.csv b/inst/Examples/src/bestdose_prior.csv
new file mode 100644
index 00000000..bead3957
--- /dev/null
+++ b/inst/Examples/src/bestdose_prior.csv
@@ -0,0 +1,50 @@
+ke,v,prob
+0.08736658442020416,103.9818823337555,0.06405457820397967
+1.0041342623233795,188.76309752464294,0.019607843140974805
+0.10507606186866761,140.8398061990738,0.019607843150332413
+0.3220049516916275,68.28039824962616,0.019607843051963586
+0.3493796042442322,75.56532084941864,0.01960746473012483
+0.04363290762901306,85.5460512638092,0.019607220383661656
+0.3312444769859314,92.81326055526733,0.019607147791086518
+0.09591740503311158,70.4040265083313,0.019606957258885387
+0.3085237291574478,92.42967247962952,0.019739312212632944
+0.3319534166574478,118.09373497962952,0.01963231510029087
+0.06250557525157929,107.03847944736481,0.019671831895535216
+0.06885235497951508,76.48277580738068,0.019704372081019293
+0.104416601395607,91.46958649158478,0.0196439196378291
+0.11730292952060699,78.46177399158478,0.019596592693258374
+0.06250557525157929,95.56875288486481,0.019558583620767178
+0.013994081234931946,178.73696267604828,0.01973719361931262
+0.07190197534561157,66.9762921333313,0.01951179533204779
+0.3406021231412888,99.86292123794556,0.0195017935173522
+0.02050767450332642,140.53142726421356,0.019496909863886087
+0.28171865940093993,140.80943822860718,0.02006966089395463
+0.08736658442020416,127.7562963962555,0.019975556860205728
+0.283145404958725,125.71949124336243,0.021760196860831894
+0.2902425238609314,94.35134649276733,0.020471376913953784
+0.3010685295820236,128.3288460969925,0.019431925066036485
+0.044804392004013066,99.5206606388092,0.021064943923371642
+0.09146677973270416,92.7318823337555,0.01938595576577497
+0.2925854926109314,128.32107305526733,0.018981650341481797
+0.28490263152122497,92.49683499336243,0.018729335263653696
+0.09826037378311157,111.1413311958313,0.022513324580852458
+0.08170614981651307,97.3233950138092,0.022071606402980712
+0.2925854926109314,105.42556524276733,0.02188055586035111
+0.036757444095611574,109.3395733833313,0.01813223949068482
+0.09240295190811158,113.7780499458313,0.01634642552944466
+0.1311608946800232,96.91233038902283,0.020745654510324375
+0.3540655417442322,86.77137553691864,0.03911312316643254
+0.3132096666574478,115.76463341712952,0.030954606713449555
+0.07658791284561157,113.1188702583313,0.00868552053488122
+0.2908282660484314,112.85232305526733,0.014634851770971591
+0.05615650746822358,130.38859486579895,0.0012110863756026463
+0.08678084223270416,87.3705542087555,0.00043536953019302594
+0.057327991843223576,130.87199330329895,0.004552363447102815
+0.08736658442020416,87.4144995212555,0.0002790654992634474
+0.07658791284561157,113.1628155708313,0.026588971195136306
+0.3540655417442322,86.81532084941864,0.0001850932799277153
+0.3132096666574478,115.80857872962952,0.008825717688183198
+0.2908282660484314,112.80837774276733,0.03942289250196689
+0.1311608946800232,96.86838507652283,0.018470031787835278
+0.05615650746822357,130.43254017829895,0.05302517206894798
+0.08678084223270416,87.4144995212555,0.0349542088912624
diff --git a/inst/Examples/src/bestdose_prior_pmcore.csv b/inst/Examples/src/bestdose_prior_pmcore.csv
new file mode 100644
index 00000000..bead3957
--- /dev/null
+++ b/inst/Examples/src/bestdose_prior_pmcore.csv
@@ -0,0 +1,50 @@
+ke,v,prob
+0.08736658442020416,103.9818823337555,0.06405457820397967
+1.0041342623233795,188.76309752464294,0.019607843140974805
+0.10507606186866761,140.8398061990738,0.019607843150332413
+0.3220049516916275,68.28039824962616,0.019607843051963586
+0.3493796042442322,75.56532084941864,0.01960746473012483
+0.04363290762901306,85.5460512638092,0.019607220383661656
+0.3312444769859314,92.81326055526733,0.019607147791086518
+0.09591740503311158,70.4040265083313,0.019606957258885387
+0.3085237291574478,92.42967247962952,0.019739312212632944
+0.3319534166574478,118.09373497962952,0.01963231510029087
+0.06250557525157929,107.03847944736481,0.019671831895535216
+0.06885235497951508,76.48277580738068,0.019704372081019293
+0.104416601395607,91.46958649158478,0.0196439196378291
+0.11730292952060699,78.46177399158478,0.019596592693258374
+0.06250557525157929,95.56875288486481,0.019558583620767178
+0.013994081234931946,178.73696267604828,0.01973719361931262
+0.07190197534561157,66.9762921333313,0.01951179533204779
+0.3406021231412888,99.86292123794556,0.0195017935173522
+0.02050767450332642,140.53142726421356,0.019496909863886087
+0.28171865940093993,140.80943822860718,0.02006966089395463
+0.08736658442020416,127.7562963962555,0.019975556860205728
+0.283145404958725,125.71949124336243,0.021760196860831894
+0.2902425238609314,94.35134649276733,0.020471376913953784
+0.3010685295820236,128.3288460969925,0.019431925066036485
+0.044804392004013066,99.5206606388092,0.021064943923371642
+0.09146677973270416,92.7318823337555,0.01938595576577497
+0.2925854926109314,128.32107305526733,0.018981650341481797
+0.28490263152122497,92.49683499336243,0.018729335263653696
+0.09826037378311157,111.1413311958313,0.022513324580852458
+0.08170614981651307,97.3233950138092,0.022071606402980712
+0.2925854926109314,105.42556524276733,0.02188055586035111
+0.036757444095611574,109.3395733833313,0.01813223949068482
+0.09240295190811158,113.7780499458313,0.01634642552944466
+0.1311608946800232,96.91233038902283,0.020745654510324375
+0.3540655417442322,86.77137553691864,0.03911312316643254
+0.3132096666574478,115.76463341712952,0.030954606713449555
+0.07658791284561157,113.1188702583313,0.00868552053488122
+0.2908282660484314,112.85232305526733,0.014634851770971591
+0.05615650746822358,130.38859486579895,0.0012110863756026463
+0.08678084223270416,87.3705542087555,0.00043536953019302594
+0.057327991843223576,130.87199330329895,0.004552363447102815
+0.08736658442020416,87.4144995212555,0.0002790654992634474
+0.07658791284561157,113.1628155708313,0.026588971195136306
+0.3540655417442322,86.81532084941864,0.0001850932799277153
+0.3132096666574478,115.80857872962952,0.008825717688183198
+0.2908282660484314,112.80837774276733,0.03942289250196689
+0.1311608946800232,96.86838507652283,0.018470031787835278
+0.05615650746822357,130.43254017829895,0.05302517206894798
+0.08678084223270416,87.4144995212555,0.0349542088912624
diff --git a/inst/Examples/src/bestdose_target.csv b/inst/Examples/src/bestdose_target.csv
new file mode 100644
index 00000000..79883dbe
--- /dev/null
+++ b/inst/Examples/src/bestdose_target.csv
@@ -0,0 +1,9 @@
+id,evid,time,dur,dose,addl,ii,input,out,outeq,c0,c1,c2,c3
+1,1,0,0,0,0,0,1,.,.,.,.,.,.
+1,0,2,.,.,.,.,.,0.759050697604428,1,.,.,.,.
+1,0,4,.,.,.,.,.,0.384085169721793,1,.,.,.,.
+1,0,6,.,.,.,.,.,0.194349887386702,1,.,.,.,.
+1,1,12,0,0,0,0,1,.,.,.,.,.,.
+1,0,14,.,.,.,.,.,0.392266577540038,1,.,.,.,.
+1,0,16,.,.,.,.,.,0.198489739204705,1,.,.,.,.
+1,0,18,.,.,.,.,.,0.100437250648841,1,.,.,.,.
diff --git a/inst/Examples/src/bestdose_target_pmcore.csv b/inst/Examples/src/bestdose_target_pmcore.csv
new file mode 100644
index 00000000..79883dbe
--- /dev/null
+++ b/inst/Examples/src/bestdose_target_pmcore.csv
@@ -0,0 +1,9 @@
+id,evid,time,dur,dose,addl,ii,input,out,outeq,c0,c1,c2,c3
+1,1,0,0,0,0,0,1,.,.,.,.,.,.
+1,0,2,.,.,.,.,.,0.759050697604428,1,.,.,.,.
+1,0,4,.,.,.,.,.,0.384085169721793,1,.,.,.,.
+1,0,6,.,.,.,.,.,0.194349887386702,1,.,.,.,.
+1,1,12,0,0,0,0,1,.,.,.,.,.,.
+1,0,14,.,.,.,.,.,0.392266577540038,1,.,.,.,.
+1,0,16,.,.,.,.,.,0.198489739204705,1,.,.,.,.
+1,0,18,.,.,.,.,.,0.100437250648841,1,.,.,.,.
diff --git a/man/bd.Rd b/man/bd.Rd
new file mode 100644
index 00000000..66fdf525
--- /dev/null
+++ b/man/bd.Rd
@@ -0,0 +1,243 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PM_bestdose.R
+\name{bd}
+\alias{bd}
+\title{Object to contain BestDose optimization results}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
+
+This object is created after a successful BestDose optimization run.
+BestDose finds optimal dosing regimens to achieve target drug concentrations
+or AUC values using Bayesian optimization.
+}
+\section{Public fields}{
+\if{html}{\out{}}
+\describe{
+\item{\code{past}}{PM_data object containing past patient data used in the optimization (if applicable)}
+
+\item{\code{past_pred}}{PM_sim object containing model predictions for the past data}
+
+\item{\code{future}}{PM_data object containing future data used in the optimization}
+
+\item{\code{future_pred}}{PM_sim object containing model predictions for the future data}
+
+\item{\code{result}}{List containing optimization results, including optimal doses, predictions, and objective function value}
+
+\item{\code{posterior}}{The \code{bd_post} object used to compute the posterior distribution (if applicable)}
+
+\item{\code{prior_weight}}{The prior weight (lambda) used in the optimization (if applicable)}
+
+\item{\code{start}}{Start specification provided by user (numeric hours or datetime string)}
+
+\item{\code{start_offset}}{Resolved numeric offset (hours) used internally; NULL means use target times as-is}
+
+\item{\code{future_requested_doses}}{Numeric vector of originally requested future doses (0 = optimize, non-zero = fixed)}
+
+\item{\code{future_target_info}}{List with target_type, target, and target_time from the original future specification}
+}
+\if{html}{\out{
}}
+}
+\section{Active bindings}{
+\if{html}{\out{}}
+\describe{
+\item{\code{doses}}{Optimal doses found by BestDose optimization}
+
+\item{\code{objf}}{Objective function value at the optimal doses}
+
+\item{\code{method}}{Optimization method used}
+
+\item{\code{status}}{Status of the optimization run}
+}
+\if{html}{\out{
}}
+}
+\section{Methods}{
+\subsection{Public methods}{
+\itemize{
+\item \href{#method-bd-new}{\code{bd$new()}}
+\item \href{#method-bd-print}{\code{bd$print()}}
+\item \href{#method-bd-save}{\code{bd$save()}}
+\item \href{#method-bd-plot}{\code{bd$plot()}}
+\item \href{#method-bd-report}{\code{bd$report()}}
+\item \href{#method-bd-clone}{\code{bd$clone()}}
+}
+}
+\if{html}{\out{
}}
+\if{html}{\out{}}
+\if{latex}{\out{\hypertarget{method-bd-new}{}}}
+\subsection{Method \code{new()}}{
+Initialize a \code{bd} object by running a one-shot BestDose optimization.
+Creates the posterior internally, then optimizes. For reusing the posterior
+across multiple optimizations, use \code{bd_post$new()} followed by \verb{$optimize()}.
+\subsection{Usage}{
+\if{html}{\out{}}\preformatted{bd$new(
+ model = NULL,
+ prior = NULL,
+ past_data = NULL,
+ future = NULL,
+ dose_range = list(min = 0, max = 1000),
+ prior_weight = 0.5,
+ start = 0,
+ max_cycles = 500,
+ settings = NULL,
+ posterior = NULL,
+ quiet = FALSE,
+ simulate = FALSE
+)}\if{html}{\out{
}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{}}
+\describe{
+\item{\code{model}}{PM_model object or path to compiled model}
+
+\item{\code{prior}}{Prior parameter value distributions for the model, can be a \link{PM_result}, \link{PM_final}, or path to \code{theta.csv} file.}
+
+\item{\code{past_data}}{PM_data object or path to CSV file with past patient data (optional). Omitting this will plan a new regimen for the patient
+based on the population prior without any individualization other than inclusion of any covariates which are in the model.}
+
+\item{\code{future}}{The future regimen to optimize. It can be one of three things:
+\itemize{
+\item A list describing future dosing/target setup with elements:
+\itemize{
+\item \code{dose} The amount to administer (required). Use \code{0} to optimize that dose to achieve the target; any other value is treated as fixed.
+Can be a scalar or a numeric vector. If a vector is provided, each element corresponds to each dose event in order.
+\item \code{frequency} How often to administer dose(s) (default 24 hours). Can be a scalar or numeric vector; when vectorized,
+element \code{i} is the interval after dose \code{i} and before dose \code{i+1}.
+\item \code{route} Either 0 for oral/bolus, or a value >0 indicating infusion duration. Can be a scalar or numeric vector;
+when vectorized, element \code{i} applies to dose \code{i}.
+\item \code{number} The total number of doses to administer in the future (default 1).
+If \code{number} is smaller than \code{max(length(dose), length(frequency), length(route))}, it is increased to that maximum and an informational message is shown.
+If any of \code{dose}, \code{frequency}, or \code{route} are shorter than \code{number}, their last value is recycled to length \code{number}.
+\item \code{target} The target value to achieve after every dose (required).
+\item \code{target_type} The target type (default \code{"concentration"}). Must be one of \code{"concentration"}, \code{"auc"}, or \code{"time"}.
+\code{"concentration"} targets a concentration at \code{target_time}; \code{"auc"} targets an AUC from dose time to \code{target_time};
+\code{"time"} uses \code{target_time} as a proportion between 0 and 1 of each dose interval.
+\item \code{target_time} For \code{target_type = "concentration"} or \code{"auc"}, the time after each dose at which to evaluate the target,
+default 24 hours if unspecified. For \code{target_type = "time"}, this must be a single numeric proportion between 0 and 1,
+and each target observation time is set to \code{round(frequency[i] * target_time)} hours after dose \code{i}.
+\item \code{covariates} Named list of covariates in the model with values to use. If there are no covariates in the model, this can be omitted.
+}
+\item A \code{PM_data} object containing the future regimen and target observations. Dose events should be included with \code{evid=1} and target observations should be included with \code{evid=0} and \code{out} equal to the target value.
+When \code{future} is supplied this way, concentration targets are assumed.
+\item A path to a CSV file containing the future regimen and target observations in the same format as described for the \code{PM_data} object above.
+When \code{future} is supplied this way, concentration targets are assumed.
+}}
+
+\item{\code{dose_range}}{List with 'min' and 'max' elements defining the dose search range (default: 0 to 1000)}
+
+\item{\code{prior_weight}}{Numeric between 0 and 1 indicating the balance between the prior and the posterior in the optimization (default: 0.5).
+If 1, the model prior parameter
+value distribution will be used as the sole basis for optimization. If 0, only the Bayesian posterior parameter distribution will be used,
+which incorporates the patient history. Values between 0 and 1 will weight the contribution of the prior vs posterior parameter distributions in the optimization.
+Choose values closer to 1 when the patient history is sparse or not believed to be very informative, and values closer to 0 when the patient history is rich and believed to be highly informative.}
+
+\item{\code{start}}{Start time for the future regimen. Can be either:
+\itemize{
+\item \code{NULL}: use the supplied future/target times as-is, without appending them after past data.
+\item Numeric hours (default \code{0}): relative to the last event date/time in past data, or if no past data,
+relative to the next local clock hour.
+\item Date-time character string: interpreted as an absolute start date/time and converted to an internal hour offset.
+Ambiguous formats (e.g. \code{01/02/26}) are resolved using the \code{date_format} Pmetrics option
+(set via \link{setPMoptions}), which defaults to \verb{\%m/\%d/\%y} for US locales and \verb{\%d/\%m/\%y} otherwise.
+Accepted formats include:
+\itemize{
+\item \verb{YYYY-mm-dd HH:MM[:SS]}
+\item \verb{mm/dd/[YY]YY HH:MM[:SS]} or \verb{mm-dd-[YY]YY HH:MM[:SS]}
+\item \verb{dd/mm/[YY]YY HH:MM[:SS]} or \verb{dd-mm-[YY]YY HH:MM[:SS]}
+}
+}}
+
+\item{\code{max_cycles}}{Maximum number of optimization cycles for computing the posterior (default: 500)}
+
+\item{\code{settings}}{List of additional settings for posterior computation (optional)}
+
+\item{\code{posterior}}{\code{bd_post} object to use instead of computing a new one (optional).
+When provided, \code{prior}, \code{model}, \code{past_data}, \code{max_cycles}, and \code{settings} are ignored, as they are already included in the \code{posterior} object.}
+
+\item{\code{quiet}}{Logical indicating whether to suppress verbose simulation output (default: FALSE)}
+}
+\if{html}{\out{
}}
+}
+\subsection{Returns}{
+A \code{bd} object containing the optimization results and associated information
+}
+}
+\if{html}{\out{
}}
+\if{html}{\out{}}
+\if{latex}{\out{\hypertarget{method-bd-print}{}}}
+\subsection{Method \code{print()}}{
+Print summary of BestDose results
+\subsection{Usage}{
+\if{html}{\out{}}\preformatted{bd$print()}\if{html}{\out{
}}
+}
+
+}
+\if{html}{\out{
}}
+\if{html}{\out{}}
+\if{latex}{\out{\hypertarget{method-bd-save}{}}}
+\subsection{Method \code{save()}}{
+Save results to RDS file
+\subsection{Usage}{
+\if{html}{\out{}}\preformatted{bd$save(filename = "bestdose_result.rds")}\if{html}{\out{
}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{}}
+\describe{
+\item{\code{filename}}{Path to save file. Default: "bestdose_result.rds"}
+}
+\if{html}{\out{
}}
+}
+}
+\if{html}{\out{
}}
+\if{html}{\out{}}
+\if{latex}{\out{\hypertarget{method-bd-plot}{}}}
+\subsection{Method \code{plot()}}{
+Plot observed and predicted concentrations over time for both past and future data
+\subsection{Usage}{
+\if{html}{\out{}}\preformatted{bd$plot(...)}\if{html}{\out{
}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{}}
+\describe{
+\item{\code{...}}{Additional arguments passed to plot.bd function}
+}
+\if{html}{\out{
}}
+}
+}
+\if{html}{\out{
}}
+\if{html}{\out{}}
+\if{latex}{\out{\hypertarget{method-bd-report}{}}}
+\subsection{Method \code{report()}}{
+Generate an HTML report for a BestDose result
+\subsection{Usage}{
+\if{html}{\out{}}\preformatted{bd$report(...)}\if{html}{\out{
}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{}}
+\describe{
+\item{\code{...}}{Parameters passed to the internal BestDose report generator.}
+}
+\if{html}{\out{
}}
+}
+}
+\if{html}{\out{
}}
+\if{html}{\out{}}
+\if{latex}{\out{\hypertarget{method-bd-clone}{}}}
+\subsection{Method \code{clone()}}{
+The objects of this class are cloneable with this method.
+\subsection{Usage}{
+\if{html}{\out{}}\preformatted{bd$clone(deep = FALSE)}\if{html}{\out{
}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{}}
+\describe{
+\item{\code{deep}}{Whether to make a deep clone.}
+}
+\if{html}{\out{
}}
+}
+}
+}
diff --git a/man/bd_post.Rd b/man/bd_post.Rd
new file mode 100644
index 00000000..35b20322
--- /dev/null
+++ b/man/bd_post.Rd
@@ -0,0 +1,141 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PM_bestdose.R
+\name{bd_post}
+\alias{bd_post}
+\title{Compute a reusable BestDose posterior}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}
+
+Use \code{bd_post} to compute the Bayesian posterior once from
+prior population data and patient history, then call \verb{$optimize()} multiple
+times with different targets, dose ranges, or bias weights.
+}
+\section{Public fields}{
+\if{html}{\out{}}
+\describe{
+\item{\code{handle}}{Memory pointer to the computed posterior (opaque to users)}
+
+\item{\code{theta}}{Matrix of support points in the posterior distribution}
+
+\item{\code{theta_dim}}{Dimensions of the theta matrix}
+
+\item{\code{param_names}}{Names of the parameters in the posterior}
+
+\item{\code{posterior_weights}}{Weights of the posterior support points}
+
+\item{\code{population_weights}}{Weights of the population support points}
+
+\item{\code{past}}{PM_data object containing past patient data used in the posterior computation (if applicable)}
+
+\item{\code{model_info}}{Information about the model used}
+
+\item{\code{settings}}{Settings used for the posterior computation}
+}
+\if{html}{\out{
}}
+}
+\section{Methods}{
+\subsection{Public methods}{
+\itemize{
+\item \href{#method-bd_post-new}{\code{bd_post$new()}}
+\item \href{#method-bd_post-optimize}{\code{bd_post$optimize()}}
+\item \href{#method-bd_post-clone}{\code{bd_post$clone()}}
+}
+}
+\if{html}{\out{
}}
+\if{html}{\out{}}
+\if{latex}{\out{\hypertarget{method-bd_post-new}{}}}
+\subsection{Method \code{new()}}{
+Initialize the \code{bd_post} object by computing the posterior distribution from the given prior, model, and past data
+\subsection{Usage}{
+\if{html}{\out{}}\preformatted{bd_post$new(
+ prior,
+ model,
+ past_data = NULL,
+ max_cycles = 500,
+ settings = NULL,
+ quiet = FALSE
+)}\if{html}{\out{
}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{}}
+\describe{
+\item{\code{prior}}{Prior information for the model, can be a PM_result, PM_final, or path to theta.csv}
+
+\item{\code{model}}{PM_model object or path to compiled model}
+
+\item{\code{past_data}}{PM_data object or path to CSV file with past patient data (optional)}
+
+\item{\code{max_cycles}}{Maximum number of optimization cycles for computing the posterior (default: 500)}
+
+\item{\code{settings}}{List of additional settings for posterior computation (optional)}
+
+\item{\code{quiet}}{Logical indicating whether to suppress verbose simulation output (default: FALSE)}
+}
+\if{html}{\out{
}}
+}
+\subsection{Returns}{
+A \code{bd_post} object containing the computed posterior distribution and associated information
+}
+}
+\if{html}{\out{
}}
+\if{html}{\out{}}
+\if{latex}{\out{\hypertarget{method-bd_post-optimize}{}}}
+\subsection{Method \code{optimize()}}{
+Run optimization and return a \code{bd} result object
+\subsection{Usage}{
+\if{html}{\out{}}\preformatted{bd_post$optimize(
+ target,
+ dose_range = list(min = 0, max = 1000),
+ prior_weight = 0.5,
+ start = 0,
+ quiet = NULL
+)}\if{html}{\out{
}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{}}
+\describe{
+\item{\code{target}}{Future regimen/target specification. Accepts the same forms as the \code{future} argument to \code{bd$new()}:
+a list with elements such as \code{dose}, \code{frequency}, \code{route}, \code{number}, \code{target}, \code{target_type}, and \code{target_time};
+a \code{PM_data} object; or a path to a CSV file.}
+
+\item{\code{dose_range}}{List with 'min' and 'max' elements defining the dose search range (default: 0 to 1000)}
+
+\item{\code{prior_weight}}{Numeric between 0 and 1 indicating the weight of bias in the optimization (default: 0.5)}
+
+\item{\code{start}}{Start time for future regimen. Can be either \code{NULL}, numeric hours (default: \code{0}), or a date-time
+character string. Ambiguous formats are resolved using the \code{date_format} Pmetrics option
+(see \link{setPMoptions}). Accepted string formats:
+\itemize{
+\item \verb{YYYY-mm-dd HH:MM[:SS]}
+\item \verb{mm/dd/[YY]YY HH:MM[:SS]} or \verb{mm-dd-[YY]YY HH:MM[:SS]}
+\item \verb{dd/mm/[YY]YY HH:MM[:SS]} or \verb{dd-mm-[YY]YY HH:MM[:SS]}
+}}
+
+\item{\code{quiet}}{Logical indicating whether to suppress verbose simulation output. If NULL, uses the quiet setting from posterior computation (default: NULL)}
+}
+\if{html}{\out{
}}
+}
+\subsection{Returns}{
+A \code{bd} object containing the optimization results
+}
+}
+\if{html}{\out{
}}
+\if{html}{\out{}}
+\if{latex}{\out{\hypertarget{method-bd_post-clone}{}}}
+\subsection{Method \code{clone()}}{
+The objects of this class are cloneable with this method.
+\subsection{Usage}{
+\if{html}{\out{}}\preformatted{bd_post$clone(deep = FALSE)}\if{html}{\out{
}}
+}
+
+\subsection{Arguments}{
+\if{html}{\out{}}
+\describe{
+\item{\code{deep}}{Whether to make a deep clone.}
+}
+\if{html}{\out{
}}
+}
+}
+}
diff --git a/man/compile_model.Rd b/man/compile_model.Rd
index 93d82703..d0727715 100644
--- a/man/compile_model.Rd
+++ b/man/compile_model.Rd
@@ -23,4 +23,3 @@ Result of the compilation process.
\description{
Compiles the text representation of a model into a binary file.
}
-\keyword{internal}
diff --git a/man/dummy_compile.Rd b/man/dummy_compile.Rd
index f28ad59a..cbaf7edc 100644
--- a/man/dummy_compile.Rd
+++ b/man/dummy_compile.Rd
@@ -15,4 +15,3 @@ Path to the build directory.
\description{
Dummy function to cache compilation artifacts.
}
-\keyword{internal}
diff --git a/man/fit.Rd b/man/fit.Rd
index 287a8ef2..443b982f 100644
--- a/man/fit.Rd
+++ b/man/fit.Rd
@@ -23,4 +23,3 @@ Result of the fitting process.
\description{
Fits the model at the given path to the data at the given path using the provided parameters.
}
-\keyword{internal}
diff --git a/man/is_cargo_installed.Rd b/man/is_cargo_installed.Rd
index ae0fb2b1..9367a427 100644
--- a/man/is_cargo_installed.Rd
+++ b/man/is_cargo_installed.Rd
@@ -12,4 +12,3 @@ TRUE if Cargo is installed, FALSE otherwise.
\description{
Checks if Cargo is installed on the system.
}
-\keyword{internal}
diff --git a/man/model_parameters.Rd b/man/model_parameters.Rd
index 414e4ba8..46475849 100644
--- a/man/model_parameters.Rd
+++ b/man/model_parameters.Rd
@@ -17,4 +17,3 @@ List of model parameters.
\description{
Retrieves the model parameters from the compiled model at the given path.
}
-\keyword{internal}
diff --git a/man/plot.bd.Rd b/man/plot.bd.Rd
new file mode 100644
index 00000000..56cfb65a
--- /dev/null
+++ b/man/plot.bd.Rd
@@ -0,0 +1,63 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PM_bestdose.R
+\name{plot.bd}
+\alias{plot.bd}
+\title{Plot BestDose predictions}
+\usage{
+\method{plot}{bd}(
+ x,
+ mult = 1,
+ outeq = 1,
+ quiet = FALSE,
+ legend = TRUE,
+ log = FALSE,
+ grid = FALSE,
+ future_region = FALSE,
+ xlab = "Time",
+ ylab = "Concentration",
+ title = NULL,
+ xlim = NULL,
+ ylim = NULL,
+ print = TRUE,
+ ...
+)
+}
+\arguments{
+\item{x}{A \code{bd} object containing the best dose predictions.}
+
+\item{mult}{Numeric multiplier to apply to the concentrations (default: 1, no scaling).}
+
+\item{outeq}{Numeric value of outeq to filter observations for plotting (default: 1).}
+
+\item{quiet}{Logical indicating whether to suppress messages (default: FALSE).}
+
+\item{legend}{Logical indicating whether to display a legend (default: TRUE).}
+
+\item{log}{Logical indicating whether to use a logarithmic scale for the y-axis (default: FALSE).}
+
+\item{grid}{Logical indicating whether to display a grid (default: FALSE).}
+
+\item{future_region}{Logical indicating whether to display the divider line and shaded region for future data (default: FALSE).}
+
+\item{xlab}{Label for the x-axis.}
+
+\item{ylab}{Label for the y-axis.}
+
+\item{title}{Title of the plot.}
+
+\item{xlim}{Limits for the x-axis.}
+
+\item{ylim}{Limits for the y-axis.}
+
+\item{print}{Logical indicating whether to print the plot (default: TRUE).}
+
+\item{...}{Additional arguments passed to the plotting function.}
+}
+\value{
+Invisibly returns a plotly object.
+}
+\description{
+Plot observed and predicted concentrations over time for both past and future data,
+with options to apply a multiplier to the concentrations, and customize the plot appearance. The top 5 posterior support
+points are highlighted with distinct colors.
+}
diff --git a/man/simulate_all.Rd b/man/simulate_all.Rd
index 2a6378b9..e7ebda22 100644
--- a/man/simulate_all.Rd
+++ b/man/simulate_all.Rd
@@ -21,4 +21,3 @@ Simulation results.
\description{
Simulates all subjects in the data set using the model at the given path.
}
-\keyword{internal}
diff --git a/man/simulate_one.Rd b/man/simulate_one.Rd
index 3cd97646..80c95c9f 100644
--- a/man/simulate_one.Rd
+++ b/man/simulate_one.Rd
@@ -21,4 +21,3 @@ Simulation results.
\description{
Simulates the first subject in the data set using the model at the given path.
}
-\keyword{internal}
diff --git a/man/temporary_path.Rd b/man/temporary_path.Rd
index f3437de8..ba55af0f 100644
--- a/man/temporary_path.Rd
+++ b/man/temporary_path.Rd
@@ -12,4 +12,3 @@ Temporary build path.
\description{
Retrieves the temporary path used for building models.
}
-\keyword{internal}
diff --git a/output.txt b/output.txt
new file mode 100644
index 00000000..82ceca1f
--- /dev/null
+++ b/output.txt
@@ -0,0 +1,8 @@
+ M Cargo.lock
+ M src/rust/Cargo.toml
+?? inst/Examples/Rscript/sylvain.R
+?? inst/Examples/Rscript/sylvain_sim.R
+?? inst/Examples/Runs/
+?? inst/Examples/src/Dapto_Lyon.csv
+?? output.txt
+{"additions":202,"author":{"id":"MDQ6VXNlcjEwMjMwMDY=","is_bot":false,"login":"Siel","name":"Julian Otalvaro"},"baseRefName":"main","body":"","changedFiles":6,"commits":[{"authoredDate":"2026-04-17T12:03:14Z","authors":[{"email":"juliandavid347@gmail.com","id":"MDQ6VXNlcjEwMjMwMDY=","login":"Siel","name":"Julián D. Otálvaro"}],"committedDate":"2026-04-17T12:03:14Z","messageBody":"","messageHeadline":"Fix analytical public indexing compatibility","oid":"6183707da86fc76e3e3ed461dd4eae0e3e4a8c17"},{"authoredDate":"2026-04-17T12:22:02Z","authors":[{"email":"juliandavid347@gmail.com","id":"MDQ6VXNlcjEwMjMwMDY=","login":"Siel","name":"Julián D. Otálvaro"}],"committedDate":"2026-04-17T12:22:02Z","messageBody":"","messageHeadline":"cleanup","oid":"7395ce4a34ed8f4c2ca39b9755ea4d0f2e9e6f40"},{"authoredDate":"2026-04-20T17:34:06Z","authors":[{"email":"juliandavid347@gmail.com","id":"MDQ6VXNlcjEwMjMwMDY=","login":"Siel","name":"Julián D. Otálvaro"}],"committedDate":"2026-04-20T17:34:06Z","messageBody":"","messageHeadline":"chore: update PMcore to 0.25.2","oid":"d05c575ac7efd0c4d5b0e78c3f281b2fe9657672"},{"authoredDate":"2026-04-20T17:41:13Z","authors":[{"email":"juliandavid347@gmail.com","id":"MDQ6VXNlcjEwMjMwMDY=","login":"Siel","name":"Julián D. Otálvaro"}],"committedDate":"2026-04-20T17:41:13Z","messageBody":"","messageHeadline":"chore: pin extendr-api to 0.9.0","oid":"9fcc5c9e3fae635c8812184737b2e8768028ff11"},{"authoredDate":"2026-04-20T17:49:52Z","authors":[{"email":"markushh@uio.no","id":"MDQ6VXNlcjY2MDU4NjQy","login":"mhovd","name":"Markus"}],"committedDate":"2026-04-20T17:49:52Z","messageBody":"","messageHeadline":"chore: Update Cargo.toml","oid":"bea323a88bbe3945791637db66132ae911e0705a"}],"deletions":165,"files":[{"path":"Cargo.lock","additions":171,"deletions":133},{"path":"R/PM_model.R","additions":2,"deletions":2},{"path":"R/model_transpiler.R","additions":24,"deletions":20},{"path":"src/rust/Cargo.toml","additions":3,"deletions":5},{"path":"src/rust/src/lib.rs","additions":1,"deletions":4},{"path":"tests/testthat/test-analytical-indexing.R","additions":1,"deletions":1}],"headRefName":"fix/pmetrics-analytical","mergeCommit":null,"title":"fix: Analytical public indexing compatibility"}
diff --git a/src/rust/Cargo.toml b/src/rust/Cargo.toml
index c86e0cbf..bc263609 100755
--- a/src/rust/Cargo.toml
+++ b/src/rust/Cargo.toml
@@ -9,7 +9,8 @@ name = 'pm_rs'
[dependencies]
extendr-api = "=0.9.0"
-pmcore = { version = "=0.25.2", features = ["exa"] }
+# pmcore = { version = "=0.25.2", features = ["exa"] }
+pmcore = { git = "https://github.com/LAPKB/PMcore.git", branch="fix/pmetrics-analytical", features = ["exa"] }
libloading = "0.9"
rayon = "1.10.0"
diff --git a/src/rust/src/bestdose_executor.rs b/src/rust/src/bestdose_executor.rs
new file mode 100644
index 00000000..4fbf8ad6
--- /dev/null
+++ b/src/rust/src/bestdose_executor.rs
@@ -0,0 +1,441 @@
+use crate::{logs::RFormatLayer, settings::settings};
+use extendr_api::prelude::*;
+use pmcore::bestdose::{BestDosePosterior, BestDoseResult, DoseRange, Target};
+use pmcore::prelude::{data, ODE};
+use pmcore::routines::initialization::parse_prior;
+use std::{
+ any::Any,
+ path::PathBuf,
+ sync::{Arc, Mutex},
+};
+
+/// Helper to parse target type from string
+pub(crate) fn parse_target_type(target_str: &str) -> std::result::Result {
+ match target_str.to_lowercase().as_str() {
+ "concentration" => Ok(Target::Concentration),
+ "auc_from_zero" | "auc" => Ok(Target::AUCFromZero),
+ "auc_from_last_dose" | "auc_interval" => Ok(Target::AUCFromLastDose),
+ _ => Err(format!(
+ "Invalid target type: {}. Must be 'concentration', 'auc_from_zero', or 'auc_from_last_dose'",
+ target_str
+ )),
+ }
+}
+
+/// R-compatible prediction row for BestDose output
+#[derive(Debug, IntoDataFrameRow)]
+pub struct BestDosePredictionRow {
+ id: String,
+ time: f64,
+ observed: f64,
+ pop_mean: f64,
+ pop_median: f64,
+ post_mean: f64,
+ post_median: f64,
+ outeq: usize,
+}
+
+impl BestDosePredictionRow {
+ 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(),
+ }
+ }
+}
+
+/// R-compatible AUC prediction row
+#[derive(Debug, IntoDataFrameRow)]
+pub struct BestDoseAucRow {
+ time: f64,
+ auc: f64,
+}
+
+/// Convert BestDoseResult to R-compatible list structure
+pub(crate) fn convert_bestdose_result_to_r(
+ result: BestDoseResult,
+) -> std::result::Result {
+ // Extract doses
+ let doses: Vec = result.doses();
+
+ // Objective function
+ let objf = result.objf();
+
+ // Status
+ let status_str = format!("{:?}", result.status());
+
+ // Predictions as data frame
+ let pred_rows: Vec = result
+ .predictions()
+ .predictions()
+ .iter()
+ .map(|p| BestDosePredictionRow::from_np_prediction(p, "subject_1"))
+ .collect();
+ let pred_df = pred_rows
+ .into_dataframe()
+ .map_err(|e| format!("Failed to create predictions dataframe: {:?}", e))?;
+
+ // AUC predictions (if available)
+ let auc_val = if let Some(auc_preds) = result.auc_predictions() {
+ let auc_rows: Vec = auc_preds
+ .iter()
+ .map(|(time, auc)| BestDoseAucRow {
+ time: *time,
+ auc: *auc,
+ })
+ .collect();
+ let auc_df = auc_rows
+ .into_dataframe()
+ .map_err(|e| format!("Failed to create AUC dataframe: {:?}", e))?;
+ Robj::from(auc_df)
+ } else {
+ Robj::from(()) // NULL for no AUC
+ };
+
+ // Optimization method
+ let method_str = format!("{}", result.optimization_method());
+
+ // Build the list using list! macro
+ let output = list!(
+ doses = doses,
+ objf = objf,
+ status = status_str,
+ predictions = pred_df,
+ auc_predictions = auc_val,
+ method = method_str
+ );
+
+ Ok(output.into())
+}
+/// Opaque handle that keeps the dynamic model library alive while reusing the
+/// prepared `BestDosePosterior` for multiple optimization runs.
+pub struct BestDosePosteriorHandle {
+ posterior: BestDosePosterior,
+ #[allow(dead_code)]
+ library: libloading::Library,
+}
+
+impl BestDosePosteriorHandle {
+ pub fn new(
+ model_path: PathBuf,
+ prior_path: PathBuf,
+ past_data_path: Option,
+ params: List,
+ ) -> std::result::Result {
+ let (library, (eq, meta)) =
+ unsafe { pmcore::prelude::pharmsol::exa::load::load::(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) =
+ parse_prior(&prior_path.to_str().unwrap().to_string(), &settings)
+ .map_err(|e| format!("Failed to parse prior: {}", e))?;
+
+ let population_weights = prior_weights
+ .ok_or_else(|| "Prior file must contain a 'prob' column with weights".to_string())?;
+
+ let past_data = if let Some(path) = past_data_path {
+ let data = data::read_pmetrics(path.to_str().unwrap())
+ .map_err(|e| format!("Failed to read past data: {}", e))?;
+ let subjects = data.subjects();
+ if subjects.is_empty() {
+ return Err("Past data file contains no subjects".to_string());
+ }
+ Some(subjects[0].clone())
+ } else {
+ None
+ };
+
+ let posterior = BestDosePosterior::compute(
+ &population_theta,
+ &population_weights,
+ past_data,
+ eq,
+ settings,
+ )
+ .map_err(|e| format!("Failed to compute BestDose posterior: {}", e))?;
+
+ Ok(Self { posterior, library })
+ }
+
+ #[allow(clippy::too_many_arguments)]
+ pub fn optimize(
+ &self,
+ target_data_path: PathBuf,
+ time_offset: Option,
+ dose_min: f64,
+ dose_max: f64,
+ bias_weight: f64,
+ target_type: &str,
+ ) -> std::result::Result {
+ let target_data = {
+ let data = data::read_pmetrics(target_data_path.to_str().unwrap())
+ .map_err(|e| format!("Failed to read target data: {}", e))?;
+ let subjects = data.subjects();
+ if subjects.is_empty() {
+ return Err("Target data file contains no subjects".to_string());
+ }
+ subjects[0].clone()
+ };
+
+ let target_enum = parse_target_type(target_type)?;
+ let dose_range = DoseRange::new(dose_min, dose_max);
+
+ self.posterior
+ .optimize(
+ target_data,
+ time_offset,
+ dose_range,
+ bias_weight,
+ target_enum,
+ )
+ .map_err(|e| format!("Optimization failed: {}", e))
+ }
+
+ pub fn posterior(&self) -> &BestDosePosterior {
+ &self.posterior
+ }
+}
+
+pub(crate) struct PosteriorSummary {
+ theta_values: Vec,
+ theta_dim: (i32, i32),
+ param_names: Vec,
+ posterior_weights: Vec,
+ population_weights: Vec,
+}
+
+fn summarize_handle(handle: &BestDosePosteriorHandle) -> PosteriorSummary {
+ let theta = handle.posterior().theta();
+ let matrix = theta.matrix();
+ let nrows = matrix.nrows() as i32;
+ let ncols = matrix.ncols() as i32;
+ let mut theta_values = vec![0.0; (nrows * ncols) as usize];
+
+ for col in 0..ncols as usize {
+ for row in 0..nrows as usize {
+ theta_values[row + col * nrows as usize] = *matrix.get(row, col);
+ }
+ }
+
+ PosteriorSummary {
+ theta_values,
+ theta_dim: (nrows, ncols),
+ param_names: theta.param_names(),
+ posterior_weights: handle.posterior().posterior_weights().to_vec(),
+ population_weights: handle.posterior().population_weights().to_vec(),
+ }
+}
+
+fn vec_to_doubles(values: Vec, label: &str) -> std::result::Result {
+ Doubles::try_from(values)
+ .map_err(|e| format!("Failed to convert {} to doubles: {:?}", label, e))
+}
+
+fn dims_to_integers(dim: (i32, i32)) -> std::result::Result {
+ Integers::try_from(vec![dim.0, dim.1])
+ .map_err(|e| format!("Failed to convert dims to integers: {:?}", e))
+}
+
+fn names_to_strings(names: &[String]) -> Strings {
+ Strings::from_values(names.iter().map(|s| s.as_str()))
+}
+
+fn panic_payload_to_string(payload: &(dyn Any + Send)) -> String {
+ if let Some(message) = payload.downcast_ref::<&str>() {
+ (*message).to_string()
+ } else if let Some(message) = payload.downcast_ref::() {
+ message.clone()
+ } else {
+ "panic with non-string payload".to_string()
+ }
+}
+
+fn panic_info_to_string(info: &std::panic::PanicHookInfo<'_>) -> String {
+ let location = info
+ .location()
+ .map(|location| format!("{}:{}", location.file(), location.line()))
+ .unwrap_or_else(|| "unknown location".to_string());
+
+ format!(
+ "{} at {}",
+ panic_payload_to_string(info.payload()),
+ location
+ )
+}
+
+fn run_with_panic_hook(context: &'static str, operation: F) -> std::result::Result
+where
+ F: FnOnce() -> std::result::Result,
+{
+ let panic_message = Arc::new(Mutex::new(None::));
+ let panic_message_for_hook = Arc::clone(&panic_message);
+ let previous_hook = std::panic::take_hook();
+
+ std::panic::set_hook(Box::new(move |info| {
+ let formatted = panic_info_to_string(info);
+ if let Ok(mut slot) = panic_message_for_hook.lock() {
+ *slot = Some(formatted.clone());
+ }
+ tracing::error!(context = context, panic = %formatted, "panic intercepted");
+ }));
+
+ let result = std::panic::catch_unwind(std::panic::AssertUnwindSafe(operation));
+ std::panic::set_hook(previous_hook);
+
+ match result {
+ Ok(result) => result,
+ Err(payload) => {
+ let panic_message = panic_message
+ .lock()
+ .ok()
+ .and_then(|slot| slot.clone())
+ .unwrap_or_else(|| panic_payload_to_string(payload.as_ref()));
+ Err(format!("{} panicked: {}", context, panic_message))
+ }
+ }
+}
+
+pub(crate) fn prepare_bestdose_posterior(
+ model_path: PathBuf,
+ prior_path: PathBuf,
+ past_data_path: Option,
+ params: List,
+) -> std::result::Result<(BestDosePosteriorHandle, PosteriorSummary), String> {
+ let handle = BestDosePosteriorHandle::new(model_path, prior_path, past_data_path, params)?;
+
+ let summary = summarize_handle(&handle);
+ Ok((handle, summary))
+}
+
+pub(crate) fn bestdose_prepare_internal(
+ model_path: &str,
+ prior_path: &str,
+ past_data_path: Nullable,
+ params: List,
+ kind: &str,
+) -> Robj {
+ RFormatLayer::reset_global_timer();
+ let _ = crate::setup_logs();
+
+ let past_path = past_data_path.into_option().map(PathBuf::from);
+
+ let preparation = match kind {
+ "ode" => prepare_bestdose_posterior(
+ PathBuf::from(model_path),
+ PathBuf::from(prior_path),
+ past_path,
+ params.clone(),
+ ),
+ "analytical" => Err("BestDose for analytical models is not yet supported".to_string()),
+ other => Err(format!("{} is not a supported model type", other)),
+ };
+
+ match preparation {
+ Ok((handle, summary)) => {
+ let theta_values = match vec_to_doubles(summary.theta_values, "theta_values") {
+ Ok(values) => values,
+ Err(e) => return Robj::from(e),
+ };
+ let theta_dim = match dims_to_integers(summary.theta_dim) {
+ Ok(dim) => dim,
+ Err(e) => return Robj::from(e),
+ };
+ let posterior_weights =
+ match vec_to_doubles(summary.posterior_weights, "posterior_weights") {
+ Ok(values) => values,
+ Err(e) => return Robj::from(e),
+ };
+ let population_weights =
+ match vec_to_doubles(summary.population_weights, "population_weights") {
+ Ok(values) => values,
+ Err(e) => return Robj::from(e),
+ };
+ let param_names = names_to_strings(&summary.param_names);
+ let handle_ptr = ExternalPtr::new(handle);
+
+ let output = list!(
+ handle = handle_ptr,
+ theta_values = theta_values,
+ theta_dim = theta_dim,
+ param_names = param_names,
+ posterior_weights = posterior_weights,
+ population_weights = population_weights,
+ nspp = summary.theta_dim.0,
+ n_parameters = summary.theta_dim.1
+ );
+
+ output.into()
+ }
+ Err(e) => Robj::from(format!("BestDose prepare failed: {}", e)),
+ }
+}
+
+pub(crate) fn bestdose_optimize_internal(
+ handle: ExternalPtr,
+ target_data_path: &str,
+ time_offset: Nullable,
+ dose_min: f64,
+ dose_max: f64,
+ bias_weight: f64,
+ target_type: &str,
+) -> Robj {
+ let time_offset = time_offset.into_option();
+
+ tracing::info!(
+ target_data_path = target_data_path,
+ time_offset = ?time_offset,
+ dose_min = dose_min,
+ dose_max = dose_max,
+ bias_weight = bias_weight,
+ target_type = target_type,
+ "bestdose_optimize_internal parameters"
+ );
+
+ match run_with_panic_hook("bestdose_optimize_internal", || match handle.try_addr() {
+ Ok(inner) => {
+ tracing::info!(
+ handle = %format!("{:p}", inner as *const BestDosePosteriorHandle),
+ target_data_path = target_data_path,
+ time_offset = ?time_offset,
+ dose_min = dose_min,
+ dose_max = dose_max,
+ bias_weight = bias_weight,
+ target_type = target_type,
+ "bestdose_optimize_internal invoking optimize"
+ );
+
+ match inner.optimize(
+ PathBuf::from(target_data_path),
+ time_offset,
+ dose_min,
+ dose_max,
+ bias_weight,
+ target_type,
+ ) {
+ Ok(result) => match convert_bestdose_result_to_r(result) {
+ Ok(robj) => Ok(robj),
+ Err(e) => Ok(Robj::from(format!("Failed to convert result: {}", e))),
+ },
+ Err(e) => Ok(Robj::from(format!("BestDose optimization failed: {}", e))),
+ }
+ }
+ Err(e) => {
+ tracing::error!(handle_error = %e, "bestdose_optimize_internal received an invalid handle");
+ Ok(Robj::from(format!("Invalid BestDose handle: {}", e)))
+ }
+ }) {
+ Ok(robj) => robj,
+ Err(e) => Robj::from(format!("BestDose optimization panicked: {}", e)),
+ }
+}
diff --git a/src/rust/src/lib.rs b/src/rust/src/lib.rs
index a2e19e85..75a77c62 100755
--- a/src/rust/src/lib.rs
+++ b/src/rust/src/lib.rs
@@ -1,5 +1,5 @@
// mod build;
-
+mod bestdose_executor;
mod executor;
mod logs;
mod settings;
@@ -16,6 +16,7 @@ use simulation::SimulationRow;
use std::process::Command;
use tracing_subscriber::layer::SubscriberExt;
+use crate::bestdose_executor::BestDosePosteriorHandle;
use crate::logs::RFormatLayer;
fn validate_paths(data_path: &str, model_path: &str) -> Result<()> {
@@ -266,7 +267,7 @@ fn model_parameters(model_path: &str, kind: &str) -> Result> {
/// @return Temporary build path.
/// @export
#[extendr]
-fn temporary_path() -> String {
+pub fn temporary_path() -> String {
build::temp_path().to_string_lossy().to_string()
}
@@ -291,6 +292,47 @@ fn setup_logs() -> anyhow::Result<()> {
Ok(())
}
+// Macro to generate exports.
+// This ensures exported functions are registered with R.
+// See corresponding C code in `entrypoint.c`.
+#[extendr]
+fn bestdose_prepare(
+ model_path: &str,
+ prior_path: &str,
+ past_data_path: Nullable,
+ params: List,
+ kind: &str,
+) -> Robj {
+ bestdose_executor::bestdose_prepare_internal(
+ model_path,
+ prior_path,
+ past_data_path,
+ params,
+ kind,
+ )
+}
+
+#[extendr]
+fn bestdose_optimize(
+ handle: ExternalPtr,
+ target_data_path: &str,
+ time_offset: Nullable,
+ dose_min: f64,
+ dose_max: f64,
+ bias_weight: f64,
+ target_type: &str,
+) -> Robj {
+ bestdose_executor::bestdose_optimize_internal(
+ handle,
+ target_data_path,
+ time_offset,
+ dose_min,
+ dose_max,
+ bias_weight,
+ target_type,
+ )
+}
+
extendr_module! {
mod Pmetrics;
fn simulate_one;
@@ -302,7 +344,8 @@ extendr_module! {
fn model_parameters;
fn temporary_path;
fn setup_logs;
-
+ fn bestdose_prepare;
+ fn bestdose_optimize;
}
// To generate the exported function in R, run the following command: