diff --git a/Cargo.lock b/Cargo.lock index 8787dd928..bd0cced3a 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1435,9 +1435,9 @@ dependencies = [ [[package]] name = "pharmsol" -version = "0.29.0" +version = "0.29.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d956a0a3ae74b4a02c16cd3dff3153a7b958f3f8c48b397f469df6151c1fc339" +checksum = "99dde6038703fc17c65bee8758f889ded76cc1d5c75b101ab399cb35c82e3854" dependencies = [ "ahash", "argmin", @@ -1466,9 +1466,9 @@ dependencies = [ [[package]] name = "pharmsol-dsl" -version = "0.29.0" +version = "0.29.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "80ff36016934559d648d3d3dcdf99003c2f314ff7a00e97d8bc471c609e0db0f" +checksum = "afaa658ecea67c8d44784bd8e4447da8c7a218699fa30e52f907d423b8b622ca" dependencies = [ "serde", "serde_json", @@ -1477,9 +1477,9 @@ dependencies = [ [[package]] name = "pharmsol-macros" -version = "0.29.0" +version = "0.29.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "44344f94af37484c5f09f75af566b23baf9643f06076aa27ef741f5e657b47ba" +checksum = "3673a31de5a0b739892d6739a3fc2639340667c7a19a8549db122f508022337b" dependencies = [ "pharmsol-dsl", "proc-macro-crate", @@ -1508,9 +1508,9 @@ dependencies = [ [[package]] name = "pmcore" -version = "0.28.0" +version = "0.29.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "732c974a86e0970239056ed726e0a7269856a957d3095c024461813cd65bd4da" +checksum = "941b66ea8169eb15284e24ac63b51da18b81c5abb501d396ce28b0750a49b390" dependencies = [ "anyhow", "argmin", diff --git a/R/PM_model.R b/R/PM_model.R index 4a0ecac25..0ad33f20c 100644 --- a/R/PM_model.R +++ b/R/PM_model.R @@ -252,19 +252,30 @@ PM_model <- R6::R6Class( #' - `dx[i]` for the change in amount with respect to time (i.e., \eqn{dx/dt}), #' where `i` is the compartment number, #' - `x[i]` for the compartment amount, where `i` is the compartment number. - #' - `rateiv[j]` for the infusion rate of input `j`, where `j` is the input number - #' in the data corresponding to doses for that input. - #' - Bolus doses are indicated by `DUR = 0` for dose events in the - #' data. Currently only one bolus input is allowed, which goes into compartment 1 - #' and is not modifiable. It does not appear in the differential equations. + #' - `rateiv[j]` or `R[j]` for the infusion rate of input `j`, where `j` is + #' the input number in the data. An infusion may be a standalone additive term + #' or one exact product with a scale, such as `R[1] * iv_scale`. The scale is + #' evaluated continuously in the differential equation. + #' - `B[j]` or `bolus[j]` for a bolus input. A bolus may be a standalone + #' additive term or one exact product with a scale, such as `B[1] * v`. + #' The scale is applied to the bolus dose and does not affect an infusion with + #' the same input number. #' #' For example, #' ``` #' eqn = function() { - #' dx[1] = -ka * x[1] - #' dx[2] = rateiv[1] + ka * x[1] - ke * x[2] + #' dx[1] = B[1] * fa1 - ka * x[1] + #' dx[2] = rateiv[1] * iv_scale + ka * x[1] - ke * x[2] #' } #' ``` + #' * **Expressions** are ordinary R arithmetic. Conditions may appear anywhere, + #' not only as a whole equation: `2 * if (wt > 70) v1 else v2` and + #' `if (a > 0) (if (b > 0) 1 else 2) else 3` both work. `ifelse(cond, a, b)`, + #' `pmin`, `pmax`, and `max`/`min` with any number of arguments are accepted; + #' `log(x, base)`, `trunc`, `sign` and `round(x, digits)` follow R's rules + #' (including R's round-half-to-even). Functions R provides but the DSL does + #' not - such as `pnorm`, `log1p` or `gamma` - are reported when the model is + #' defined, naming the call and the supported alternatives. #' * **Additional equations** in R code can be defined in this block, which are similar to #' the `sec` block, but will only be available within the `eqn` block as opposed #' to global availability when defined in `sec`. They can be added to either @@ -286,7 +297,24 @@ PM_model <- R6::R6Class( #' As for `eqn`, additional equations in R code can be defined in this block, #' but will only be available within the `lag` block. #' @param fa A function defining the bioavailability (fraction absorbed) equations, - #' similar to `lag`. + #' similar to `lag`. The `lag` and `fa` blocks also apply to analytical bolus + #' models; they are not valid for analytical infusion-only models. + #' + #' Effect functions for combination drug action are also available, anywhere an + #' expression is allowed. They compute an effect from supplied coefficients and + #' do not fit those coefficients: + #' + #' - `estimate_effect_2(u, v, alpha, h1, h2)` for two sites, and + #' - `estimate_effect_3(a, b, c, alpha12, alpha13, alpha23, alpha123, h1, h2, h3)` + #' for three. + #' + #' `u` and `v` are normalised exposures (concentration divided by the site's + #' `e50`), `alpha` is the interaction coefficient and `h1`/`h2` are the Hill + #' exponents. The interaction term is computed internally as + #' \eqn{w = \alpha \cdot u \cdot v}; the retired six-argument `get_e2(a, b, w, h1, + #' h2, alpha_s)` call must be rewritten as `estimate_effect_2(a, b, alpha_s, h1, + #' h2)`. These functions are evaluated at run time, on every derivative or + #' output evaluation, and cannot be used in the `pri` block. #' #' Example: #' ``` @@ -308,7 +336,8 @@ PM_model <- R6::R6Class( #' This sets the initial amount of drug in compartment 2 to the value #' of a covariate `init2` multiplied by the volume of the compartment, #' `V`, assuming `V` is either a primary parameter or defined in the - #' `sec` block. + #' `sec` block. Initial conditions are evaluated before the first event of + #' every occasion, including occasions started by an `EVID = 4` row. #' #' As for `eqn`, additional equations in R code can be defined in this block, #' but will only be available within the `ini` block. @@ -471,7 +500,10 @@ PM_model <- R6::R6Class( } }) self$arg_list$x <- NULL - self$dsl <- x$dsl + # `dsl` is rendered from `arg_list`, so it is never inherited. A + # PM_model saved before a DSL change would otherwise replay stale + # source instead of rendering the model it defines. + self$dsl <- NULL } else { cli::cli_abort(c( "x" = "Non supported input for {.arg x}: {typeof(x)}", @@ -1447,20 +1479,10 @@ PM_model <- R6::R6Class( fit_call() } }, - error = function(e) { - msg_txt <- conditionMessage(e) - # The DSL is generated from the model definition, so the backend's - # line numbers refer to generated text, not the user's R script. - if (grepl("compile model|parse DSL|DSL[0-9]{4}", msg_txt, ignore.case = TRUE)) { - hint <- paste( - "The model source is generated from your model definition;", - "any reported line numbers refer to that generated text, not your R script.", - "Check the secondary, output, and derivative equations for unsupported constructs." - ) - stop(paste0(msg_txt, "\n\n", hint), call. = FALSE) - } - stop(e) - } + # The model was already validated when it was compiled, so any DSL + # diagnostic here is reported as it comes from the backend, with the + # R statement that produced it already named. + error = function(e) stop(e) ) # The Rust backend writes the estimation artifacts (theta.csv, @@ -1653,16 +1675,18 @@ PM_model <- R6::R6Class( } if (!quiet) cli::cli_inform(c("i" = "Preparing model...")) - tryCatch( - { - self$dsl <- model_to_dsl(self) - }, + rendered <- tryCatch( + model_to_dsl_traced(self), error = function(e) { cli::cli_abort( c("x" = "Model preparation failed: {conditionMessage(e)}", "i" = "Please check the model definition and try again.") ) } ) + self$dsl <- rendered$text + # Validate as soon as the model is rendered, so a model that cannot be + # compiled fails where it is defined instead of at the start of a run. + private$validate_dsl(rendered) return(invisible(self)) }, # end compile method @@ -1849,6 +1873,18 @@ PM_model <- R6::R6Class( } # end copy ), # end public list private = list( + # Ask the backend to parse and analyze the rendered source. Diagnostics are + # reported against the R statement that produced the offending line. + validate_dsl = function(rendered) { + tryCatch( + model_metadata(rendered$text, self$model_list$solver), + error = function(e) { + cli::cli_abort(dsl_error_bullets(conditionMessage(e), rendered)) + } + ) + invisible(NULL) + }, + from_dsl = function(dsl, pri, err, solver, quiet) { if (!is.character(dsl) || length(dsl) != 1 || is.na(dsl) || !nzchar(dsl)) { cli::cli_abort("{.arg dsl} must be one non-empty character string.") diff --git a/R/model_dsl.R b/R/model_dsl.R index 13dd509ee..b8e026f9c 100644 --- a/R/model_dsl.R +++ b/R/model_dsl.R @@ -17,8 +17,8 @@ # Naming conventions used by the emitter: # * States `x[i]` -> `x{i}` (declared in `states = ...`) # * Outputs `Y[i]` / `y[i]` -> `outeq_{i}` (declared in `outputs = ...`) -# * Inputs `b[j]`/`bolus[j]` -> `bolus(input_{j}) -> x{k}` (route) -# `rateiv[j]`/`r[j]`-> `infusion(input_{j}) -> x{k}` (route) +# * ODE inputs `b[j]`/`bolus[j]` -> `bolus(input_{j})` (RHS term) +# `rateiv[j]`/`r[j]`-> `infusion(input_{j})` (RHS term) # * Params / covariates keep their (lower-cased) names. # # pharmsol identifies routes and outputs by label and rejects bare numeric @@ -27,10 +27,9 @@ # renumbered: the data's `INPUT`/`OUTEQ` values are canonicalised to the same # labels before the backend reads them. # -# Infusions and boluses are declared as routes into the compartment in which they -# appear in the derivative equations; the corresponding `rateiv[j]`/`b[j]` terms -# are stripped from the derivative because the DSL runtime injects them -# automatically from the route declaration. +# ODE inputs and scales remain on the derivative RHS. The DSL infers their routes +# and lowers them into the existing simulator machinery. Analytical models keep +# explicit route declarations and route properties. # --------------------------------------------------------------------------- # The set of DSL intrinsic functions and how R names map to them. R names that @@ -46,10 +45,315 @@ dsl_function_map <- function() { ) } +# Runtime-only effect functions provided by pharmsol, with their required +# argument counts. Pmetrics validates the arity here so a wrong call is +# reported against the R expression rather than against generated DSL text. +dsl_effect_functions <- function() { + c(estimate_effect_2 = 5L, estimate_effect_3 = 10L) +} + +# The documented signature of each runtime effect function, used in messages. +dsl_effect_signature <- function() { + c( + estimate_effect_2 = "estimate_effect_2(u, v, alpha, h1, h2)", + estimate_effect_3 = "estimate_effect_3(a, b, c, alpha12, alpha13, alpha23, alpha123, h1, h2, h3)" + ) +} + +# Function names that used to exist with different arguments. They get a +# migration message instead of a generic "unsupported function" error, because +# the replacement changes the *meaning* of the call rather than just its name. +dsl_legacy_functions <- function() { + c( + get_e2 = paste( + "Replace {.code get_e2(a, b, w, h1, h2, alpha_s)} with", + "{.code estimate_effect_2(a, b, alpha_s, h1, h2)}." + ), + get_e3 = paste( + "Replace the six-argument {.code get_e3} call with", + "{.code estimate_effect_3(a, b, c, alpha12, alpha13, alpha23, alpha123, h1, h2, h3)}." + ) + ) +} + +# The name of the function a call invokes, or NULL when the head is not a plain +# symbol. `base::exp(x)` is a call whose *head* is another call, and +# `as.character()` on that head returns several elements, so comparing it to a +# string raised a base R error instead of a model diagnostic. +dsl_call_head <- function(expr) { + if (!is.call(expr)) { + return(NULL) + } + head <- expr[[1]] + if (is.symbol(head)) as.character(head) else NULL +} + +# One-line rendering of an R expression, used to quote the user's own code in +# diagnostics instead of generated DSL text. +dsl_expr_label <- function(expr) { + paste(deparse(expr, width.cutoff = 500L), collapse = " ") +} + +# The numeric value of a literal, or NULL when the expression is not one. A +# negative literal is a call to `-` in R's parse tree, so `round(x, -1)` is +# recognised here rather than being rejected as a non-literal. +dsl_literal_number <- function(expr) { + if (is.numeric(expr) && length(expr) == 1 && is.finite(expr)) { + return(as.numeric(expr)) + } + if (is.call(expr) && length(expr) == 2L && identical(expr[[1]], as.name("-")) && + is.numeric(expr[[2]]) && length(expr[[2]]) == 1 && is.finite(expr[[2]])) { + return(-as.numeric(expr[[2]])) + } + NULL +} + +# Abort with a message that names the unsupported call and what to write +# instead. Every rejection in this file goes through this helper so the wording +# stays consistent and always carries a fix. `reason` and `fix` are plain text: +# braces in them are escaped so cli shows them literally. +dsl_literal <- function(x) { + gsub("}", "}}", gsub("{", "{{", x, fixed = TRUE), fixed = TRUE) +} + +dsl_unsupported <- function(expr, reason, fix = NULL) { + bullets <- c("x" = paste0(dsl_literal(reason), " {.code {dsl_expr_label(expr)}}")) + if (!is.null(fix)) bullets <- c(bullets, "i" = dsl_literal(fix)) + cli::cli_abort(bullets, call = NULL) +} + +# Turn a backend DSL diagnostic into one that points at the user's model. The +# backend reports positions in generated text; `rendered$src` records the R +# statement behind each generated line, so the message can quote that instead. +dsl_error_bullets <- function(message, rendered) { + core <- sub("^.*?error\\[DSL[0-9]+\\]:\\s*", "", message) + # Keep the first line only: the rest is the generated-source caret diagram. + core <- trimws(sub("\\n.*$", "", core)) + if (!nzchar(core)) { + core <- trimws(message) + } + + bullets <- c("x" = dsl_literal(core)) + + line <- NA_integer_ + found <- regmatches(message, regexpr("line ([0-9]+)", message)) + if (length(found) == 1L && nzchar(found)) { + line <- as.integer(sub("^line ", "", found)) + } + origin <- NA_character_ + if (!is.na(line) && line >= 1L && line <= length(rendered$src)) { + origin <- rendered$src[[line]] + } + + if (!is.na(origin) && nzchar(origin)) { + bullets <- c(bullets, "i" = dsl_literal(sprintf("From your model, %s", origin))) + } + c( + bullets, + "i" = "Check the block this line came from for an undeclared name, a duplicate assignment, or a value that is not a number." + ) +} + +# --------------------------------------------------------------------------- +# Render context +# --------------------------------------------------------------------------- +# Rendering carries a context so an expression can introduce helper derived +# values when the DSL cannot express an R construct in place. Helpers are +# emitted as ordinary derived assignments immediately before the statement that +# needs them, which keeps R's meaning without teaching the backend new syntax. + +dsl_ctx <- function(used = character(0)) { + ctx <- new.env(parent = emptyenv()) + ctx$used <- unique(c(used, dsl_language_words())) + ctx$index <- 0L + ctx$hoist <- list() + ctx +} + +# Words the DSL treats specially, so a generated helper can never shadow them. +dsl_language_words <- function() { + c( + "model", "kind", "ode", "analytical", "sde", "parameters", "constants", + "covariates", "states", "routes", "derive", "dynamics", "outputs", "init", + "drift", "diffusion", "particles", "if", "else", "for", "in", "let", + "true", "false" + ) +} + +# A helper name that cannot collide with anything the model already uses. +dsl_fresh_temp <- function(ctx) { + repeat { + ctx$index <- ctx$index + 1L + name <- sprintf("_pmt%d", ctx$index) + if (!name %in% ctx$used) { + ctx$used <- c(ctx$used, name) + return(name) + } + } +} + +# Emit `code` as a derived helper and return the helper's name. +dsl_hoist_code <- function(ctx, code) { + name <- dsl_fresh_temp(ctx) + ctx$hoist[[length(ctx$hoist) + 1L]] <- list(name = name, code = code) + name +} + +# Render `expr` into a helper derived value and return its name. +dsl_hoist <- function(ctx, expr, allow_if = TRUE) { + dsl_hoist_code(ctx, expr_to_dsl(expr, ctx, allow_if = allow_if)) +} + +# Consume the helpers accumulated while rendering one statement, labelled with +# the statement they belong to. +dsl_take_hoists <- function(ctx, src = NA_character_) { + hoists <- ctx$hoist + ctx$hoist <- list() + dsl_lines( + vapply(hoists, function(h) sprintf("%s = %s", h$name, h$code), character(1)), + rep(src, length(hoists)) + ) +} + +# --------------------------------------------------------------------------- +# Emitted lines and provenance +# --------------------------------------------------------------------------- +# Generated source is carried as a small value type that keeps every line and +# the R statement it came from together, so they cannot drift apart. Emitters +# concatenate parts with `dsl_c()` instead of `c()`; plain strings are accepted +# and count as having no origin (headers, separators, generated declarations). + +dsl_lines <- function(text = character(0), src = NA_character_) { + text <- as.character(text) + src <- if (length(src) == 1L) rep(as.character(src), length(text)) else as.character(src) + if (length(src) != length(text)) { + cli::cli_abort( + "internal error: {length(text)} generated line{?s} but {length(src)} source label{?s}", + call = NULL + ) + } + structure(list(text = text, src = src), class = "dsl_lines") +} + +dsl_text <- function(x) { + if (inherits(x, "dsl_lines")) x$text else as.character(x) +} + +dsl_src <- function(x) { + if (inherits(x, "dsl_lines")) x$src else rep(NA_character_, length(dsl_text(x))) +} + +dsl_c <- function(...) { + parts <- list(...) + parts <- parts[!vapply(parts, is.null, logical(1))] + if (length(parts) == 0L) { + return(dsl_lines()) + } + dsl_lines( + unlist(lapply(parts, dsl_text), use.names = FALSE), + unlist(lapply(parts, dsl_src), use.names = FALSE) + ) +} + +# Human-readable origin of a statement, used in diagnostics. `where` is the +# "(file:line)" suffix when R kept source references for the model. +dsl_statement_label <- function(block, expr, where = "") { + sprintf("%s%s: %s", block, where, dsl_expr_label(expr)) +} + +# "(file.R:12)" for a source reference, or "" when R did not keep one. A plain +# `Rscript` run has keep.source = FALSE, so models defined there have no line +# information; the diagnostic is still complete without it. +dsl_source_where <- function(ref, srcfile = NULL) { + if (is.null(ref)) { + return("") + } + file <- attr(ref, "srcfile") + if (is.null(file)) { + file <- srcfile + } + line <- suppressWarnings(as.integer(ref[[1L]])) + if (is.null(file) || length(line) != 1L || is.na(line)) { + return("") + } + name <- tryCatch(file$filename, error = function(e) NULL) + if (is.null(name) || length(name) != 1L || !nzchar(name) || identical(name, "")) { + return("") + } + sprintf(" (%s:%d)", basename(name), line) +} + +# Does this expression reference a dose input anywhere? +dsl_expr_has_route <- function(expr) { + length(eqn_route_inputs(expr, "bolus")) > 0L || + length(eqn_route_inputs(expr, "infusion")) > 0L +} + +# A dose input inside a conditional has no meaning: the route is inferred from +# the derivative RHS, and conditional dosing is not a supported event model. +# Report it against the user's expression instead of letting the helper hoist +# the input out of the derivative where it would silently stop being a dose. +dsl_check_conditional_routes <- function(expr) { + if (is.call(expr)) { + head <- dsl_call_head(expr) + if (!is.null(head) && head %in% c("if", "ifelse") && dsl_expr_has_route(expr)) { + cli::cli_abort( + c( + "x" = "A dose input cannot be inside a conditional: {.code {dsl_expr_label(expr)}}.", + "i" = "Put the conditional part in a derived value, then use one {.code bolus(input) * scale} term outside the conditional." + ), + call = NULL + ) + } + lapply(as.list(expr), dsl_check_conditional_routes) + } + invisible(NULL) +} + +# Fold n-ary rendered arguments into the binary form the DSL provides. +dsl_fold_binary <- function(fn, codes) { + acc <- codes[[1]] + for (i in seq_along(codes)[-1]) { + acc <- sprintf("%s(%s, %s)", fn, acc, codes[[i]]) + } + acc +} + +# R rounds halves to even; the DSL's `round` rounds halves away from zero. Emit +# R's rule explicitly rather than accept a different answer for exact halves. +dsl_round_helpers <- function(ctx, code, digits) { + scaled <- if (digits == 0L) code else sprintf("(%s * 10.0^%d.0)", code, digits) + value <- dsl_hoist_code(ctx, scaled) + base <- dsl_hoist_code(ctx, sprintf("floor(%s)", value)) + rest <- dsl_hoist_code(ctx, sprintf("%s - %s", value, base)) + # The conditional has to be a statement of its own: the DSL rejects a + # conditional wrapped in parentheses, so it cannot be scaled in place. + rounded <- dsl_hoist_code( + ctx, + sprintf( + "if (%s < 0.5) %s else if (%s > 0.5) %s + 1.0 else if (floor(%s / 2.0) * 2.0 == %s) %s else %s + 1.0", + rest, base, rest, base, base, base, base, base + ) + ) + if (digits == 0L) { + return(rounded) + } + dsl_hoist_code(ctx, sprintf("(%s / 10.0^%d.0)", rounded, digits)) +} + # Convert a pure R expression (no assignments or blocks) to a DSL expression. -expr_to_dsl <- function(expr, allow_if = TRUE) { +expr_to_dsl <- function(expr, ctx = NULL, allow_if = TRUE) { # Numeric literals: emit integers with a trailing `.0` so the DSL treats them # as floating point, matching the Rust transpiler's behaviour. + if (is.numeric(expr) && length(expr) == 1 && is.na(expr)) { + dsl_unsupported( + expr, "Model expressions must be finite numbers, but this one is missing:", + "Remove the missing value or replace it with a number." + ) + } + if (is.numeric(expr) && length(expr) == 1 && !is.finite(expr)) { + dsl_unsupported(expr, "Model expressions must be finite numbers, but this one is not:", "Use a large finite value such as 1e30 instead of Inf or NaN.") + } if (is.numeric(expr) && length(expr) == 1) { val <- expr if (is.finite(val) && val == floor(val) && abs(val) <= .Machine$integer.max) { @@ -57,14 +361,31 @@ expr_to_dsl <- function(expr, allow_if = TRUE) { } return(as.character(val)) } + # R's logical constants map to the DSL's boolean literals, which are valid in + # conditions. They are not numbers, so the DSL rejects them as values. + if (is.logical(expr) && length(expr) == 1) { + if (is.na(expr)) { + dsl_unsupported( + expr, "Missing values cannot be used in a model:", + "Remove the missing value or replace it with a number." + ) + } + return(if (expr) "true" else "false") + } if (is.symbol(expr)) { return(tolower(as.character(expr))) } if (!is.call(expr)) { - cli::cli_abort("Unable to convert model expression to DSL: unsupported token.") + dsl_unsupported(expr, "This token cannot be used in a model expression:") } - op <- as.character(expr[[1]]) + op <- dsl_call_head(expr) + if (is.null(op)) { + dsl_unsupported( + expr, "This call cannot be translated:", + "Use a plain function such as exp(x), and write arithmetic with operators." + ) + } # Indexing: x[i] -> x{i}. Only literal, positive integer indices are allowed. if (op == "[") { @@ -77,46 +398,147 @@ expr_to_dsl <- function(expr, allow_if = TRUE) { )) } idx <- as.integer(idx_raw) - if (var %in% c("b", "bolus", "rateiv", "r")) { + # Indexing anything else used to be rewritten into a different identifier + # (`p1[1]` became `p11`), which then failed as an unknown name. + if (!var %in% c("x", "b", "bolus", "rateiv", "r")) { + cli::cli_abort( + c( + "x" = "{.code {var}[{idx}]} cannot be used as a value here.", + "i" = "Inside an expression, only {.code x[i]} (a state) and {.code b[i]}/{.code bolus[i]}, {.code r[i]}/{.code rateiv[i]} (dose inputs) can be indexed.", + "i" = "Assign to {.code lag[i]}, {.code fa[i]}, {.code dx[i]} or {.code y[i]} with a statement of its own." + ), + call = NULL + ) + } + if (var %in% c("b", "bolus")) { cli::cli_abort(c( - "x" = "Bolus/infusion inputs may only be used as standalone additive terms in derivative equations.", - "i" = "Write, for example, {.code dx[1] = -ke * x[1] + rateiv[1]}, not inside a product." + "x" = "Bolus inputs may only be standalone additive terms or one exact product with a scale.", + "i" = "Use {.code B[1]}, {.code B[1] * scale}, or {.code scale * B[1]}." )) } - return(sprintf("%s%d", var, idx)) - } - - # Conditionals map to `if (cond) then else else`. The pharmsol DSL authoring - # surface only accepts a conditional as an entire equation right-hand side or - # as the `else` branch of another conditional (an `else if` chain). Any other - # position (nested in an operator/function, or in the `then` branch) is - # rejected by the DSL parser, so we catch it here to give an R-level error - # that points at the model instead of a cryptic parse error on generated code. - if (op == "if") { - if (!allow_if) { + if (var %in% c("rateiv", "r")) { cli::cli_abort(c( - "x" = "A conditional {.code if (...) ... else ...} can only be a whole equation right-hand side.", - "i" = "It cannot be nested inside another expression (e.g. {.code 2 * if (...)}) or in the {.code then} branch.", - "i" = "Assign it to a secondary variable first, e.g. {.code tmp = if (cond) a else b}, then use {.code tmp}." + "x" = "Infusion inputs may only be standalone additive terms or one exact product with a scale.", + "i" = "Use {.code R[1]}, {.code R[1] * scale}, or {.code scale * R[1]}." )) } + return(sprintf("%s%d", var, idx)) + } + + # R conditionals are expressions, but the DSL can only place one as a whole + # assignment right-hand side. Anywhere else the value is hoisted into a + # derived helper, so nested conditionals work exactly as they do in R. + if (op %in% c("if", "ifelse")) { args <- as.list(expr[-1]) if (length(args) != 3) { - cli::cli_abort(c( - "x" = "Conditional expressions in the DSL must include an `else` branch.", - "i" = "Write {.code if (cond) a else b}." - )) + dsl_unsupported( + expr, "Conditional expressions must include an else branch:", + "Write if (cond) a else b." + ) } - # Only the `else` branch may itself be a conditional (right-associative - # `else if` chains); the condition and `then` branch may not. - cond <- expr_to_dsl(args[[1]], allow_if = FALSE) - then_code <- expr_to_dsl(args[[2]], allow_if = FALSE) - else_code <- expr_to_dsl(args[[3]], allow_if = TRUE) + dsl_check_conditional_routes(expr) + if (!allow_if) { + if (is.null(ctx)) { + dsl_unsupported( + expr, "A conditional cannot be nested inside another expression:", + "Assign it to a secondary variable first, for example tmp = if (cond) a else b." + ) + } + return(dsl_hoist(ctx, expr, allow_if = TRUE)) + } + # Only the else branch may itself be a conditional (right-associative + # else-if chains); nested conditionals in the condition or the then branch + # are hoisted into helpers. + cond <- expr_to_dsl(args[[1]], ctx, allow_if = FALSE) + then_code <- expr_to_dsl(args[[2]], ctx, allow_if = FALSE) + else_code <- expr_to_dsl(args[[3]], ctx, allow_if = TRUE) return(sprintf("if (%s) %s else %s", cond, then_code, else_code)) } args <- as.list(expr[-1]) - a <- lapply(args, function(x) expr_to_dsl(x, allow_if = FALSE)) + a <- lapply(args, function(x) expr_to_dsl(x, ctx, allow_if = FALSE)) + + # Legacy spellings carry a migration message, because the replacement changes + # what the call means and not only its name. + legacy <- dsl_legacy_functions() + if (op %in% names(legacy)) { + cli::cli_abort( + c( + "x" = "Unsupported function {.val {op}} in {.code {dsl_expr_label(expr)}}.", + "i" = legacy[[op]], + "i" = "The runtime computes w = alpha * u * v, so check that the value you used to pass as w followed the same convention." + ), + call = NULL + ) + } + + # Runtime-only effect functions. Arity is checked here so a wrong call is + # reported against the R expression rather than generated DSL text. + effects <- dsl_effect_functions() + if (op %in% names(effects)) { + expected <- unname(effects[[op]]) + if (length(a) != expected) { + cli::cli_abort( + c( + "x" = "{.fn {op}} takes {expected} arguments, but {.code {dsl_expr_label(expr)}} has {length(a)}.", + "i" = "Write {.code {dsl_effect_signature()[[op]]}}.", + "i" = "u and v are normalised exposures, alpha is the interaction coefficient, and h1/h2 are the Hill exponents." + ), + call = NULL + ) + } + return(sprintf("%s(%s)", op, paste(unlist(a), collapse = ", "))) + } + + # pmin/pmax are elementwise in R; for the scalar expressions a model uses, + # that is min/max. The DSL's min/max are binary while R's are variadic, and + # trunc/sign are not DSL intrinsics at all. + if (op %in% c("pmin", "pmax")) { + if (length(a) < 2) { + dsl_unsupported(expr, stringr::str_glue("{op} needs at least two arguments:"), stringr::str_glue("Write {op}(a, b).")) + } + return(dsl_fold_binary(if (op == "pmin") "min" else "max", a)) + } + if (op %in% c("max", "min") && length(a) != 2) { + if (length(a) == 1) { + return(a[[1]]) + } + if (length(a) == 0) { + dsl_unsupported(expr, stringr::str_glue("{op} needs at least one argument:")) + } + return(dsl_fold_binary(op, a)) + } + if (op == "trunc") { + if (length(a) != 1) dsl_unsupported(expr, "trunc takes one argument:") + value <- dsl_hoist_code(ctx, a[[1]]) + return(dsl_hoist_code(ctx, sprintf("if (%s < 0.0) ceil(%s) else floor(%s)", value, value, value))) + } + if (op == "sign") { + if (length(a) != 1) dsl_unsupported(expr, "sign takes one argument:") + value <- dsl_hoist_code(ctx, a[[1]]) + # Multiplying by zero keeps NaN propagating, as sign() does in R. + return(dsl_hoist_code(ctx, sprintf( + "if (%s < 0.0) -1.0 else if (%s > 0.0) 1.0 else (%s * 0.0)", value, value, value + ))) + } + if (op == "round") { + if (length(a) == 1) { + return(dsl_round_helpers(ctx, a[[1]], 0L)) + } + digits <- if (length(args) >= 2) dsl_literal_number(args[[2]]) else NULL + if (is.null(digits) || digits != floor(digits) || abs(digits) > 15) { + dsl_unsupported( + expr, "round needs a whole number of digits between -15 and 15:", + "Write the number of digits literally, for example round(x, 2)." + ) + } + if (length(a) != 2) dsl_unsupported(expr, "round takes a value and an optional number of digits:") + return(dsl_round_helpers(ctx, a[[1]], as.integer(digits))) + } + if (op == "log" && length(a) == 2) { + # R's log(x, base); the DSL has no two-argument logarithm. + return(sprintf("(ln(%s) / ln(%s))", a[[1]], a[[2]])) + } fmap <- dsl_function_map() @@ -140,13 +562,16 @@ expr_to_dsl <- function(expr, allow_if = TRUE) { "!" = sprintf("!(%s)", a[[1]]), { # Function call: look up in the DSL intrinsic map. - if (!is.null(fmap[[op]])) { + if (op %in% names(fmap)) { sprintf("%s(%s)", fmap[[op]], paste(unlist(a), collapse = ", ")) } else { - cli::cli_abort(c( - "x" = "Unsupported function {.val {op}} for the DSL backend.", - "i" = "Supported functions are: {.val {sort(unique(unname(dsl_function_map())))}}." - )) + dsl_unsupported( + expr, stringr::str_glue("Unsupported function '{op}':"), + stringr::str_glue( + "Supported functions are {paste(sort(unique(unname(dsl_function_map()))), collapse = ', ')}, " + , "plus the runtime effect functions {paste(names(dsl_effect_functions()), collapse = ', ')}." + ) + ) } } ) @@ -156,11 +581,22 @@ expr_to_dsl <- function(expr, allow_if = TRUE) { # Return the top-level statements of a model block function body. dsl_body_stmts <- function(fun) { b <- body(fun) - if (is.call(b) && as.character(b[[1]]) == "{") { - as.list(b[-1]) + refs <- attr(b, "srcref") + srcfile <- attr(b, "srcfile") + if (identical(dsl_call_head(b), "{")) { + stmts <- as.list(b[-1]) + # The source references of a braced body include the brace itself. + refs <- if (length(refs) >= length(stmts) + 1L) refs[seq_along(stmts) + 1L] else NULL } else { - list(b) + stmts <- list(b) + refs <- if (length(refs) >= 1L) refs[1] else NULL } + lapply(seq_along(stmts), function(i) { + list( + expr = stmts[[i]], + where = if (is.null(refs)) "" else dsl_source_where(refs[[i]], srcfile) + ) + }) } # Is `expr` an assignment (`<-` or `=`)? @@ -189,7 +625,7 @@ dsl_flatten_add <- function(expr, sign = 1) { # If `expr` is a route reference (`b[j]`, `bolus[j]`, `rateiv[j]`, `r[j]`), # return `list(kind = "bolus"|"infusion", input = j)`, else NULL. dsl_route_of <- function(expr) { - if (is.call(expr) && as.character(expr[[1]]) == "[") { + if (identical(dsl_call_head(expr), "[")) { v <- tolower(as.character(expr[[2]])) idx <- expr[[3]] if (v %in% c("b", "bolus", "rateiv", "r") && is.numeric(idx) && length(idx) == 1) { @@ -200,6 +636,36 @@ dsl_route_of <- function(expr) { NULL } +# Recognize one exact linear route product, with the route on either side. +dsl_scaled_route_of <- function(expr) { + if (!is.call(expr) || !identical(expr[[1]], as.name("*")) || length(expr) != 3) { + return(NULL) + } + + left <- dsl_route_of(expr[[2]]) + right <- dsl_route_of(expr[[3]]) + if (is.null(left) == is.null(right)) { + return(NULL) + } + + route <- if (!is.null(left)) left else right + route$scale <- if (!is.null(left)) expr[[3]] else expr[[2]] + if (route$kind == "bolus") { + uses_state_or_input <- function(value) { + if (!is.call(value)) return(FALSE) + if (identical(value[[1]], as.name("["))) { + name <- tolower(as.character(value[[2]])) + if (name %in% c("x", "dx", "b", "bolus", "r", "rateiv")) return(TRUE) + } + any(vapply(as.list(value)[-1], uses_state_or_input, logical(1))) + } + if (uses_state_or_input(route$scale)) { + cli::cli_abort("A bolus scale cannot depend directly on state or dose inputs.") + } + } + route +} + # Return the sorted, unique data input indices referenced by bolus (`b[]` / # `bolus[]`) or infusion (`rateiv[]` / `r[]`) terms in a model equation function. # Used to validate that the model represents every dose input present in the data. @@ -222,20 +688,44 @@ eqn_route_inputs <- function(fun, kind = c("bolus", "infusion")) { sort(unique(found)) } +# Every identifier a block assigns or reads as a plain variable, so generated +# helper names can never collide with the model's own names. +dsl_assigned_names <- function(fun) { + if (is.null(fun) || !is.function(fun)) { + return(character(0)) + } + found <- character(0) + walk <- function(expr) { + if (is.call(expr)) { + if (dsl_is_assign(expr)) { + lhs <- expr[[2]] + if (is.symbol(lhs)) { + found <<- c(found, tolower(as.character(lhs))) + } else if (is.call(lhs) && identical(lhs[[1]], as.name("["))) { + found <<- c(found, tolower(as.character(lhs[[2]]))) + } + } + lapply(as.list(expr), walk) + } + invisible(NULL) + } + walk(body(fun)) + unique(found) +} + # Build a DSL expression string from a list of signed terms. -dsl_join_terms <- function(terms) { +dsl_join_terms <- function(terms, ctx) { if (length(terms) == 0) { return("0.0") } - # A conditional is only valid as a whole equation right-hand side, so it may - # appear only as a lone, positively-signed term (emitted bare). When summed - # with other terms (or negated) it must raise the R-level guard error instead - # of emitting DSL the backend will reject. + # A conditional may be emitted in place only when it is the whole right-hand + # side. Anywhere else it is hoisted into a derived helper, so summing or + # negating a conditional is fine. allow_if <- length(terms) == 1L && terms[[1]]$sign > 0 pieces <- character(0) for (i in seq_along(terms)) { t <- terms[[i]] - es <- expr_to_dsl(t$expr, allow_if = allow_if) + es <- if (!is.null(t$dsl)) t$dsl else expr_to_dsl(t$expr, ctx, allow_if = allow_if) if (i == 1) { pieces <- if (t$sign < 0) sprintf("-(%s)", es) else es } else { @@ -245,14 +735,18 @@ dsl_join_terms <- function(terms) { pieces } -# Convert a derivative RHS into (routes, stripped DSL expression) for the given +# Convert a derivative RHS into (routes, explicit-input DSL expression) for the given # destination compartment index `comp`. -dsl_extract_routes <- function(rhs, comp) { +dsl_extract_routes <- function(rhs, comp, ctx, fa_values = NULL) { + dsl_check_conditional_routes(rhs) terms <- dsl_flatten_add(rhs) routes <- list() kept <- list() for (t in terms) { route <- dsl_route_of(t$expr) + if (is.null(route)) { + route <- dsl_scaled_route_of(t$expr) + } if (!is.null(route)) { if (t$sign < 0) { cli::cli_abort(c( @@ -262,42 +756,69 @@ dsl_extract_routes <- function(rhs, comp) { } route$comp <- comp routes[[length(routes) + 1]] <- route + input <- sprintf("%s(%s)", route$kind, pm_input_label(route$input)) + # The fa block scales the bolus dose. Applying it while the term is + # rendered keeps generated text produced in a single pass; patching it + # into the finished text afterwards depended on the term appearing once. + if (route$kind == "bolus" && !is.null(fa_values)) { + factor <- fa_values[[as.character(route$input)]] + if (!is.null(factor)) { + input <- sprintf("%s * (%s)", input, factor) + } + } + if (!is.null(route$scale)) { + input <- sprintf("%s * (%s)", input, expr_to_dsl(route$scale, ctx, allow_if = FALSE)) + } + kept[[length(kept) + 1]] <- list(sign = t$sign, dsl = input) } else { kept[[length(kept) + 1]] <- t } } - list(routes = routes, expr = dsl_join_terms(kept)) + list(routes = routes, expr = dsl_join_terms(kept, ctx)) } # Emit the ODE equation block: returns routes, derived assignments, and dx lines. -dsl_eqn_block <- function(fun) { +dsl_eqn_block <- function(fun, ctx, fa_values = NULL) { exprs <- dsl_body_stmts(fun) routes <- list() - derived <- character(0) - dx_lines <- character(0) - - for (e in exprs) { + derived <- dsl_lines() + dx_lines <- dsl_lines() + + for (stmt in exprs) { + e <- stmt$expr + where <- stmt$where + if (identical(dsl_call_head(e), "if")) { + derived <- dsl_c(derived, dsl_if_statement(e, ctx, "equation block")) + next + } if (!dsl_is_assign(e)) { - cli::cli_abort(c( - "x" = "Only assignments are supported in the equation block for the DSL backend.", - "i" = "Control-flow constructs such as {.code for} loops are not supported." - )) + dsl_unsupported( + e, "Only assignments are supported here:", + "Loops are not part of the model language; write the assignment for each case instead." + ) } + label <- dsl_statement_label("equation block", e, where) lhs <- e[[2]] rhs <- e[[3]] - if (is.call(lhs) && as.character(lhs[[1]]) == "[") { + if (identical(dsl_call_head(lhs), "[")) { tgt <- tolower(as.character(lhs[[2]])) idx <- as.integer(lhs[[3]]) if (tgt != "dx") { - cli::cli_abort("Unexpected indexed assignment to {.code {tgt}[{idx}]} in equation block.") + dsl_unsupported(e, sprintf("Unexpected assignment to %s[%s] in the equation block:", tgt, idx)) } - res <- dsl_extract_routes(rhs, idx) + res <- dsl_extract_routes(rhs, idx, ctx, fa_values) routes <- c(routes, res$routes) - dx_lines <- c(dx_lines, sprintf("dx(x%d) = %s", idx, res$expr)) + derived <- dsl_c(derived, dsl_take_hoists(ctx, label)) + dx_lines <- dsl_c(dx_lines, dsl_lines(sprintf("dx(x%d) = %s", idx, res$expr), label)) } else { # Scalar (secondary/derived) assignment. name <- tolower(as.character(lhs)) - derived <- c(derived, sprintf("%s = %s", name, expr_to_dsl(rhs))) + code <- expr_to_dsl(rhs, ctx, allow_if = TRUE) + derived <- dsl_c( + derived, + dsl_take_hoists(ctx, label), + dsl_lines(sprintf("%s = %s", name, code), label) + ) } } @@ -305,49 +826,143 @@ dsl_eqn_block <- function(fun) { } # Emit the output block: returns derived assignments and out() lines. -dsl_out_block <- function(fun) { +dsl_out_block <- function(fun, ctx) { exprs <- dsl_body_stmts(fun) - derived <- character(0) - out_lines <- character(0) - - for (e in exprs) { + derived <- dsl_lines() + out_lines <- dsl_lines() + + for (stmt in exprs) { + e <- stmt$expr + where <- stmt$where + if (identical(dsl_call_head(e), "if")) { + derived <- dsl_c(derived, dsl_if_statement(e, ctx, "output block")) + next + } if (!dsl_is_assign(e)) { - cli::cli_abort("Only assignments are supported in the output block for the DSL backend.") + dsl_unsupported(e, "Only assignments are supported in the output block:") } + label <- dsl_statement_label("output block", e, where) lhs <- e[[2]] rhs <- e[[3]] - if (is.call(lhs) && as.character(lhs[[1]]) == "[") { + if (identical(dsl_call_head(lhs), "[")) { tgt <- tolower(as.character(lhs[[2]])) idx <- as.integer(lhs[[3]]) if (tgt != "y") { - cli::cli_abort("Unexpected indexed assignment to {.code {tgt}[{idx}]} in output block.") + dsl_unsupported(e, sprintf("Unexpected assignment to %s[%s] in the output block:", tgt, idx)) } # Pmetrics data uses numeric OUTEQ identifiers, so output `y[i]` maps to # the canonical pharmsol output label for `i`. - out_lines <- c(out_lines, sprintf("out(%s) = %s", pm_output_label(idx), expr_to_dsl(rhs))) + code <- expr_to_dsl(rhs, ctx, allow_if = TRUE) + derived <- dsl_c(derived, dsl_take_hoists(ctx, label)) + out_lines <- dsl_c( + out_lines, + dsl_lines(sprintf("out(%s) = %s", pm_output_label(idx), code), label) + ) } else { name <- tolower(as.character(lhs)) - derived <- c(derived, sprintf("%s = %s", name, expr_to_dsl(rhs))) + code <- expr_to_dsl(rhs, ctx, allow_if = TRUE) + derived <- dsl_c( + derived, + dsl_take_hoists(ctx, label), + dsl_lines(sprintf("%s = %s", name, code), label) + ) } } list(derived = derived, out = out_lines) } +# Translate a braced R `if` statement. +# +# R users write a conditional block to choose between alternative values, and +# the DSL accepts a plain-assignment `if` statement, so the block is passed +# through. A branch that assigns something other than a derived scalar (a +# derivative, an output, an initial condition) cannot live inside an `if` +# statement, and the message says what to write instead. +dsl_if_statement <- function(expr, ctx, block_label) { + args <- as.list(expr[-1]) + if (length(args) < 2L || length(args) > 3L) { + dsl_unsupported( + expr, "A conditional needs a condition and a body:", + "Write if (cond) { ... } else { ... }." + ) + } + label <- dsl_statement_label(block_label, expr) + condition <- expr_to_dsl(args[[1]], ctx, allow_if = FALSE) + + branch_lines <- function(branch) { + stmts <- if (identical(dsl_call_head(branch), "{")) as.list(branch[-1]) else list(branch) + out <- dsl_lines() + for (s in stmts) { + if (!dsl_is_assign(s)) { + dsl_unsupported( + s, sprintf("Only assignments are supported inside a conditional in the %s:", block_label), + "Write the assignment for each case instead of using a loop." + ) + } + lhs <- s[[2]] + if (!is.symbol(lhs)) { + dsl_unsupported( + s, sprintf("A conditional in the %s here can only choose between scalar values:", block_label), + "Make the equation itself conditional, for example dx[1] <- if (cond) a else b." + ) + } + stmt_label <- dsl_statement_label(block_label, s) + code <- expr_to_dsl(s[[3]], ctx, allow_if = TRUE) + out <- dsl_c( + out, + dsl_take_hoists(ctx, stmt_label), + dsl_lines(sprintf("%s = %s", tolower(as.character(lhs)), code), stmt_label) + ) + } + out + } + + then_lines <- branch_lines(args[[2]]) + else_lines <- if (length(args) == 3L) branch_lines(args[[3]]) else dsl_lines() + + indent <- function(x) ifelse(nzchar(x), paste0(" ", x), x) + body <- dsl_c(dsl_lines(sprintf("if (%s) {", condition), label)) + body <- dsl_c(body, dsl_lines(indent(dsl_text(then_lines)), dsl_src(then_lines))) + if (length(dsl_text(else_lines)) > 0L) { + body <- dsl_c( + body, + dsl_lines("} else {", label), + dsl_lines(indent(dsl_text(else_lines)), dsl_src(else_lines)) + ) + } + dsl_c(body, dsl_lines("}", label)) +} + # Emit `derive` assignments from a secondary-equation block. -dsl_sec_block <- function(fun) { +dsl_sec_block <- function(fun, ctx) { exprs <- dsl_body_stmts(fun) - derived <- character(0) - for (e in exprs) { + derived <- dsl_lines() + for (stmt in exprs) { + e <- stmt$expr + where <- stmt$where + if (identical(dsl_call_head(e), "if")) { + derived <- dsl_c(derived, dsl_if_statement(e, ctx, "secondary block")) + next + } if (!dsl_is_assign(e)) { - cli::cli_abort("Only assignments are supported in the secondary block for the DSL backend.") + dsl_unsupported( + e, "Only assignments are supported in the secondary block:", + "Loops are not part of the model language; write one assignment per value." + ) } lhs <- e[[2]] rhs <- e[[3]] if (!is.symbol(lhs)) { - cli::cli_abort("Secondary equations must assign to scalar variables.") + dsl_unsupported(e, "Secondary equations must assign to a scalar variable:") } - derived <- c(derived, sprintf("%s = %s", tolower(as.character(lhs)), expr_to_dsl(rhs))) + label <- dsl_statement_label("secondary block", e, where) + code <- expr_to_dsl(rhs, ctx, allow_if = TRUE) + derived <- dsl_c( + derived, + dsl_take_hoists(ctx, label), + dsl_lines(sprintf("%s = %s", tolower(as.character(lhs)), code), label) + ) } derived } @@ -355,50 +970,89 @@ dsl_sec_block <- function(fun) { # Emit route-property modifiers (`lag(...)` / `fa(...)`) from a lag/fa block. # `target` is the DSL property name ("lag" or "fa"); the R block assigns to # `lag[j]` / `fa[j]` where `j` is the 1-based input index. -dsl_route_property_block <- function(fun, target) { - exprs <- dsl_body_stmts(fun) - derived <- character(0) - lines <- character(0) - for (e in exprs) { +dsl_route_property_block <- function(fun, target, ctx) { + exprs <- if (is.null(fun)) list() else dsl_body_stmts(fun) + derived <- dsl_lines() + lines <- dsl_lines() + values <- list() + seen_inputs <- character(0) + for (stmt in exprs) { + e <- stmt$expr + where <- stmt$where if (!dsl_is_assign(e)) { cli::cli_abort("Only assignments are supported in the {target} block for the DSL backend.") } + label <- dsl_statement_label(paste(target, "block"), e, where) lhs <- e[[2]] rhs <- e[[3]] - if (is.call(lhs) && as.character(lhs[[1]]) == "[") { + if (identical(dsl_call_head(lhs), "[")) { tgt <- tolower(as.character(lhs[[2]])) idx <- as.integer(lhs[[3]]) if (tgt != target) { cli::cli_abort("Unexpected indexed assignment to {.code {tgt}[{idx}]} in {target} block.") } - lines <- c(lines, sprintf("%s(%s) = %s", target, pm_input_label(idx), expr_to_dsl(rhs))) + key <- as.character(idx) + if (key %in% seen_inputs) { + cli::cli_abort("The {target} block assigns input {idx} more than once.") + } + seen_inputs <- c(seen_inputs, key) + # A route property cannot take a conditional directly, so any conditional + # is hoisted into an event-safe derived value first. + values[[key]] <- expr_to_dsl(rhs, ctx, allow_if = FALSE) + derived <- dsl_c(derived, dsl_take_hoists(ctx, label)) + lines <- dsl_c( + lines, + dsl_lines(sprintf("%s(%s) = %s", target, pm_input_label(idx), values[[key]]), label) + ) } else { - derived <- c(derived, sprintf("%s = %s", tolower(as.character(lhs)), expr_to_dsl(rhs))) + name <- tolower(as.character(lhs)) + code <- expr_to_dsl(rhs, ctx, allow_if = TRUE) + derived <- dsl_c( + derived, + dsl_take_hoists(ctx, label), + dsl_lines(sprintf("%s = %s", name, code), label) + ) } } - list(derived = derived, lines = lines) + + list(derived = derived, lines = lines, values = values) } # Emit `init(...)` statements from an initial-conditions block. -dsl_ini_block <- function(fun) { +dsl_ini_block <- function(fun, ctx) { exprs <- dsl_body_stmts(fun) - derived <- character(0) - lines <- character(0) - for (e in exprs) { + derived <- dsl_lines() + lines <- dsl_lines() + for (stmt in exprs) { + e <- stmt$expr + where <- stmt$where + if (identical(dsl_call_head(e), "if")) { + derived <- dsl_c(derived, dsl_if_statement(e, ctx, "initial-conditions block")) + next + } if (!dsl_is_assign(e)) { cli::cli_abort("Only assignments are supported in the initial-conditions block for the DSL backend.") } + label <- dsl_statement_label("initial-conditions block", e, where) lhs <- e[[2]] rhs <- e[[3]] - if (is.call(lhs) && as.character(lhs[[1]]) == "[") { + if (identical(dsl_call_head(lhs), "[")) { tgt <- tolower(as.character(lhs[[2]])) idx <- as.integer(lhs[[3]]) if (tgt != "x") { cli::cli_abort("Unexpected indexed assignment to {.code {tgt}[{idx}]} in initial-conditions block.") } - lines <- c(lines, sprintf("init(x%d) = %s", idx, expr_to_dsl(rhs))) + code <- expr_to_dsl(rhs, ctx, allow_if = TRUE) + derived <- dsl_c(derived, dsl_take_hoists(ctx, label)) + lines <- dsl_c(lines, dsl_lines(sprintf("init(x%d) = %s", idx, code), label)) } else { - derived <- c(derived, sprintf("%s = %s", tolower(as.character(lhs)), expr_to_dsl(rhs))) + name <- tolower(as.character(lhs)) + code <- expr_to_dsl(rhs, ctx, allow_if = TRUE) + derived <- dsl_c( + derived, + dsl_take_hoists(ctx, label), + dsl_lines(sprintf("%s = %s", name, code), label) + ) } } list(derived = derived, lines = lines) @@ -406,12 +1060,12 @@ dsl_ini_block <- function(fun) { # Finalize route usages into concrete DSL routes. # -# Each usage is `list(kind, input, comp)` where `input` is the identifier used in -# the model equations (`b[j]` / `rateiv[j]`). Route labels are unique *per kind* +# Each usage is `list(kind, input, comp, scale)` where `scale` is present for an +# exact route-input product. Route labels are unique *per kind* # in the pharmsol DSL, so the same label may declare both a bolus and an # infusion; no relabelling of the data is ever required. # -# Returns a list of `{kind, label, comp}` route declarations. +# Returns a list of `{kind, input, label, comp, scale}` route declarations. dsl_finalize_routes <- function(routes) { by_input <- list() seen_inputs <- integer(0) @@ -430,10 +1084,8 @@ dsl_finalize_routes <- function(routes) { grp <- by_input[[as.character(inp)]] for (kd in c("bolus", "infusion")) { - comps <- unique(vapply( - Filter(function(r) identical(r$kind, kd), grp), - function(r) as.integer(r$comp), integer(1) - )) + kind_routes <- Filter(function(r) identical(r$kind, kd), grp) + comps <- unique(vapply(kind_routes, function(r) as.integer(r$comp), integer(1))) if (length(comps) == 0) next if (length(comps) > 1) { cli::cli_abort(c( @@ -442,8 +1094,17 @@ dsl_finalize_routes <- function(routes) { )) } + scaled <- Filter(function(r) !is.null(r$scale), kind_routes) + if ((kd == "bolus" || length(scaled) > 0) && length(kind_routes) > 1) { + cli::cli_abort(c( + "x" = "The {kd} input {inp} is used more than once.", + "i" = "Use one input term with a combined scale." + )) + } + scale <- if (length(scaled) == 1) scaled[[1]]$scale else NULL + final_routes[[length(final_routes) + 1L]] <- list( - kind = kd, label = pm_input_label(inp), comp = comps[[1]] + kind = kd, input = inp, label = pm_input_label(inp), comp = comps[[1]], scale = scale ) } } @@ -544,6 +1205,17 @@ dsl_analytical_param_map <- function(structure) { # Assemble the full pharmsol DSL text for a PM_model object. model_to_dsl <- function(model) { + dsl_render(model)$text +} + +# The same rendering, plus the R statement behind each emitted line. Used to +# report backend diagnostics against the user's model instead of generated text. +model_to_dsl_traced <- function(model) { + dsl_render(model) +} + +# Render a PM_model to DSL text, keeping per-line provenance. +dsl_render <- function(model) { arg_list <- model$arg_list model_list <- model$model_list if (is.null(arg_list) || is.null(model_list)) { @@ -554,6 +1226,16 @@ model_to_dsl <- function(model) { parameters <- tolower(names(arg_list$pri)) covariate_names <- tolower(names(arg_list$cov)) + # Generated helper names must never collide with a name the model uses. + ctx <- dsl_ctx(used = c( + parameters, + covariate_names, + unlist(lapply( + list(arg_list$sec, arg_list$eqn, arg_list$out, arg_list$ini, arg_list$lag, arg_list$fa), + dsl_assigned_names + )) + )) + header <- character(0) header <- c(header, sprintf("name = %s", if (is.null(model_list$name)) "user" else model_list$name)) header <- c(header, sprintf("kind = %s", if (type == "Analytical") "analytical" else "ode")) @@ -572,41 +1254,56 @@ model_to_dsl <- function(model) { } # Derived (secondary) equations shared across blocks. - derived <- character(0) + derived <- dsl_lines() if (!is.null(arg_list$sec)) { - derived <- c(derived, dsl_sec_block(arg_list$sec)) + derived <- dsl_c(derived, dsl_sec_block(arg_list$sec, ctx)) } if (type == "Analytical") { - return(dsl_analytical(model, header, derived, parameters)) + return(dsl_analytical(model, header, derived, parameters, ctx)) } # ---- ODE model ---- - eqn <- dsl_eqn_block(arg_list$eqn) - derived <- c(derived, eqn$derived) + # The fa block is read first because its values multiply the bolus terms + # while the derivatives are rendered, rather than being patched into the + # generated text afterwards. + fa_values <- NULL + if (!is.null(arg_list$fa)) { + fa <- dsl_route_property_block(arg_list$fa, "fa", ctx) + derived <- dsl_c(derived, fa$derived) + fa_values <- fa$values + } + + eqn <- dsl_eqn_block(arg_list$eqn, ctx, fa_values) + derived <- dsl_c(derived, eqn$derived) + routes <- dsl_finalize_routes(eqn$routes) - out <- dsl_out_block(arg_list$out) - derived <- c(derived, out$derived) + out <- dsl_out_block(arg_list$out, ctx) + derived <- dsl_c(derived, out$derived) - init_lines <- character(0) + init_lines <- dsl_lines() if (!is.null(arg_list$ini)) { - ini <- dsl_ini_block(arg_list$ini) - derived <- c(derived, ini$derived) + ini <- dsl_ini_block(arg_list$ini, ctx) + derived <- dsl_c(derived, ini$derived) init_lines <- ini$lines } - lag_lines <- character(0) + lag_lines <- dsl_lines() if (!is.null(arg_list$lag)) { - lag <- dsl_route_property_block(arg_list$lag, "lag") - derived <- c(derived, lag$derived) + lag <- dsl_route_property_block(arg_list$lag, "lag", ctx) + derived <- dsl_c(derived, lag$derived) lag_lines <- lag$lines } - fa_lines <- character(0) - if (!is.null(arg_list$fa)) { - fa <- dsl_route_property_block(arg_list$fa, "fa") - derived <- c(derived, fa$derived) - fa_lines <- fa$lines + # Every fa value must multiply a declared bolus term. + if (!is.null(fa_values)) { + bolus_inputs <- vapply(Filter(function(r) r$kind == "bolus", routes), + function(r) as.character(r$input), character(1)) + for (key in names(fa_values)) { + if (!key %in% bolus_inputs) { + cli::cli_abort("The fa block references input {key}, which has no bolus term.") + } + } } # Number of states and outputs. @@ -621,38 +1318,35 @@ model_to_dsl <- function(model) { states <- paste0("x", seq_len(n_states)) outputs <- pm_output_label(seq_len(n_out)) - routes <- dsl_finalize_routes(eqn$routes) - route_lines <- vapply(routes, function(r) { - sprintf("%s(%s) -> x%d", r$kind, r$label, r$comp) - }, character(1)) - - # Assemble the DSL text in an order that respects definite assignment: - # declarations, routes, derived values, route properties, initial conditions, - # derivatives, and finally outputs. - lines <- c( - header, - sprintf("states = %s", paste(states, collapse = ", ")), - sprintf("outputs = %s", paste(outputs, collapse = ", ")), - "", - route_lines, - if (length(route_lines) > 0) "" else NULL, + # Inputs remain explicit in the derivatives; the DSL infers their routes. + # Derived values precede lag, initial conditions, and derivatives. + body <- dsl_c( + dsl_lines(header, "header"), + dsl_lines(sprintf("states = %s", paste(states, collapse = ", "))), + dsl_lines(sprintf("outputs = %s", paste(outputs, collapse = ", "))), + dsl_lines(""), derived, - if (length(derived) > 0) "" else NULL, + if (length(dsl_text(derived)) > 0) dsl_lines("") else NULL, lag_lines, - fa_lines, - if (length(lag_lines) > 0 || length(fa_lines) > 0) "" else NULL, + if (length(dsl_text(lag_lines)) > 0) dsl_lines("") else NULL, init_lines, - if (length(init_lines) > 0) "" else NULL, + if (length(dsl_text(init_lines)) > 0) dsl_lines("") else NULL, eqn$dx, - "", + dsl_lines(""), out$out ) - paste(lines, collapse = "\n") + dsl_render_result(body) +} + +# The shared return shape of a render: text plus the origin of each line. +dsl_render_result <- function(body) { + lines <- dsl_text(body) + list(text = paste(lines, collapse = "\n"), lines = lines, src = dsl_src(body)) } # Assemble DSL text for an analytical (library-structure) model. -dsl_analytical <- function(model, header, derived, parameters) { +dsl_analytical <- function(model, header, derived, parameters, ctx) { arg_list <- model$arg_list template <- model$arg_list$tem if (is.null(template)) { @@ -667,8 +1361,26 @@ dsl_analytical <- function(model, header, derived, parameters) { )) } - out <- dsl_out_block(arg_list$out) - derived <- c(derived, out$derived) + out <- dsl_out_block(arg_list$out, ctx) + derived <- dsl_c(derived, out$derived) + + has_absorption <- stringr::str_detect(structure, "absorption") + if (!has_absorption && (!is.null(arg_list$lag) || !is.null(arg_list$fa))) { + cli::cli_abort("The `lag` and `fa` blocks can only be used with analytical bolus models.") + } + + lag_lines <- dsl_lines() + fa_lines <- dsl_lines() + if (has_absorption && !is.null(arg_list$lag)) { + lag <- dsl_route_property_block(arg_list$lag, "lag", ctx) + derived <- dsl_c(derived, lag$derived) + lag_lines <- lag$lines + } + if (has_absorption && !is.null(arg_list$fa)) { + fa <- dsl_route_property_block(arg_list$fa, "fa", ctx) + derived <- dsl_c(derived, fa$derived) + fa_lines <- fa$lines + } # The DSL analytical structures require specific derived-parameter names (e.g. # `kcp`, `kpc`, `vc`). The Pmetrics model-library templates use their own @@ -677,7 +1389,7 @@ dsl_analytical <- function(model, header, derived, parameters) { param_aliases <- dsl_analytical_param_map(structure) if (length(param_aliases) > 0) { alias_lines <- paste0(names(param_aliases), " = ", unname(param_aliases)) - derived <- c(alias_lines, derived) + derived <- dsl_c(dsl_lines(alias_lines), derived) } # Determine the number of compartments the structure requires. @@ -689,26 +1401,29 @@ dsl_analytical <- function(model, header, derived, parameters) { # Declare the dose route. Absorption ("bolus") templates receive a bolus into # the depot (x1); IV templates receive an infusion into the central # compartment (x1). - route_line <- if (stringr::str_detect(structure, "absorption")) { + route_line <- if (has_absorption) { sprintf("bolus(%s) -> x1", pm_input_label(1)) } else { sprintf("infusion(%s) -> x1", pm_input_label(1)) } - lines <- c( - header, - sprintf("structure = %s", structure), - sprintf("states = %s", paste(states, collapse = ", ")), - sprintf("outputs = %s", paste(outputs, collapse = ", ")), - "", - route_line, - "", + body <- dsl_c( + dsl_lines(header, "header"), + dsl_lines(sprintf("structure = %s", structure)), + dsl_lines(sprintf("states = %s", paste(states, collapse = ", "))), + dsl_lines(sprintf("outputs = %s", paste(outputs, collapse = ", "))), + dsl_lines(""), + dsl_lines(route_line), + dsl_lines(""), derived, - if (length(derived) > 0) "" else NULL, + if (length(dsl_text(derived)) > 0) dsl_lines("") else NULL, + lag_lines, + fa_lines, + if (length(dsl_text(lag_lines)) > 0 || length(dsl_text(fa_lines)) > 0) dsl_lines("") else NULL, out$out ) - paste(lines, collapse = "\n") + dsl_render_result(body) } # Number of states (compartments) for a DSL analytical structure. diff --git a/data/NPex.rda b/data/NPex.rda index be0103957..b1e8482e2 100755 Binary files a/data/NPex.rda and b/data/NPex.rda differ diff --git a/data/badData.rda b/data/badData.rda index 39d67dbca..79b6f7ff1 100755 Binary files a/data/badData.rda and b/data/badData.rda differ diff --git a/data/dataEx.rda b/data/dataEx.rda index 0fe3decaa..6da32b7de 100755 Binary files a/data/dataEx.rda and b/data/dataEx.rda differ diff --git a/data/modEx.rda b/data/modEx.rda index 98583c3be..1aaabcc6a 100755 Binary files a/data/modEx.rda and b/data/modEx.rda differ diff --git a/data/simEx.rda b/data/simEx.rda index de2f20c36..97a54a8a1 100755 Binary files a/data/simEx.rda and b/data/simEx.rda differ diff --git a/man/PM_model.Rd b/man/PM_model.Rd index b944b7ba0..56cafa699 100755 --- a/man/PM_model.Rd +++ b/man/PM_model.Rd @@ -301,20 +301,31 @@ Use the following notation in equations: \item \code{dx[i]} for the change in amount with respect to time (i.e., \eqn{dx/dt}), where \code{i} is the compartment number, \item \code{x[i]} for the compartment amount, where \code{i} is the compartment number. -\item \code{rateiv[j]} for the infusion rate of input \code{j}, where \code{j} is the input number -in the data corresponding to doses for that input. -\item Bolus doses are indicated by \code{DUR = 0} for dose events in the -data. Currently only one bolus input is allowed, which goes into compartment 1 -and is not modifiable. It does not appear in the differential equations. +\item \code{rateiv[j]} or \code{R[j]} for the infusion rate of input \code{j}, where \code{j} is +the input number in the data. An infusion may be a standalone additive term +or one exact product with a scale, such as \code{R[1] * iv_scale}. The scale is +evaluated continuously in the differential equation. +\item \code{B[j]} or \code{bolus[j]} for a bolus input. A bolus may be a standalone +additive term or one exact product with a scale, such as \code{B[1] * v}. +The scale is applied to the bolus dose and does not affect an infusion with +the same input number. For example, \if{html}{\out{
}}\preformatted{eqn = function() \{ - dx[1] = -ka * x[1] - dx[2] = rateiv[1] + ka * x[1] - ke * x[2] + dx[1] = B[1] * fa1 - ka * x[1] + dx[2] = rateiv[1] * iv_scale + ka * x[1] - ke * x[2] \} }\if{html}{\out{
}} } +\item \strong{Expressions} are ordinary R arithmetic. Conditions may appear anywhere, +not only as a whole equation: \code{2 * if (wt > 70) v1 else v2} and +\code{if (a > 0) (if (b > 0) 1 else 2) else 3} both work. \code{ifelse(cond, a, b)}, +\code{pmin}, \code{pmax}, and \code{max}/\code{min} with any number of arguments are accepted; +\code{log(x, base)}, \code{trunc}, \code{sign} and \code{round(x, digits)} follow R's rules +(including R's round-half-to-even). Functions R provides but the DSL does +not - such as \code{pnorm}, \code{log1p} or \code{gamma} - are reported when the model is +defined, naming the call and the supported alternatives. \item \strong{Additional equations} in R code can be defined in this block, which are similar to the \code{sec} block, but will only be available within the \code{eqn} block as opposed to global availability when defined in \code{sec}. They can be added to either @@ -338,7 +349,23 @@ this code could be used to model delayed absorption if an antacid is present. As for \code{eqn}, additional equations in R code can be defined in this block, but will only be available within the \code{lag} block.} \item{\code{fa}}{A function defining the bioavailability (fraction absorbed) equations, -similar to \code{lag}. +similar to \code{lag}. The \code{lag} and \code{fa} blocks also apply to analytical bolus +models; they are not valid for analytical infusion-only models. + +Effect functions for combination drug action are also available, anywhere an +expression is allowed. They compute an effect from supplied coefficients and +do not fit those coefficients: +\itemize{ +\item \code{estimate_effect_2(u, v, alpha, h1, h2)} for two sites, and +\item \code{estimate_effect_3(a, b, c, alpha12, alpha13, alpha23, alpha123, h1, h2, h3)} +for three. +} + +\code{u} and \code{v} are normalised exposures (concentration divided by the site's +\code{e50}), \code{alpha} is the interaction coefficient and \code{h1}/\code{h2} are the Hill +exponents. The interaction term is computed internally as +\eqn{w = \alpha \cdot u \cdot v}; the retired six-argument \code{get_e2(a, b, w, h1, h2, alpha_s)} call must be rewritten as \code{estimate_effect_2(a, b, alpha_s, h1, h2)}. These functions are evaluated at run time, on every derivative or +output evaluation, and cannot be used in the \code{pri} block. Example: @@ -362,7 +389,8 @@ Example: This sets the initial amount of drug in compartment 2 to the value of a covariate \code{init2} multiplied by the volume of the compartment, \code{V}, assuming \code{V} is either a primary parameter or defined in the -\code{sec} block. +\code{sec} block. Initial conditions are evaluated before the first event of +every occasion, including occasions started by an \code{EVID = 4} row. As for \code{eqn}, additional equations in R code can be defined in this block, but will only be available within the \code{ini} block.} diff --git a/src/rust/Cargo.toml b/src/rust/Cargo.toml index 43df0b1bb..fa592aab0 100755 --- a/src/rust/Cargo.toml +++ b/src/rust/Cargo.toml @@ -9,7 +9,7 @@ name = 'pm_rs' [dependencies] extendr-api = "=0.9.0" -pmcore = { version = "=0.28.0", features = ["dsl"] } +pmcore = { version = "=0.29.0", features = ["dsl"] } rayon = "1.12.0" anyhow = "1.0.103" tracing = "0.1.44" diff --git a/tests/testthat/test-data-currency.R b/tests/testthat/test-data-currency.R index 0af10ad16..8bc64149c 100644 --- a/tests/testthat/test-data-currency.R +++ b/tests/testthat/test-data-currency.R @@ -97,6 +97,13 @@ test_that("shipped example objects expose the current R6 members", { }) test_that("shipped example objects embed the current R6 method bodies", { + # covr::package_coverage() instruments every R6 method body in place + # (covr:::traverse_R6), so the live generator can never match the + # uninstrumented bodies frozen in data/*.rda. The member-set test above + # is unaffected because instrumentation preserves method names, and this + # check still runs for real under R CMD check. + skip_on_covr() + objects <- shipped_r6_objects() problems <- character() diff --git a/tests/testthat/test-io-labels.R b/tests/testthat/test-io-labels.R index 4cfc5b1f2..76e8b5f83 100644 --- a/tests/testthat/test-io-labels.R +++ b/tests/testthat/test-io-labels.R @@ -60,8 +60,9 @@ testthat::test_that("An input driving both a bolus and an infusion keeps one lab ) mod$compile(quiet = TRUE) - testthat::expect_match(mod$dsl, "bolus(input_1) -> x1", fixed = TRUE) - testthat::expect_match(mod$dsl, "infusion(input_1) -> x1", fixed = TRUE) + testthat::expect_match(mod$dsl, "bolus(input_1)", fixed = TRUE) + testthat::expect_match(mod$dsl, "infusion(input_1)", fixed = TRUE) + testthat::expect_false(grepl("->", mod$dsl, fixed = TRUE)) testthat::expect_false(grepl("input_2", mod$dsl, fixed = TRUE)) # The shared label must be accepted by the pharmsol backend. testthat::expect_no_error(model_parameters(mod$dsl)) diff --git a/tests/testthat/test-model-dsl-dose-scaling.R b/tests/testthat/test-model-dsl-dose-scaling.R new file mode 100644 index 000000000..9b647751f --- /dev/null +++ b/tests/testthat/test-model-dsl-dose-scaling.R @@ -0,0 +1,155 @@ +testthat::skip_on_cran() + +model_to_dsl <- getFromNamespace("model_to_dsl", "Pmetrics") + +scaled_bolus_model <- function(eqn_block, fa_block = NULL) { + PM_model$new( + list( + pri = list( + kes = ab(0, 5), + vs = ab(0, 5), + ff = ab(0, 1) + ), + cov = list(wt = interp()), + sec = function() { + v <- vs * (wt / 70) + ke <- kes * (wt / 70)^-0.25 + }, + fa = fa_block, + eqn = eqn_block, + out = function() { + y[1] <- x[1] / v + y[2] <- x[2] + }, + err = list( + additive(1, c(1, 0, 0, 0), outeq = 1), + additive(1, c(1, 0, 0, 0), outeq = 2) + ) + ), + compile = FALSE + ) +} + +testthat::test_that("ODE dose inputs and derived scales stay on the RHS", { + mod <- scaled_bolus_model(function() { + dx[1] <- B[1] * v + R[1] - ke * x[1] + dx[2] <- B[2] + x[1] / v * ff + }) + + dsl <- model_to_dsl(mod) + testthat::expect_match(dsl, "bolus(input_1) * (v)", fixed = TRUE) + testthat::expect_match(dsl, "infusion(input_1)", fixed = TRUE) + testthat::expect_false(grepl("->", dsl, fixed = TRUE)) + testthat::expect_false(grepl("fa(input_1)", dsl, fixed = TRUE)) + testthat::expect_match(dsl, "v = vs * (wt / 70.0)", fixed = TRUE) +}) + +testthat::test_that("The bolus may be on either side of the scale", { + mod <- scaled_bolus_model(function() { + dx[1] <- v * B[1] + R[1] - ke * x[1] + dx[2] <- B[2] + x[1] / v * ff + }) + + testthat::expect_match(model_to_dsl(mod), "bolus(input_1) * (v)", fixed = TRUE) +}) + +testthat::test_that("Equation scaling and an explicit fa block are multiplied", { + mod <- scaled_bolus_model( + eqn_block = function() { + dx[1] <- B[1] * v + R[1] - ke * x[1] + dx[2] <- B[2] + x[1] / v * ff + }, + fa_block = function() { + fa[1] <- ff + } + ) + + dsl <- model_to_dsl(mod) + testthat::expect_match(dsl, "bolus(input_1) * (ff) * (v)", fixed = TRUE) + testthat::expect_false(grepl("fa(input_1)", dsl, fixed = TRUE)) + matches <- gregexpr("bolus(input_1)", dsl, fixed = TRUE)[[1]] + testthat::expect_equal(sum(matches > 0), 1) +}) + +testthat::test_that("A linear infusion scale stays explicit in the derivative", { + right <- scaled_bolus_model(function() { + dx[1] <- R[1] * v - ke * x[1] + dx[2] <- B[2] + x[1] / v * ff + }) + left <- scaled_bolus_model(function() { + dx[1] <- v * R[1] - ke * x[1] + dx[2] <- B[2] + x[1] / v * ff + }) + + for (mod in list(right, left)) { + dsl <- model_to_dsl(mod) + testthat::expect_false(grepl("->", dsl, fixed = TRUE)) + testthat::expect_match(dsl, "infusion(input_1) * (v)", fixed = TRUE) + testthat::expect_false(grepl("fa(input_1)", dsl, fixed = TRUE)) + } +}) + +testthat::test_that("Complex route expressions stay unsupported", { + nested_infusion <- scaled_bolus_model(function() { + dx[1] <- R[1] * v * ff - ke * x[1] + dx[2] <- B[2] + x[1] / v * ff + }) + testthat::expect_error(model_to_dsl(nested_infusion), "Infusion inputs may only") + + nested_bolus <- scaled_bolus_model(function() { + dx[1] <- B[1] * v * ff + R[1] - ke * x[1] + dx[2] <- B[2] + x[1] / v * ff + }) + testthat::expect_error(model_to_dsl(nested_bolus), "Bolus inputs may only") + + state_dependent <- scaled_bolus_model(function() { + dx[1] <- B[1] * x[1] + R[1] - ke * x[1] + dx[2] <- B[2] + x[1] / v * ff + }) + testthat::expect_error(model_to_dsl(state_dependent), "cannot depend directly on state") +}) + +testthat::test_that("Analytical bolus models emit lag and fa", { + mod <- PM_model$new( + pri = list(Ka = ab(0.5, 6), Ke = ab(0.1, 1.5), V = ab(25, 120)), + eqn = function() { + one_comp_bolus + }, + lag = function() { + lag[1] <- 2 + }, + fa = function() { + fa[1] <- 0.5 + }, + out = function() { + y[1] <- x[2] / V + }, + err = list(proportional(5, c(0.1, 0.15, 0, 0))), + compile = FALSE + ) + + dsl <- model_to_dsl(mod) + testthat::expect_match(dsl, "bolus(input_1) -> x1", fixed = TRUE) + testthat::expect_match(dsl, "lag(input_1) = 2.0", fixed = TRUE) + testthat::expect_match(dsl, "fa(input_1) = 0.5", fixed = TRUE) + testthat::expect_no_error(model_parameters(dsl)) +}) + +testthat::test_that("Analytical infusion models reject lag and fa", { + mod <- PM_model$new( + pri = list(Ke = ab(0.1, 1.5), V = ab(25, 120)), + eqn = function() { + one_comp_iv + }, + fa = function() { + fa[1] <- 0.5 + }, + out = function() { + y[1] <- x[1] / V + }, + err = list(proportional(5, c(0.1, 0.15, 0, 0))), + compile = FALSE + ) + + testthat::expect_error(model_to_dsl(mod), "analytical bolus models") +}) diff --git a/tests/testthat/test-model-dsl-language.R b/tests/testthat/test-model-dsl-language.R new file mode 100644 index 000000000..2571c4d1b --- /dev/null +++ b/tests/testthat/test-model-dsl-language.R @@ -0,0 +1,480 @@ +testthat::skip_on_cran() + +# The R-to-DSL emitter should read like R: constructs R users expect are either +# translated exactly or rejected with a message that names the fix. These tests +# pin both halves of that contract, plus the numbers the translated code +# actually produces. + +model_to_dsl <- getFromNamespace("model_to_dsl", "Pmetrics") +model_to_dsl_traced <- getFromNamespace("model_to_dsl_traced", "Pmetrics") +dsl_is_valid <- function(dsl) { + is.null(tryCatch({ model_metadata(dsl, NULL); NULL }, error = function(e) e)) +} + +make_fn <- function(text) { + eval(call("function", as.pairlist(alist()), parse(text = text)[[1]])) +} + +# A model whose single output is the expression under test. State 1 is frozen so +# the output is a pure function of the parameters; state 2 exists so the dose +# events in the template data resolve to a declared route. +dose_model <- function(out_expr, x_init = 3.5, sec_lines = NULL, ini_lines = NULL, cov_names = NULL, compile = FALSE) { + fn <- function(lines) make_fn(sprintf("{\n%s\n}", paste(lines, collapse = "\n"))) + PM_model$new( + pri = list(p1 = ab(-10, 10), p2 = ab(-10, 10), p3 = ab(-10, 10)), + cov = if (is.null(cov_names)) NULL else stats::setNames(lapply(cov_names, function(x) interp()), cov_names), + sec = if (is.null(sec_lines)) NULL else fn(sec_lines), + ini = fn(if (is.null(ini_lines)) c(sprintf("x[1] <- %s", format(x_init, digits = 15)), "x[2] <- 0.0") else ini_lines), + eqn = fn(c("dx[1] <- 0.0", "dx[2] <- b[1]")), + out = fn(sprintf("y[1] <- %s", out_expr)), + err = list(additive(1, c(0, 0.1, 0, 0))), + compile = compile + ) +} + +# Run the expression through the compiled DSL and compare with R's own answer. +dsl_value <- function(out_expr, p1 = 1, p2 = 1, p3 = 1, x_init = 3.5, sec_lines = NULL, ini_lines = NULL, cov_names = NULL) { + mod <- dose_model(out_expr, x_init = x_init, sec_lines = sec_lines, ini_lines = ini_lines, cov_names = cov_names) + res <- mod$sim(data = dataEx, theta = matrix(c(p1, p2, p3), nrow = 1), quiet = TRUE) + res$out[1] +} + +r_value <- function(out_expr, p1 = 1, p2 = 1, p3 = 1, x_init = 3.5, cov_values = list(), sec = NULL) { + x <- c(x_init) + env <- list2env(c(list(p1 = p1, p2 = p2, p3 = p3, x = x), cov_values), parent = baseenv()) + if (!is.null(sec)) { + # Evaluate the secondary block in R so the reference value is R's answer for + # the whole model, not a hand-computed stand-in. + eval(parse(text = paste(c("{", sec, "}"), collapse = "\n")), envir = env) + } + eval(parse(text = out_expr)[[1]], envir = env) +} + +testthat::test_that("an unsupported function is reported against the R expression", { + for (bad in c("pnorm(p1)", "qnorm(p1)", "foo(p1, p2)")) { + message <- tryCatch(model_to_dsl(dose_model(bad)), error = function(e) conditionMessage(e)) + testthat::expect_match(message, "Unsupported function") + testthat::expect_match(message, "Supported functions are") + testthat::expect_false(grepl("subscript out of bounds", message, fixed = TRUE)) + } +}) + +testthat::test_that("values that are not finite numbers are rejected with a fix", { + for (bad in c("Inf * p1", "NaN + p1", "NA * p1")) { + message <- tryCatch(model_to_dsl(dose_model(bad)), error = function(e) conditionMessage(e)) + testthat::expect_match(message, "finite numbers|Missing values") + testthat::expect_false(grepl("subscript out of bounds", message, fixed = TRUE)) + } +}) + +testthat::test_that("the retired effect name is reported as a migration", { + message <- tryCatch( + model_to_dsl(dose_model("get_e2(p1, p2, p3, 1, 1, 0.5)")), + error = function(e) conditionMessage(e) + ) + testthat::expect_match(message, "estimate_effect_2") + testthat::expect_match(message, "w = alpha \\* u \\* v") +}) + +testthat::test_that("effect functions are available with a checked signature", { + two <- model_to_dsl(dose_model("estimate_effect_2(p1, p2, p3, 1, 1)")) + testthat::expect_match(two, "out(outeq_1) = estimate_effect_2(p1, p2, p3, 1.0, 1.0)", fixed = TRUE) + testthat::expect_true(dsl_is_valid(two)) + + three <- model_to_dsl(dose_model("estimate_effect_3(p1, p2, p3, 0.1, 0.1, 0.1, 0.1, 1, 1, 1)")) + testthat::expect_match(three, "estimate_effect_3(", fixed = TRUE) + testthat::expect_true(dsl_is_valid(three)) + + message <- tryCatch( + model_to_dsl(dose_model("estimate_effect_2(p1, p2)")), + error = function(e) conditionMessage(e) + ) + testthat::expect_match(message, "takes 5 arguments") + testthat::expect_match(message, "estimate_effect_2\\(u, v, alpha, h1, h2\\)") +}) + +testthat::test_that("conditionals are accepted in any position", { + # The DSL can only put a conditional in a whole right-hand side; everything + # else is hoisted into a derived helper. + cases <- c( + "2.0 * if (p1 > 0) p2 else p3", + "if (p1 > 0) (if (p2 > 0) 1 else 2) else 3", + "if ((if (p1 > 0) 1 else 2) > 1.5) p2 else p3", + "exp(if (p1 > 0) 1 else 2)", + "p1 + 2.0 * if (p2 > 0) 1 else 3", + "-(2.0 * if (p1 > 0) p2 else p3)" + ) + for (case in cases) { + dsl <- model_to_dsl(dose_model(case)) + testthat::expect_true(dsl_is_valid(dsl), info = case) + } +}) + +testthat::test_that("conditional dose inputs are refused with the supported alternative", { + mod <- PM_model$new( + pri = list(p1 = ab(0, 1), p2 = ab(0, 1), p3 = ab(0, 1)), + eqn = make_fn("{\ndx[1] <- if (p1 > 0.0) b[1] else 0.0\ndx[2] <- 0.0\n}"), + out = make_fn("y[1] <- x[1]"), + err = list(additive(1, c(0, 0.1, 0, 0))), + compile = FALSE + ) + message <- tryCatch({ mod$compile(quiet = TRUE); "no error" }, error = function(e) conditionMessage(e)) + testthat::expect_match(message, "dose input cannot be inside a conditional") + testthat::expect_match(message, "bolus\\(input\\) \\* scale") +}) + +testthat::test_that("the function names R has that the DSL does not are translated", { + translated <- list( + "pmin(p1, p2)" = "min(p1, p2)", + "pmax(p1, p2, p3)" = "max(max(p1, p2), p3)", + "max(p1, p2, p3)" = "max(max(p1, p2), p3)", + "min(p1)" = "p1", + "log(p1, 2)" = "(ln(p1) / ln(2.0))" + ) + for (case in names(translated)) { + dsl <- model_to_dsl(dose_model(case)) + testthat::expect_true(grepl(translated[[case]], dsl, fixed = TRUE), info = case) + testthat::expect_true(dsl_is_valid(dsl), info = case) + } + + # trunc and sign need a helper, so check the pieces rather than one string. + trunc_dsl <- model_to_dsl(dose_model("trunc(p1)")) + testthat::expect_match(trunc_dsl, "ceil(", fixed = TRUE) + testthat::expect_match(trunc_dsl, "floor(", fixed = TRUE) + testthat::expect_true(dsl_is_valid(trunc_dsl)) + + sign_dsl <- model_to_dsl(dose_model("sign(p1)")) + testthat::expect_match(sign_dsl, "if (", fixed = TRUE) + testthat::expect_true(dsl_is_valid(sign_dsl)) + + # R rounds halves to even; the DSL's round does not, so it is emitted as R's + # rule and not as the intrinsic. + round_dsl <- model_to_dsl(dose_model("round(p1)")) + testthat::expect_false(grepl("round(", round_dsl, fixed = TRUE)) + testthat::expect_true(dsl_is_valid(round_dsl)) +}) + +testthat::test_that("translated expressions compute what R computes", { + # One entry per construct. Anything the model must declare beyond the three + # parameters (a covariate, a secondary value) is passed alongside, so a case + # is a complete model rather than only an expression. + case <- function(expr, p1 = 1, p2 = 1, p3 = 1, x = 3.5, + cov_names = NULL, cov_values = list(), sec = NULL) { + list(expr = expr, p1 = p1, p2 = p2, p3 = p3, x = x, + cov_names = cov_names, cov_values = cov_values, sec = sec) + } + + cases <- list( + # --- R functions the DSL does not have ------------------------------- + case("pmin(p1, p2)", p1 = 2, p2 = 5), + case("pmin(p1, p2)", p1 = -3, p2 = -7), + case("pmax(p1, p2)", p1 = -3, p2 = -7), + case("pmax(p1, p2, p3)", p1 = 2, p2 = 5, p3 = 7.5), + case("max(p1, p2, p3)", p1 = 2, p2 = 5, p3 = 7.5), + case("min(p1, p2)", p1 = -2, p2 = 5), + case("max(p1)", p1 = 7), + case("min(p1, p2, p3, 0)", p1 = 1, p2 = 2, p3 = 3), + case("log(p1, 2)", p1 = 8), + case("log(p1, 10)", p1 = 1000), + case("log(p1)", p1 = exp(3)), + case("exp(log(p1, 2))", p1 = 8), + case("trunc(p1 * x[1])", p1 = -1.7), + case("trunc(p1 * x[1])", p1 = 1.7), + case("trunc(p1)", p1 = -0.5), + case("trunc(p1)", p1 = 0.5), + case("sign(p1)", p1 = -2), + case("sign(p1)", p1 = 0), + case("sign(p1)", p1 = 2), + case("sign(p1)", p1 = -0.0001), + # --- round: R rounds halves to even, the DSL intrinsic does not ------ + case("round(x[1])", x = 2.5), + case("round(x[1])", x = 3.5), + case("round(x[1])", x = -2.5), + case("round(x[1])", x = -3.5), + case("round(x[1], 0)", x = 2.5), + case("round(p1 * x[1], 2)", p1 = 1.234, x = 3), + case("round(p1 * x[1], 1)", p1 = 0.25, x = 5), + case("round(x[1], -1)", x = 25), + # --- conditionals ---------------------------------------------------- + case("if (p1 > 0) 1 else 2", p1 = 1), + case("if (p1 > 0) 1 else 2", p1 = -1), + case("if (p1 > 0) 2 else 4", p1 = 1), + case("2.0 * if (p1 > 0) p2 else p3", p1 = 1, p2 = 3, p3 = 9), + case("2.0 * if (p1 > 0) p2 else p3", p1 = -1, p2 = 3, p3 = 9), + case("2.0 * if (p1 > 0) 1 else 2", p1 = 1), + case("2.0 * if (p1 > 0) 1 else 2", p1 = -1), + case("p1 / if (p2 > 0) 2 else 4", p1 = 8, p2 = 1), + case("ifelse(p1 > 0, p2, p3)", p1 = -1, p2 = 3, p3 = 9), + case("if (p1 > 0) (if (p2 > 0) p2 else p3) else p1", p1 = 1, p2 = -2, p3 = 9), + case("if (p1 > 0) p1 else if (p2 > 0) p2 else p3", p1 = 0, p2 = 5, p3 = 9), + case("if (p1 > 0) p1 else if (p2 > 0) p2 else p3", p1 = 0, p2 = 0, p3 = 9), + case("(if (p1 > 0) 1 else 2) * (if (p2 > 0) 3 else 4)", p1 = 1, p2 = 1), + case("(if (p1 > 0) 1 else 2) * (if (p2 > 0) 3 else 4)", p1 = -1, p2 = 1), + case("if (p1 > 0) 1 else 2 + if (p2 > 0) 3 else 4", p1 = 1, p2 = 1), + case("if (p1 > 0) 1 else 2 + if (p2 > 0) 3 else 4", p1 = -1, p2 = 1), + case("exp(if (p1 > 0) 1 else 2)", p1 = 1), + case("exp(if (p1 > 0) 1 else 2)", p1 = -1), + case("if ((if (p1 > 0) 1 else 2) > 1.5) p2 else p3", p1 = -1, p2 = 3, p3 = 9), + case("p1 + 2.0 * if (p2 > 0) 1 else 3", p1 = 1, p2 = -1), + case("-(2.0 * if (p1 > 0) p2 else p3)", p1 = 1, p2 = 3, p3 = 9), + # --- conditionals on a covariate, and a helper-name clash ------------ + case("if (wt > 40) p1 else p2", p1 = 7, p2 = 9, cov_names = "wt", cov_values = list(wt = 46.7)), + case("if (wt > 50) p1 else p2", p1 = 7, p2 = 9, cov_names = "wt", cov_values = list(wt = 46.7)), + case("z1 + z1", sec = "z1 <- if (p1 > 0) 2 else 3", p1 = 1), + # A helper name can only collide through a backticked assignment, which the + # emitter must notice when it picks names for generated helpers. + case("z1 + `_pmt1`", sec = c("`_pmt1` <- 100", "z1 <- 2 * if (p1 > 0) 1 else 2"), p1 = 1) + ) + + for (cs in cases) { + expected <- r_value(cs$expr, cs$p1, cs$p2, cs$p3, cs$x, cs$cov_values, cs$sec) + got <- dsl_value(cs$expr, cs$p1, cs$p2, cs$p3, cs$x, + sec_lines = cs$sec, cov_names = cs$cov_names + ) + testthat::expect_equal(got, expected, tolerance = 1e-9, + info = sprintf( + "%s with p=(%s, %s, %s), x=%s%s", cs$expr, cs$p1, cs$p2, cs$p3, cs$x, + if (is.null(cs$sec)) "" else paste0(", sec=", paste(cs$sec, collapse = "; ")) + ) + ) + } +}) + +testthat::test_that("derived values keep their value through the DSL", { + # A derived value used by an output must survive the DSL's f64 derived buffer, + # whether it comes from a parameter, a constant or a conditional. + testthat::expect_equal(dsl_value("z1", sec_lines = "z1 <- 1.0"), 1) + testthat::expect_equal(dsl_value("z1", sec_lines = "z1 <- 1.5"), 1.5) + testthat::expect_equal(dsl_value("z1", sec_lines = "z1 <- p1", p1 = 4.25), 4.25) + testthat::expect_equal(dsl_value("z1", sec_lines = "z1 <- if (p1 > 0) 1.0 else 2.0", p1 = 1), 1) +}) + +testthat::test_that("a backend diagnostic names the R statement that produced it", { + traced <- model_to_dsl_traced(dose_model("z1 * undeclared_value")) + testthat::expect_equal(length(traced$lines), length(traced$src)) + testthat::expect_match(traced$src[[which(grepl("^out\\(", traced$lines))]], "output block: y\\[1\\]") + + mod <- dose_model("undeclared_value * x[1]") + message <- tryCatch({ mod$compile(quiet = TRUE); "no error" }, error = function(e) conditionMessage(e)) + testthat::expect_match(message, "unknown identifier") + testthat::expect_match(message, "output block: y\\[1\\]") +}) + +testthat::test_that("a model that cannot compile fails where it is defined", { + # The constructor validates by default, so an undeclared name is reported when + # the model is created rather than at the start of a fit. + message <- tryCatch( + { + dose_model("undeclared_value * p1", compile = TRUE) + "no error" + }, + error = function(e) conditionMessage(e) + ) + testthat::expect_match(message, "unknown identifier") + + # A model built with compile = FALSE is validated when it is compiled. + mod <- dose_model("undeclared_value * p1") + testthat::expect_null(mod$dsl) + message <- tryCatch( + { + mod$compile(quiet = TRUE) + "no error" + }, + error = function(e) conditionMessage(e) + ) + testthat::expect_match(message, "unknown identifier") +}) + +testthat::test_that("every rejection is a structured error that names the fix", { + # A model that cannot be translated must fail as a cli condition with an + # explanation. Raw R errors used to escape - every unknown function raised + # `subscript out of bounds`, because the lookup that produced the helpful + # message could never be reached. + cases <- list( + list("pnorm(p1)", "Unsupported function"), + list("qnorm(p1)", "Unsupported function"), + list("dnorm(p1)", "Unsupported function"), + list("gamma(p1)", "Unsupported function"), + list("lgamma(p1)", "Unsupported function"), + list("log1p(p1)", "Unsupported function"), + list("expm1(p1)", "Unsupported function"), + list("c(1, 2)", "Unsupported function"), + list("base::exp(p1)", "cannot be translated"), + list("get_e2(p1, p2, p3, 1, 1, 0.5)", "estimate_effect_2"), + list("get_e3(p1, p2, p3, 1, 1, 1, 1, 1, 1, 1, 1)", "estimate_effect_3"), + list("estimate_effect_2(p1, p2)", "takes 5 arguments"), + list("estimate_effect_2(p1, p2, p3, 1, 1, 1)", "takes 5 arguments"), + list("Inf * p1", "finite numbers"), + list("NaN + p1", "finite numbers"), + list("NA * p1", "Missing values"), + list("if (p1 > 0) 1", "else branch"), + list("\"text\"", "token"), + list("p1[1]", "cannot be used as a value"), + list("y[1] * 2", "cannot be used as a value"), + list("dx[1] * 2", "cannot be used as a value"), + list("x[i]", "Dynamic") + ) + + for (cs in cases) { + err <- tryCatch( + { + model_to_dsl(dose_model(cs[[1]])) + NULL + }, + error = function(e) e + ) + testthat::expect_false(is.null(err), info = cs[[1]]) + testthat::expect_true(inherits(err, "rlang_error"), info = cs[[1]]) + testthat::expect_match(conditionMessage(err), cs[[2]], info = cs[[1]]) + } + + # Control flow is refused, and the message says what to write instead. + err <- tryCatch( + { + model_to_dsl(dose_model("z1", sec_lines = "for (i in 1:3) z1 <- i")) + NULL + }, + error = function(e) e + ) + testthat::expect_true(inherits(err, "rlang_error")) + testthat::expect_match(conditionMessage(err), "Only assignments are supported") +}) + +testthat::test_that("a braced conditional block chooses between scalar values", { + dsl <- model_to_dsl( + dose_model( + "x[1] + z2", + sec_lines = c( + "if (p1 > 0) {", + " z1 <- 1", + " z2 <- 2", + "} else {", + " z1 <- 3", + " z2 <- 4", + "}" + ) + ) + ) + testthat::expect_match(dsl, "if (p1 > 0.0) {", fixed = TRUE) + testthat::expect_match(dsl, "} else {", fixed = TRUE) + testthat::expect_match(dsl, "z1 = 1.0", fixed = TRUE) + testthat::expect_match(dsl, "z2 = 4.0", fixed = TRUE) + testthat::expect_true(dsl_is_valid(dsl)) + + # The chosen value reaches the arithmetic: p1 > 0 selects the first branch. + testthat::expect_equal( + dsl_value( + "x[1] + z2", + p1 = 1, + sec_lines = c("if (p1 > 0) {", " z1 <- 1", " z2 <- 2", "} else {", " z1 <- 3", " z2 <- 4", "}") + ), + 3.5 + 2 + ) + testthat::expect_equal( + dsl_value( + "x[1] + z2", + p1 = -1, + sec_lines = c("if (p1 > 0) {", " z1 <- 1", " z2 <- 2", "} else {", " z1 <- 3", " z2 <- 4", "}") + ), + 3.5 + 4 + ) +}) + +testthat::test_that("a conditional block without an else keeps the earlier value", { + dsl <- model_to_dsl( + dose_model("z1", sec_lines = c("z1 <- 0", "if (p1 > 0) {", " z1 <- 1", "}")) + ) + testthat::expect_match(dsl, "z1 = 0.0", fixed = TRUE) + testthat::expect_match(dsl, "if (p1 > 0.0) {", fixed = TRUE) + testthat::expect_true(dsl_is_valid(dsl)) + + testthat::expect_equal( + dsl_value("z1", p1 = 1, sec_lines = c("z1 <- 0", "if (p1 > 0) {", " z1 <- 1", "}")), 1 + ) + testthat::expect_equal( + dsl_value("z1", p1 = -1, sec_lines = c("z1 <- 0", "if (p1 > 0) {", " z1 <- 1", "}")), 0 + ) +}) + +testthat::test_that("a conditional block around an equation says what to write instead", { + # The DSL accepts only plain assignments inside an `if` statement, so a + # derivative, output or initial condition inside one is refused with the + # conditional-equation form. + for (case in list( + list(c("if (p1 > 0) {", " dx[1] <- -p1 * x[1]", "} else {", " dx[1] <- -p2 * x[1]", "}"), "equation block", "dx\\[1\\] <- if \\(cond\\) a else b"), + list(c("if (p1 > 0) {", " y[1] <- 1", "} else {", " y[1] <- 2", "}"), "output block", "conditional"), + list(c("if (p1 > 0) {", " x[1] <- 1", "} else {", " x[1] <- 2", "}"), "initial-conditions block", "conditional") + )) { + lines <- case[[1]] + where <- case[[2]] + mod <- if (where == "equation block") { + PM_model$new( + pri = list(p1 = ab(0, 1), p2 = ab(0, 1), p3 = ab(0, 1)), + eqn = make_fn(sprintf("{\n%s\ndx[2] <- 0.0\n}", paste(lines, collapse = "\n"))), + out = make_fn("y[1] <- x[1]"), + err = list(additive(1, c(0, 0.1, 0, 0))), compile = FALSE + ) + } else if (where == "output block") { + PM_model$new( + pri = list(p1 = ab(0, 1), p2 = ab(0, 1), p3 = ab(0, 1)), + eqn = make_fn("{\ndx[1] <- 0.0\ndx[2] <- 0.0\n}"), + out = make_fn(sprintf("{\n%s\n}", paste(lines, collapse = "\n"))), + err = list(additive(1, c(0, 0.1, 0, 0))), compile = FALSE + ) + } else { + PM_model$new( + pri = list(p1 = ab(0, 1), p2 = ab(0, 1), p3 = ab(0, 1)), + ini = make_fn(sprintf("{\n%s\n}", paste(lines, collapse = "\n"))), + eqn = make_fn("{\ndx[1] <- 0.0\ndx[2] <- 0.0\n}"), + out = make_fn("y[1] <- x[1]"), + err = list(additive(1, c(0, 0.1, 0, 0))), compile = FALSE + ) + } + err <- tryCatch( + { + model_to_dsl(mod) + NULL + }, + error = function(e) e + ) + testthat::expect_true(inherits(err, "rlang_error"), info = where) + testthat::expect_match(conditionMessage(err), where, info = where) + testthat::expect_match(conditionMessage(err), case[[3]], info = where) + } +}) + +testthat::test_that("a diagnostic names the file and line when R kept them", { + # The model's blocks are defined here, so R keeps their source references when + # the file was sourced with keep.source = TRUE. A plain `Rscript` run does not, + # and the diagnostic is still complete without them. + mod <- PM_model$new( + pri = list(p1 = ab(0, 1), p2 = ab(0, 1), p3 = ab(0, 1)), + ini = function() { + x[1] <- 3.5 + x[2] <- 0.0 + }, + eqn = function() { + dx[1] <- 0.0 + dx[2] <- b[1] + }, + out = function() { + y[1] <- undeclared_name * p1 + }, + err = list(additive(1, c(0, 0.1, 0, 0))), compile = FALSE + ) + + if (is.null(attr(body(mod$arg_list$out), "srcref"))) { + testthat::skip("source references are not kept in this environment") + } + + err <- tryCatch( + { + mod$compile(quiet = TRUE) + NULL + }, + error = function(e) e + ) + testthat::expect_true(inherits(err, "rlang_error")) + testthat::expect_match(conditionMessage(err), "unknown identifier") + testthat::expect_match(conditionMessage(err), "test-model-dsl-language\\.R:[0-9]+") +}) diff --git a/tests/testthat/test-pm-model-from-dsl.R b/tests/testthat/test-pm-model-from-dsl.R index aec9d3428..ae34b5bde 100644 --- a/tests/testthat/test-pm-model-from-dsl.R +++ b/tests/testthat/test-pm-model-from-dsl.R @@ -4,8 +4,7 @@ kind = ode params = k, v states = central outputs = cp -infusion(input_1) -> central -dx(central) = rate(input_1) - k * central +dx(central) = infusion(input_1) - k * central out(cp) = central / v ~ continuous() " @@ -36,9 +35,7 @@ kind = ode params = k states = central outputs = cp -bolus(oral) -> central -infusion(iv) -> central -dx(central) = -k * central +dx(central) = bolus(oral) + infusion(iv) - k * central out(cp) = central ~ continuous() " diff --git a/tests/testthat/test-pmcore-pharmsol-migration.R b/tests/testthat/test-pmcore-pharmsol-migration.R index 407e72331..91eaed334 100644 --- a/tests/testthat/test-pmcore-pharmsol-migration.R +++ b/tests/testthat/test-pmcore-pharmsol-migration.R @@ -23,13 +23,14 @@ build_passthrough_ode_model <- function(solver = NULL) { ) } -test_that("ODE generation produces a DSL model with 1-based indices and an infusion route", { +test_that("ODE generation produces 1-based labels and an explicit infusion term", { mod <- build_passthrough_ode_model("TSIT45") mod$compile(quiet = TRUE) dsl <- mod$dsl testthat::expect_match(dsl, "kind = ode") - testthat::expect_match(dsl, "infusion\\(input_1\\) -> x1") + testthat::expect_match(dsl, "infusion(input_1)", fixed = TRUE) + testthat::expect_false(grepl("->", dsl, fixed = TRUE)) testthat::expect_match(dsl, "dx\\(x1\\)") testthat::expect_match(dsl, "out\\(outeq_1\\)") # Outputs are 1-based to match the Pmetrics data OUTEQ column. @@ -74,57 +75,61 @@ test_that("Inline if/else expressions emit parenthesized, unbraced DSL condition ) }) -test_that("Conditionals in unsupported positions raise a clear R-level error", { - # The DSL accepts a conditional only as a whole equation right-hand side or - # chained in the `else` branch. Other positions must fail with an actionable - # R error rather than a cryptic DSL parse error on generated code. Patterns - # match single tokens so they survive cli line-wrapping of the message. +test_that("Conditionals nested in expressions are hoisted into helpers", { + # A conditional may appear anywhere in R. The DSL can only take one as a whole + # right-hand side, so the emitter hoists the value into a derived helper. expr_to_dsl <- getFromNamespace("expr_to_dsl", "Pmetrics") + dsl_ctx <- getFromNamespace("dsl_ctx", "Pmetrics") + dsl_take_hoists <- getFromNamespace("dsl_take_hoists", "Pmetrics") + dsl_text <- getFromNamespace("dsl_text", "Pmetrics") - # Nested inside an operator or function call. - testthat::expect_error( - expr_to_dsl(quote(2 * if (c) a else b)), - "right-hand" - ) - testthat::expect_error( - expr_to_dsl(quote(exp(if (c) a else b))), - "right-hand" - ) - # In the `then` branch (only the `else` branch may chain). - testthat::expect_error( - expr_to_dsl(quote(if (a) (if (b) x else y) else z)), - "right-hand" + ctx <- dsl_ctx() + testthat::expect_equal(expr_to_dsl(quote(2 * if (c) a else b), ctx), "2.0 * _pmt1") + testthat::expect_equal(dsl_text(dsl_take_hoists(ctx)), "_pmt1 = if (c) a else b") + ctx <- dsl_ctx() + testthat::expect_equal(expr_to_dsl(quote(exp(if (c) a else b)), ctx), "exp(_pmt1)") + testthat::expect_equal(dsl_text(dsl_take_hoists(ctx)), "_pmt1 = if (c) a else b") + + # A nested conditional in the `then` branch is hoisted as well. The + # parentheses are the ones the caller wrote around the inner conditional. + ctx <- dsl_ctx() + testthat::expect_equal( + expr_to_dsl(quote(if (a) (if (b) x else y) else z), ctx), + "if (a) (_pmt1) else z" ) - # Missing `else` branch. + testthat::expect_equal(dsl_text(dsl_take_hoists(ctx)), "_pmt1 = if (b) x else y") + + # A missing `else` branch is still an error. testthat::expect_error( - expr_to_dsl(quote(if (a) b)), - "include" + expr_to_dsl(quote(if (a) b), dsl_ctx()), + "else branch" ) }) -test_that("Conditionals embedded in derivative equations raise a clear error", { - # Regression: the additive-term path for `dx` equations must not bypass the - # conditional guard. An `if` summed with other terms is rejected, but a lone - # whole-RHS conditional (which the backend accepts) still emits. +test_that("Conditionals summed into derivative equations are hoisted", { + # The additive-term path for `dx` equations must hoist a conditional rather + # than emit DSL the backend rejects. dsl_eqn_block <- getFromNamespace("dsl_eqn_block", "Pmetrics") + dsl_ctx <- getFromNamespace("dsl_ctx", "Pmetrics") + dsl_text <- getFromNamespace("dsl_text", "Pmetrics") - testthat::expect_error( - dsl_eqn_block(function() { + summed <- dsl_eqn_block( + function() { dx[1] <- a + if (c) b else d - }), - "right-hand" - ) - testthat::expect_error( - dsl_eqn_block(function() { - dx[1] <- a - if (c) b else d - }), - "right-hand" + }, + dsl_ctx() ) + testthat::expect_equal(dsl_text(summed$dx), "dx(x1) = a + (_pmt1)") + testthat::expect_equal(dsl_text(summed$derived), "_pmt1 = if (c) b else d") + # A whole-RHS conditional remains valid and is emitted bare. - whole_rhs <- dsl_eqn_block(function() { - dx[1] <- if (c) a else b - }) - testthat::expect_equal(whole_rhs$dx, "dx(x1) = if (c) a else b") + whole_rhs <- dsl_eqn_block( + function() { + dx[1] <- if (c) a else b + }, + dsl_ctx() + ) + testthat::expect_equal(dsl_text(whole_rhs$dx), "dx(x1) = if (c) a else b") }) test_that("Secondary if/else conditionals generate a parseable DSL model", {