diff --git a/NAMESPACE b/NAMESPACE index 14a43f196..3c55037ad 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -324,6 +324,8 @@ S3method(vec_restore,rxEt) export("RxODE<-") export("ini<-") export("model<-") +export("rxForcedPars<-") +export("rxParLoader<-") export("rxode2<-") export("rxode<-") export(.assertRenameErrorModelLine) @@ -629,6 +631,7 @@ export(rxExpandSens3_) export(rxExpandSens_) export(rxFixPop) export(rxFixRes) +export(rxForcedPars) export(rxForget) export(rxFromSE) export(rxFun) @@ -647,6 +650,7 @@ export(rxIndLinState) export(rxIndLinStrategy) export(rxInit) export(rxInits) +export(rxInjectedPars) export(rxIntToBase) export(rxIntToLetter) export(rxInv) @@ -676,6 +680,7 @@ export(rxOldQsDes) export(rxOmegaVarCovDeriv) export(rxOmegaVarCovDeriv_) export(rxOptExpr) +export(rxParLoader) export(rxParam) export(rxParams) export(rxParseErr) @@ -692,8 +697,10 @@ export(rxProgressStop) export(rxPrune) export(rxRateDur) export(rxRawToC) +export(rxRegisterUiPrep) export(rxReload) export(rxRemoveControl) +export(rxRemoveUiPrep) export(rxRename) export(rxRepR0_) export(rxReq) diff --git a/R/err.R b/R/err.R index f7badc543..ab72f695d 100644 --- a/R/err.R +++ b/R/err.R @@ -1442,6 +1442,23 @@ rxErrTypeCombine <- function(oldErrType, newErrType) { .i <- .i + 1L } .env$iniDf <- .env$df + # A UDF modification function (rxUdfUi) can append etas to the iniDf during + # parsing (e.g. individual neural-network weight etas). .env$eta was set from + # the ini({}) omega before those UDFs ran, so refresh it from the now-final + # iniDf (the source of truth); otherwise the appended etas are mis-classified + # as covariates downstream (mu-referencing / covariate detection). + # Keep non-id (IOV) etas out: they live in .env$level, and including them in + # .env$eta makes a mu-referenced `theta + eta + iov` parse as 2 population etas. + .etaDf <- .env$df[!is.na(.env$df$neta1) & + !(.env$df$name %in% .env$level), , drop=FALSE] + if (nrow(.etaDf) > 0L) { + .etaDf <- .etaDf[order(.etaDf$neta1), , drop=FALSE] + .env$eta <- unique(.etaDf$name) + } else { + # only IOV (or no) etas remain -> no population etas; clear the list so a + # stale omega-derived name cannot leak (matches the NULL init above). + .env$eta <- NULL + } if (is.null(.env$predDf)) { ## .env$errGlobal <- c(.env$errGlobal, ## "there must be at least one prediction in the model({}) block. Use `~` for predictions") diff --git a/R/rxsolve.R b/R/rxsolve.R index d5bff77f3..071e36d7c 100644 --- a/R/rxsolve.R +++ b/R/rxsolve.R @@ -2205,6 +2205,18 @@ rxSolve.rxUi <- function(object, params = NULL, events = NULL, inits = NULL, ... .lst$drop <- c(.lst$drop, "ipredSim") } } + ## let packages rehydrate transient C state from the ui's serializable slots + ## (e.g. NN shapes) before parameters are loaded -- keeps a reloaded ui solvable. + .rxRunUiPrepHooks(object) + ## forced (externally-owned) parameters carried on the ui -- injected into every + ## solve column at setup (rxCallParLoaders), overriding params/data/inits. + if (.rxApplyForcedPars(object, .lst[[1L]])) { + on.exit(.rxClearForcedParsC(), add = TRUE) + } + ## flag this model's parameter injector (if any) so only its loader runs + if (.rxApplyParLoader(object)) { + on.exit(.rxClearActiveParLoaderC(), add = TRUE) + } .ret <- do.call("rxSolve.default", .lst) if (.pred) { .e <- attr(class(.ret), ".rxode2.env") @@ -3361,17 +3373,53 @@ rxSolve.rxSolve <- function(object, params = NULL, events = NULL, inits = NULL, # Pass the rxSolve object directly so C++ updates the correct object # rather than the global rxCurObj, which may point to a different rxSolve. if (is.null(params)) params <- object$.params.single + params <- .rxApplyInjectedPars(params, object) if (is.null(inits)) inits <- object$inits return(rxSolve.default(object, params = params, events = events, inits = inits, ..., theta = theta, eta = eta, envir = envir)) } .model <- as.character(.model) if (is.null(params)) params <- object$.params.single + params <- .rxApplyInjectedPars(params, object) if (is.null(inits)) inits <- object$inits rxSolve(.model, params = params, events = events, inits = inits, ..., theta = theta, eta = eta, envir = envir) } +## Restore par-loader-injected parameters (e.g. trained NN weights) onto the +## params for a re-solve, so re-solving from a solved object reproduces them +## even without the injecting package. Injected values override for their names +## (and are appended if absent). Handles a single named vector, a parameter +## data.frame, or a matrix; when no params are supplied it falls back to the +## solved object's stored parameter data.frame (multi-subject/multi-sim solves +## have no single named vector). +.rxApplyInjectedPars <- function(params, object) { + .inj <- rxInjectedPars(object) + if (is.null(.inj) || length(.inj) == 0L) return(params) + if (is.null(params)) params <- object$.params.dat + if (is.null(params)) return(params) + if (is.data.frame(params)) { + for (.n in names(.inj)) params[[.n]] <- .inj[[.n]] + return(params) + } + if (is.matrix(params)) { + for (.n in names(.inj)) { + if (!is.null(colnames(params)) && .n %in% colnames(params)) { + params[, .n] <- .inj[[.n]] + } else { + params <- cbind(params, + matrix(.inj[[.n]], nrow = nrow(params), ncol = 1L, + dimnames = list(NULL, .n))) + } + } + return(params) + } + if (is.numeric(params) && is.null(dim(params)) && !is.null(names(params))) { + params[names(.inj)] <- .inj + } + params +} + #' @rdname rxSolve #' @export predict.rxode2 <- function(object, ...) { @@ -3554,9 +3602,251 @@ solve.rxEt <- solve.rxSolve dadt = get(".dadt.counter", envir = .env, inherits = FALSE), jac = get(".jac.counter", envir = .env, inherits = FALSE) ), envir = .env) + ## Parameters injected by par-loaders on this solve (e.g. trained neural-network + ## weights supplied through a loader hook rather than the params vector). Saved + ## on the object so re-solving restores them even in a session without the + ## injecting package's buffer. Indices are 0-based into the model params. + .inj <- tryCatch(.Call(`_rxode2_rxGetInjectedPars`), error = function(e) NULL) + if (!is.null(.inj) && length(.inj[[1L]]) > 0L) { + .injNames <- .pars[.inj[[1L]] + 1L] + assign(".injectedPars", stats::setNames(.inj[[2L]], .injNames), envir = .env) + } else { + assign(".injectedPars", NULL, envir = .env) + } invisible(TRUE) } +#' Parameters injected into a solve by a par-loader hook +#' +#' Returns the parameters (e.g. trained neural-network weights) that a +#' registered par-loader wrote into the solve parameter vector, saved on the +#' solved object so re-solving from it restores them. +#' +#' @param obj a solved rxode2 object. +#' @return a named numeric vector, or `NULL` if nothing was injected. +#' @export +rxInjectedPars <- function(obj) { + .cls <- attr(obj, "class") + .env <- attr(.cls, ".rxode2.env") + if (is.null(.env)) return(NULL) + .rxSolveMaterializeParams(obj, .env) # ensure materialized + if (exists(".injectedPars", envir = .env, inherits = FALSE)) { + get(".injectedPars", envir = .env, inherits = FALSE) + } else { + NULL + } +} + +#' Forced parameters carried by a model / ui +#' +#' Get or set a named-numeric vector of *forced* parameters stored in a hidden +#' slot of an rxode2 ui. Forced parameters override the values supplied in +#' `params`/data and the model initial estimates on **every** solve of that ui +#' (they are written into every subject/simulation column at solve setup, the +#' same injection point used by par-loader hooks). Because the values live on +#' the ui they travel with it -- through model piping and into any `nlmixr2` fit +#' built from it -- so the model stays self-contained and portable (e.g. a fit +#' that carries trained neural-network weights re-solves, predicts and simulates +#' with those weights with no external state). +#' +#' The value is stored directly in the ui environment (not in the printed `meta` +#' block, so it stays hidden) and registered as a `sticky` model item so it +#' survives model piping. Names must be model parameters; names that are not +#' model parameters are ignored at solve time. Set to `NULL` (or an empty +#' vector) to clear. +#' +#' @param ui an rxode2 ui / model function. +#' @param value named numeric vector of forced parameter values, or `NULL`. +#' @return the getter returns the named numeric vector (or `NULL`); the setter +#' returns the (modified) ui. +#' @export +#' @author Matthew L. Fidler +rxForcedPars <- function(ui) { + .ui <- rxUiDecompress(ui) + if (!inherits(.ui, "rxUi")) return(NULL) + if (!exists("forcedPars", envir = .ui, inherits = FALSE)) return(NULL) + get("forcedPars", envir = .ui, inherits = FALSE) +} + +#' @rdname rxForcedPars +#' @export +`rxForcedPars<-` <- function(ui, value) { + .ui <- rxUiDecompress(ui) + if (!inherits(.ui, "rxUi")) { + stop("'ui' must be an rxode2 ui/model to set forcedPars", call. = FALSE) + } + .sticky <- if (exists("sticky", envir = .ui, inherits = FALSE)) { + get("sticky", envir = .ui, inherits = FALSE) + } else character(0) + if (is.null(value) || length(value) == 0L) { + if (exists("forcedPars", envir = .ui, inherits = FALSE)) { + rm("forcedPars", envir = .ui) + } + assign("sticky", setdiff(.sticky, "forcedPars"), envir = .ui) + return(invisible(.ui)) + } + if (is.null(names(value)) || any(names(value) == "")) { + stop("forcedPars must be a fully-named numeric vector", call. = FALSE) + } + # reject non-numeric input up front; as.numeric() would otherwise coerce + # characters/factors to NA and silently force parameters to NA_real_. + if (!is.numeric(value)) { + stop("forcedPars must be numeric (got '", class(value)[1], "')", call. = FALSE) + } + ## store on the ui env (hidden -- not the printed `meta` block) and mark sticky + ## so it survives model piping. + assign("forcedPars", stats::setNames(as.numeric(value), names(value)), + envir = .ui) + assign("sticky", unique(c(.sticky, "forcedPars")), envir = .ui) + invisible(.ui) +} + +## low-level: set/clear the rxode2 forced-parameter buffer honored at solve setup +.rxSetForcedParsC <- function(idx, val) { + invisible(.Call(`_rxode2_rxSetForcedPars`, as.integer(idx), as.double(val))) +} +.rxClearForcedParsC <- function() { + invisible(.Call(`_rxode2_rxClearForcedPars`)) +} + +## resolve a ui's forcedPars names to 0-based solve-param indices and load the C +## buffer for the next solve; returns TRUE if anything was set (so the caller +## clears afterward). No-op (and clears) when the model carries no forcedPars. +## `forcedSrc` supplies the values (the ui); `solveModel` supplies the definitive +## solve-parameter order the gpars layout uses (default: forcedSrc itself). +.rxApplyForcedPars <- function(forcedSrc, solveModel = forcedSrc) { + .fp <- tryCatch(rxForcedPars(forcedSrc), error = function(e) NULL) + if (is.null(.fp) || length(.fp) == 0L) { + .rxClearForcedParsC() + return(FALSE) + } + .pars <- rxModelVars(solveModel)$params + .idx <- match(names(.fp), .pars) - 1L + .keep <- !is.na(.idx) + if (!any(.keep)) { + .rxClearForcedParsC() + return(FALSE) + } + .rxSetForcedParsC(.idx[.keep], unname(.fp)[.keep]) + TRUE +} + +#' Parameter-loader flag on a model +#' +#' A model that needs an externally-owned parameter injector -- a registered +#' par-loader, e.g. a package's neural-network weight loader -- flags that +#' injector's `":"` name here. At solve setup only the loader +#' with that name runs, so an injector never overwrites an unrelated model's +#' `par_ptr` merely because it happens to be registered. The flag is stored on the +#' ui (a sticky item, so it survives model piping) like [rxForcedPars()]; set it to +#' `NULL` to clear. +#' +#' @param ui an rxode2 ui / model function. +#' @param value a single loader name (`":"`), or `NULL`. +#' @return the getter returns the name (or `NULL`); the setter returns the ui. +#' @export +#' @author Matthew L. Fidler +rxParLoader <- function(ui) { + .ui <- rxUiDecompress(ui) + if (!inherits(.ui, "rxUi")) return(NULL) + if (!exists("parLoader", envir = .ui, inherits = FALSE)) return(NULL) + get("parLoader", envir = .ui, inherits = FALSE) +} + +#' @rdname rxParLoader +#' @export +`rxParLoader<-` <- function(ui, value) { + .ui <- rxUiDecompress(ui) + if (!inherits(.ui, "rxUi")) { + stop("'ui' must be an rxode2 ui/model to set rxParLoader", call. = FALSE) + } + .sticky <- if (exists("sticky", envir = .ui, inherits = FALSE)) { + get("sticky", envir = .ui, inherits = FALSE) + } else character(0) + if (is.null(value) || length(value) == 0L || !nzchar(value[1L])) { + if (exists("parLoader", envir = .ui, inherits = FALSE)) rm("parLoader", envir = .ui) + assign("sticky", setdiff(.sticky, "parLoader"), envir = .ui) + return(invisible(.ui)) + } + assign("parLoader", as.character(value[1L]), envir = .ui) + assign("sticky", unique(c(.sticky, "parLoader")), envir = .ui) + invisible(.ui) +} + +## low-level: set/clear the active par-loader name honored at the next solve setup. +.rxSetActiveParLoaderC <- function(name) { + invisible(.Call(`_rxode2_rxSetActiveParLoader`, as.character(name))) +} +.rxClearActiveParLoaderC <- function() { + invisible(.Call(`_rxode2_rxClearActiveParLoader`)) +} + +## bridge: flag the active par-loader from a ui's parLoader item before a solve. +## Returns TRUE if a flag was set (so the caller clears afterward). +.rxApplyParLoader <- function(ui) { + .pl <- tryCatch(rxParLoader(ui), error = function(e) NULL) + if (is.null(.pl) || length(.pl) == 0L || !nzchar(.pl[1L])) return(FALSE) + .rxSetActiveParLoaderC(.pl[1L]) + TRUE +} + +## ---- ui-solve preparation hooks ------------------------------------------- +## Registry of functions called with the ui at the start of a ui solve (before +## parameter loaders run). A package uses this to rehydrate transient C-side +## state from serializable ui slots after a ui has been saved to disk and +## reloaded in a fresh session (e.g. re-registering neural-network shapes so the +## weight values carried in `rxForcedPars()` land in a network that can stride +## them). Each hook must be cheap and a no-op for models it does not own. +.rxUiPrepHooks <- new.env(parent = emptyenv()) + +#' Register or remove a solve-time ui preparation hook +#' +#' Package developers register a function called with the ui at the start of a +#' ui solve, before parameter loaders run. Use it to rehydrate transient C-side +#' state from serializable ui slots (so a ui saved to disk and reloaded in a +#' fresh session solves correctly). The hook must be cheap and a no-op for +#' models it does not own. +#' +#' @param name Unique hook name; re-registering the same name replaces it. +#' @param fn Function of one argument (the rxUi being solved). Return value is +#' ignored; errors are downgraded to a warning so a buggy hook cannot break +#' unrelated solves. +#' @return Invisibly, the hook name (register) or `NULL` (remove). +#' @export +rxRegisterUiPrep <- function(name, fn) { + if (!is.character(name) || length(name) != 1L) { + stop("'name' must be a single string", call. = FALSE) + } + if (!is.function(fn)) { + stop("'fn' must be a function of one argument (the ui)", call. = FALSE) + } + assign(name, fn, envir = .rxUiPrepHooks) + invisible(name) +} + +#' @rdname rxRegisterUiPrep +#' @export +rxRemoveUiPrep <- function(name) { + if (exists(name, envir = .rxUiPrepHooks, inherits = FALSE)) { + rm(list = name, envir = .rxUiPrepHooks) + } + invisible(NULL) +} + +.rxRunUiPrepHooks <- function(object) { + .nm <- ls(envir = .rxUiPrepHooks, all.names = TRUE) + if (length(.nm) == 0L) return(invisible()) + for (.n in .nm) { + .fn <- get(.n, envir = .rxUiPrepHooks) + tryCatch(.fn(object), + error = function(e) { + warning("rxode2 ui-prep hook '", .n, "' failed: ", + conditionMessage(e), call. = FALSE) + }) + } + invisible() +} + .rxSolveGetInit <- function(.env, arg) { .ini <- get(".init.dat", envir = .env, inherits = FALSE) if (arg %in% c("inits", "init")) { diff --git a/inst/include/rxode2.h b/inst/include/rxode2.h index 8062b6e54..4a376aab3 100644 --- a/inst/include/rxode2.h +++ b/inst/include/rxode2.h @@ -92,6 +92,19 @@ typedef rx_solve *(*t_get_solve)(void); typedef void *(*t_assignFuns)(void); +// External parameter-block loader hook signature. A package registers a +// callback (via the rxode2 function-pointer table) that rxode2 invokes once per +// solve, after gpars is filled and before integration, so it can overwrite +// reserved par_ptr slots with externally-owned values (e.g. neural-network +// weights). gpars is laid out `npars` per column with `ncols` columns; write a +// population-constant block to every column. Defined outside the +// __RXODE2PTR_H__ guard so downstream packages see the typedef. +typedef void (*t_rxParLoader)(rx_solve *rx, double *gpars, int npars, int ncols); + +// dydt forcing hook: add forcing to state derivatives from generated model dydt. +// Defined outside the __RXODE2PTR_H__ guard so downstream packages see the typedef. +typedef void (*t_rxDydtForce)(int *neq, double t, double *y, double *dydt); + // Adjoint objective-gradient backward sweep (src/adjoint.cpp). Its address is // exported to downstream packages through the rxode2 function-pointer table // (see _rxode2_rxode2Ptr in src/init.c and rxode2ptr.h); this direct declaration @@ -110,6 +123,15 @@ void rxode2AdjointTrajSweep(double *tg, double *J, double *dP, int ns, int np, rx_solve *getRxSolve_(void); void rxSetSolveAtolRtol(double atol, double rtol); void rxGetSolveAtolRtol(double *atol, double *rtol); + +// External parameter-block loader hook register/remove (t_rxParLoader typedef +// is above, outside the guard). Direct declarations for building rxode2 +// itself; downstream packages get these as function pointers via the table. +void rxRegisterParLoader(t_rxParLoader cb); +void rxRegisterParLoaderNamed(const char* name, t_rxParLoader cb); +void rxRemoveParLoader(t_rxParLoader cb); +void rxRegisterDydtForce(t_rxDydtForce cb); +void rxRemoveDydtForce(t_rxDydtForce cb); int getIndCmt(rx_solving_options* op, rx_solving_options_ind* ind, int kk); #endif rx_solve *getRxSolve2_(void); diff --git a/inst/include/rxode2_model_shared.c b/inst/include/rxode2_model_shared.c index fc7dd8050..efb68539a 100644 --- a/inst/include/rxode2_model_shared.c +++ b/inst/include/rxode2_model_shared.c @@ -20,6 +20,7 @@ _update_par_ptr_p _update_par_ptr=NULL; _getParCov_p _getParCov=NULL; _setThreadInd_t _setThreadInd=NULL; _rxPushDose_t _rxPushDose=NULL; +_rxDydtForce_t _rxDydtForceCb=NULL; linCmtA_p linCmtA; linCmtB_p linCmtB; _rx_asgn _rxode2_rxAssignPtr =NULL; @@ -603,6 +604,7 @@ void _assignFuns0(void) { _getParCov = (_getParCov_p) R_GetCCallable("rxode2","_getParCov"); _setThreadInd = (_setThreadInd_t) R_GetCCallable("rxode2","_setThreadInd"); _rxPushDose = (_rxPushDose_t) R_GetCCallable("rxode2","_rxPushDose"); + _rxDydtForceCb= (_rxDydtForce_t) R_GetCCallable("rxode2","rxCallDydtForce"); // dynamic start linCmtA=(linCmtA_p)R_GetCCallable("rxode2", "linCmtA"); linCmtB=(linCmtB_p)R_GetCCallable("rxode2", "linCmtB"); diff --git a/inst/include/rxode2_model_shared.h b/inst/include/rxode2_model_shared.h index ab761b9f9..0a920a799 100644 --- a/inst/include/rxode2_model_shared.h +++ b/inst/include/rxode2_model_shared.h @@ -382,6 +382,11 @@ typedef double (*linCmtB_p) (rx_solve *rx, int id, typedef void (*_update_par_ptr_p)(double t, unsigned int id, rx_solve *rx, int idx); +/* dydt forcing hook: generated model calls this at the end of its RHS (dydt) so a + plugin can add forcing to state derivatives (e.g. b_j for NN-weight variational + states). neq[0]=nstate, neq[1]=cSub; y=states, dydt=DADT array to add into. */ +typedef void (*_rxDydtForce_t)(int *neq, double t, double *y, double *dydt); + typedef double (*_getParCov_p)(unsigned int id, rx_solve *rx, int parNo, int idx); typedef rx_solve *(*_getRxSolve_t)(void); diff --git a/inst/include/rxode2ptr.h b/inst/include/rxode2ptr.h index b13021a95..1f17599b0 100644 --- a/inst/include/rxode2ptr.h +++ b/inst/include/rxode2ptr.h @@ -274,6 +274,28 @@ extern "C" { typedef void (*setRxThreadId_t)(int id); extern setRxThreadId_t setRxThreadId; + // Register / remove an external parameter-block loader (t_rxParLoader from + // rxode2.h); rxode2 calls registered loaders once per solve so a package can + // fill reserved par_ptr slots with externally-owned values. + typedef void (*rxRegisterParLoader_t)(t_rxParLoader cb); + extern rxRegisterParLoader_t rxRegisterParLoader; + typedef void (*rxRemoveParLoader_t)(t_rxParLoader cb); + extern rxRemoveParLoader_t rxRemoveParLoader; + // Register a NAMED loader (":") that runs only for a model + // flagging that name (rxParLoader()), so an injector never touches an unrelated + // model's par_ptr. + typedef void (*rxRegisterParLoaderNamed_t)(const char* name, t_rxParLoader cb); + extern rxRegisterParLoaderNamed_t rxRegisterParLoaderNamed; + + // Register / remove a dydt forcing hook (t_rxDydtForce from rxode2.h); the + // generated model calls all registered hooks at the end of its RHS so a + // package can add forcing to state derivatives (e.g. NN-weight variational + // states). + typedef void (*rxRegisterDydtForce_t)(t_rxDydtForce cb); + extern rxRegisterDydtForce_t rxRegisterDydtForce; + typedef void (*rxRemoveDydtForce_t)(t_rxDydtForce cb); + extern rxRemoveDydtForce_t rxRemoveDydtForce; + static inline SEXP iniRxodePtrs0(SEXP p) { if (_rxode2_rxRmvnSEXP_ == NULL) { _rxode2_rxRmvnSEXP_ = (_rxode2_rxRmvnSEXP_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 0)); @@ -351,7 +373,7 @@ extern "C" { seedEng = (seedEng_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 72)); rxNormEng = (rxNormEng_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 73)); rxUnifEng = (rxUnifEng_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 81)); - getIndCmt = (getIndCmt_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 82)); + getIndCmt = (getIndCmt_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 87)); setIndSolvePtr = (setIndSolvePtr_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 74)); getIndSolveSave = (getIndSolveSave_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 75)); setIndSolveSave = (setIndSolveSave_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 76)); @@ -359,6 +381,11 @@ extern "C" { setIndSolveLast = (setIndSolveLast_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 78)); getIndSolveLast2 = (getIndSolveLast2_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 79)); setIndSolveLast2 = (setIndSolveLast2_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 80)); + rxRegisterParLoader = (rxRegisterParLoader_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 82)); + rxRemoveParLoader = (rxRemoveParLoader_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 83)); + rxRegisterDydtForce = (rxRegisterDydtForce_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 84)); + rxRemoveDydtForce = (rxRemoveDydtForce_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 85)); + rxRegisterParLoaderNamed = (rxRegisterParLoaderNamed_t) R_ExternalPtrAddrFn(VECTOR_ELT(p, 86)); } return R_NilValue; } @@ -447,6 +474,11 @@ extern "C" { setIndSolveLast_t setIndSolveLast = NULL; \ getIndSolveLast2_t getIndSolveLast2 = NULL; \ setIndSolveLast2_t setIndSolveLast2 = NULL; \ + rxRegisterParLoader_t rxRegisterParLoader = NULL; \ + rxRemoveParLoader_t rxRemoveParLoader = NULL; \ + rxRegisterDydtForce_t rxRegisterDydtForce = NULL; \ + rxRemoveDydtForce_t rxRemoveDydtForce = NULL; \ + rxRegisterParLoaderNamed_t rxRegisterParLoaderNamed = NULL; \ SEXP iniRxodePtrs(SEXP ptr) { \ return iniRxodePtrs0(ptr); \ } \ diff --git a/man-roxygen/rmdhunks/covariates.Rmd b/man-roxygen/rmdhunks/covariates.Rmd index 0c375abb5..a1da4df57 100644 --- a/man-roxygen/rmdhunks/covariates.Rmd +++ b/man-roxygen/rmdhunks/covariates.Rmd @@ -164,3 +164,56 @@ plot(r1,C2, ylab="Central Concentration", xlab="Time") ```{r time-varying-nocb-effect} plot(r1,eff, ylab="Effect", xlab="Time") ``` + +# Forced parameters + +Sometimes a parameter's value is owned *outside* the usual `params`/data/`ini` +channel and must be applied on **every** solve of a model -- for example +weights fit by an external optimizer that should travel with the model so it +stays self-contained. `rxForcedPars()` stores such values in a hidden slot on +the model (a named numeric vector). On each solve they are injected into every +subject/simulation column at solve setup, overriding whatever `params` or the +data supplied: + +```{r} +mod <- function() { + ini({ tcl <- -2 }) + model({ + cl <- exp(tcl) * WT + d/dt(depot) <- -cl * depot + cp <- depot + }) +} +ui <- rxode2(mod) + +ev <- et(amt = 100) |> et(seq(1, 8, by = 1)) +ev$WT <- 1 + +## baseline: WT = 1 from the data +s0 <- rxSolve(ui, ev) + +## force WT = 2 for every solve of this model +rxForcedPars(ui) <- c(WT = 2) +s1 <- rxSolve(ui, ev) + +head(data.frame(time = s0$time, cp_WT1 = s0$cp, cp_forcedWT2 = s1$cp)) +``` + +Because the value lives on the model, it is carried through model piping and +into any `nlmixr2` fit built from the model, so `predict()`/`simulate()` on that +fit reuse it with no external state. The forced values are also reported by +`rxInjectedPars()` on the solved object: + +```{r} +rxInjectedPars(s1) +``` + +Clear forcing by assigning `NULL`: + +```{r} +rxForcedPars(ui) <- NULL +``` + +Names that are not model parameters are ignored, and forcing applies only to +population-constant parameters/covariates (a single value per parameter), not +to per-record time-varying covariates. diff --git a/man/rxForcedPars.Rd b/man/rxForcedPars.Rd new file mode 100644 index 000000000..e6800392f --- /dev/null +++ b/man/rxForcedPars.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/rxsolve.R +\name{rxForcedPars} +\alias{rxForcedPars} +\alias{rxForcedPars<-} +\title{Forced parameters carried by a model / ui} +\usage{ +rxForcedPars(ui) + +rxForcedPars(ui) <- value +} +\arguments{ +\item{ui}{an rxode2 ui / model function.} + +\item{value}{named numeric vector of forced parameter values, or \code{NULL}.} +} +\value{ +the getter returns the named numeric vector (or \code{NULL}); the setter +returns the (modified) ui. +} +\description{ +Get or set a named-numeric vector of \emph{forced} parameters stored in a hidden +slot of an rxode2 ui. Forced parameters override the values supplied in +\code{params}/data and the model initial estimates on \strong{every} solve of that ui +(they are written into every subject/simulation column at solve setup, the +same injection point used by par-loader hooks). Because the values live on +the ui they travel with it -- through model piping and into any \code{nlmixr2} fit +built from it -- so the model stays self-contained and portable (e.g. a fit +that carries trained neural-network weights re-solves, predicts and simulates +with those weights with no external state). +} +\details{ +The value is stored directly in the ui environment (not in the printed \code{meta} +block, so it stays hidden) and registered as a \code{sticky} model item so it +survives model piping. Names must be model parameters; names that are not +model parameters are ignored at solve time. Set to \code{NULL} (or an empty +vector) to clear. +} +\author{ +Matthew L. Fidler +} diff --git a/man/rxInjectedPars.Rd b/man/rxInjectedPars.Rd new file mode 100644 index 000000000..8a5a84b00 --- /dev/null +++ b/man/rxInjectedPars.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/rxsolve.R +\name{rxInjectedPars} +\alias{rxInjectedPars} +\title{Parameters injected into a solve by a par-loader hook} +\usage{ +rxInjectedPars(obj) +} +\arguments{ +\item{obj}{a solved rxode2 object.} +} +\value{ +a named numeric vector, or \code{NULL} if nothing was injected. +} +\description{ +Returns the parameters (e.g. trained neural-network weights) that a +registered par-loader wrote into the solve parameter vector, saved on the +solved object so re-solving from it restores them. +} diff --git a/man/rxParLoader.Rd b/man/rxParLoader.Rd new file mode 100644 index 000000000..b774821da --- /dev/null +++ b/man/rxParLoader.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/rxsolve.R +\name{rxParLoader} +\alias{rxParLoader} +\alias{rxParLoader<-} +\title{Parameter-loader flag on a model} +\usage{ +rxParLoader(ui) + +rxParLoader(ui) <- value +} +\arguments{ +\item{ui}{an rxode2 ui / model function.} + +\item{value}{a single loader name (\code{":"}), or \code{NULL}.} +} +\value{ +the getter returns the name (or \code{NULL}); the setter returns the ui. +} +\description{ +A model that needs an externally-owned parameter injector -- a registered +par-loader, e.g. a package's neural-network weight loader -- flags that +injector's \code{":"} name here. At solve setup only the loader +with that name runs, so an injector never overwrites an unrelated model's +\code{par_ptr} merely because it happens to be registered. The flag is stored on the +ui (a sticky item, so it survives model piping) like \code{\link[=rxForcedPars]{rxForcedPars()}}; set it to +\code{NULL} to clear. +} +\author{ +Matthew L. Fidler +} diff --git a/man/rxRegisterUiPrep.Rd b/man/rxRegisterUiPrep.Rd new file mode 100644 index 000000000..e9908ff0a --- /dev/null +++ b/man/rxRegisterUiPrep.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/rxsolve.R +\name{rxRegisterUiPrep} +\alias{rxRegisterUiPrep} +\alias{rxRemoveUiPrep} +\title{Register or remove a solve-time ui preparation hook} +\usage{ +rxRegisterUiPrep(name, fn) + +rxRemoveUiPrep(name) +} +\arguments{ +\item{name}{Unique hook name; re-registering the same name replaces it.} + +\item{fn}{Function of one argument (the rxUi being solved). Return value is +ignored; errors are downgraded to a warning so a buggy hook cannot break +unrelated solves.} +} +\value{ +Invisibly, the hook name (register) or \code{NULL} (remove). +} +\description{ +Package developers register a function called with the ui at the start of a +ui solve, before parameter loaders run. Use it to rehydrate transient C-side +state from serializable ui slots (so a ui saved to disk and reloaded in a +fresh session solves correctly). The hook must be cheap and a no-op for +models it does not own. +} diff --git a/src/codegen.c b/src/codegen.c index 6abda0416..bb297bda5 100644 --- a/src/codegen.c +++ b/src/codegen.c @@ -755,6 +755,7 @@ void codegen(char *model, int show_ode, const char *prefix, const char *libname, } } if (show_ode == ode_dydt){ + sAppendN(&sbOut, " if (_rxDydtForceCb != NULL) _rxDydtForceCb(_neq, __t, __zzStateVar__, __DDtStateVar__);\n", 90); sAppendN(&sbOut, " (&_solveData->subjects[_cSub])->dadt_counter[0]++;\n}\n\n", 56); } else if (show_ode == ode_jac){ //sAppendN(&sbOut, " free(__ld_DDtStateVar__);\n"); diff --git a/src/init.c b/src/init.c index 520506d44..a58d6bfba 100644 --- a/src/init.c +++ b/src/init.c @@ -275,6 +275,14 @@ void rxOptionsIni(void); void _update_par_ptr(double t, unsigned int id, rx_solve *rx, int idx); double _getParCov(unsigned int id, rx_solve *rx, int parNo, int idx); +/* rxRegisterParLoader / rxRemoveParLoader + t_rxParLoader come from rxode2.h; + they are exported to downstream packages via the function-pointer table + (_rxode2_rxode2Ptr below), not R_RegisterCCallable. */ + +/* dydt forcing hook: rxCallDydtForce is called from generated model code and so + is resolved via R_GetCCallable; rxRegisterDydtForce / rxRemoveDydtForce are + exported to plugins via the function-pointer table (like the par loaders). */ +void rxCallDydtForce(int *neq, double t, double *y, double *dydt); int par_progress(int c, int n, int d, int cores, clock_t t0, int stop); void ind_solve(rx_solve *rx, unsigned int cid, t_dydt_liblsoda dydt_lls, @@ -522,8 +530,13 @@ SEXP _rxode2_rxode2Ptr(void) { SEXP rxode2setIndSolveLast = PROTECT(R_MakeExternalPtrFn((DL_FUNC)&setIndSolveLast, R_NilValue, R_NilValue)); pro++; SEXP rxode2getIndSolveLast2 = PROTECT(R_MakeExternalPtrFn((DL_FUNC)&getIndSolveLast2, R_NilValue, R_NilValue)); pro++; SEXP rxode2setIndSolveLast2 = PROTECT(R_MakeExternalPtrFn((DL_FUNC)&setIndSolveLast2, R_NilValue, R_NilValue)); pro++; + SEXP rxode2rxRegisterParLoader = PROTECT(R_MakeExternalPtrFn((DL_FUNC)&rxRegisterParLoader, R_NilValue, R_NilValue)); pro++; + SEXP rxode2rxRemoveParLoader = PROTECT(R_MakeExternalPtrFn((DL_FUNC)&rxRemoveParLoader, R_NilValue, R_NilValue)); pro++; + SEXP rxode2rxRegisterDydtForce = PROTECT(R_MakeExternalPtrFn((DL_FUNC)&rxRegisterDydtForce, R_NilValue, R_NilValue)); pro++; + SEXP rxode2rxRemoveDydtForce = PROTECT(R_MakeExternalPtrFn((DL_FUNC)&rxRemoveDydtForce, R_NilValue, R_NilValue)); pro++; + SEXP rxode2rxRegisterParLoaderNamed = PROTECT(R_MakeExternalPtrFn((DL_FUNC)&rxRegisterParLoaderNamed, R_NilValue, R_NilValue)); pro++; -#define nVec 83 +#define nVec 88 SEXP ret = PROTECT(Rf_allocVector(VECSXP, nVec)); pro++; SET_VECTOR_ELT(ret, 0, rxode2rxRmvnSEXP); SET_VECTOR_ELT(ret, 1, rxode2rxParProgress); @@ -607,7 +620,12 @@ SEXP _rxode2_rxode2Ptr(void) { SET_VECTOR_ELT(ret, 79, rxode2getIndSolveLast2); SET_VECTOR_ELT(ret, 80, rxode2setIndSolveLast2); SET_VECTOR_ELT(ret, 81, rxode2rxUnifEng); - SET_VECTOR_ELT(ret, 82, rxode2getIndCmt); + SET_VECTOR_ELT(ret, 82, rxode2rxRegisterParLoader); + SET_VECTOR_ELT(ret, 83, rxode2rxRemoveParLoader); + SET_VECTOR_ELT(ret, 84, rxode2rxRegisterDydtForce); + SET_VECTOR_ELT(ret, 85, rxode2rxRemoveDydtForce); + SET_VECTOR_ELT(ret, 86, rxode2rxRegisterParLoaderNamed); + SET_VECTOR_ELT(ret, 87, rxode2getIndCmt); SEXP retN = PROTECT(Rf_allocVector(STRSXP, nVec)); pro++; @@ -693,7 +711,12 @@ SEXP _rxode2_rxode2Ptr(void) { SET_STRING_ELT(retN, 79, Rf_mkChar("rxode2getIndSolveLast2")); SET_STRING_ELT(retN, 80, Rf_mkChar("rxode2setIndSolveLast2")); SET_STRING_ELT(retN, 81, Rf_mkChar("rxode2rxUnifEng")); - SET_STRING_ELT(retN, 82, Rf_mkChar("rxode2getIndCmt")); + SET_STRING_ELT(retN, 82, Rf_mkChar("rxode2rxRegisterParLoader")); + SET_STRING_ELT(retN, 83, Rf_mkChar("rxode2rxRemoveParLoader")); + SET_STRING_ELT(retN, 84, Rf_mkChar("rxode2rxRegisterDydtForce")); + SET_STRING_ELT(retN, 85, Rf_mkChar("rxode2rxRemoveDydtForce")); + SET_STRING_ELT(retN, 86, Rf_mkChar("rxode2rxRegisterParLoaderNamed")); + SET_STRING_ELT(retN, 87, Rf_mkChar("rxode2getIndCmt")); #undef nVec @@ -741,9 +764,24 @@ SEXP _rxode2_rxMemoryComponents_(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SE SEXP _rxode2_rxRamBytes_(void); SEXP _rxode2_rxSolveSetCurObj_(SEXP); +extern SEXP _rxode2_rxRegisterTestParLoaders(SEXP); +extern SEXP _rxode2_rxRemoveTestParLoaders(void); +extern SEXP _rxode2_rxGetInjectedPars(void); +extern SEXP _rxode2_rxSetForcedPars(SEXP, SEXP); +extern SEXP _rxode2_rxClearForcedPars(void); +extern SEXP _rxode2_rxSetActiveParLoader(SEXP); +extern SEXP _rxode2_rxClearActiveParLoader(void); + void R_init_rxode2(DllInfo *info){ allocExtraDosingC(); R_CallMethodDef callMethods[] = { + {"_rxode2_rxRegisterTestParLoaders", (DL_FUNC) &_rxode2_rxRegisterTestParLoaders, 1}, + {"_rxode2_rxRemoveTestParLoaders", (DL_FUNC) &_rxode2_rxRemoveTestParLoaders, 0}, + {"_rxode2_rxGetInjectedPars", (DL_FUNC) &_rxode2_rxGetInjectedPars, 0}, + {"_rxode2_rxSetForcedPars", (DL_FUNC) &_rxode2_rxSetForcedPars, 2}, + {"_rxode2_rxClearForcedPars", (DL_FUNC) &_rxode2_rxClearForcedPars, 0}, + {"_rxode2_rxSetActiveParLoader", (DL_FUNC) &_rxode2_rxSetActiveParLoader, 1}, + {"_rxode2_rxClearActiveParLoader", (DL_FUNC) &_rxode2_rxClearActiveParLoader, 0}, {"_rxode2_qsDes", (DL_FUNC) &_rxode2_qsDes, 1}, {"_rxode2_rxGetSerialType_", (DL_FUNC) &_rxode2_rxGetSerialType_, 1}, {"_rxode2_mlogit_f", (DL_FUNC) &_rxode2_mlogit_f, 2}, @@ -991,6 +1029,7 @@ void R_init_rxode2(DllInfo *info){ R_RegisterCCallable("rxode2", "ind_solve", (DL_FUNC) &ind_solve); R_RegisterCCallable("rxode2", "par_solve", (DL_FUNC) &par_solve); R_RegisterCCallable("rxode2", "_update_par_ptr", (DL_FUNC) &_update_par_ptr); + R_RegisterCCallable("rxode2", "rxCallDydtForce", (DL_FUNC) &rxCallDydtForce); R_RegisterCCallable("rxode2", "_getParCov", (DL_FUNC) &_getParCov); R_RegisterCCallable("rxode2","rxRmModelLib", (DL_FUNC) &rxRmModelLib); R_RegisterCCallable("rxode2","rxGetModelLib", (DL_FUNC) &rxGetModelLib); diff --git a/src/rxData.cpp b/src/rxData.cpp index 025911411..e1845f0bb 100644 --- a/src/rxData.cpp +++ b/src/rxData.cpp @@ -4520,6 +4520,220 @@ static inline void rxSolve_parOrder(const RObject &obj, const List &rxControl, static inline void rxSolve_assignGpars(rxSolve_t* rxSolveDat); +// ---- external parameter-block loader hooks -------------------------------- +// Packages (e.g. rxode2nn) register a callback that runs once per solve, after +// gpars is populated from the supplied parameters and before integration, so +// they can overwrite reserved par_ptr slots with externally-owned values (e.g. +// neural-network weights held in a torch module). This runs single-threaded, +// before the parallel per-subject solve, so callbacks writing gpars are safe. +// gpars is laid out npars per column with `ncols` columns (>= nsub*nsim and +// nPopPar); a population-constant block is written to every column. +typedef void (*t_rxParLoader)(rx_solve* rx, double* gpars, int npars, int ncols); +#define RX_MAX_PAR_LOADERS 32 +static t_rxParLoader _rxParLoaders[RX_MAX_PAR_LOADERS] = {NULL}; +// Each loader has a NAME (empty = "unnamed"). A model that needs a specific +// injector flags its name (via rxParLoader() -> _rxActiveParLoader, below); a +// NAMED loader then runs ONLY for a model that flags it, so an injector (e.g. a +// package's neural-network weight loader) cannot clobber an unrelated model's +// par_ptr just because it happens to be registered. Unnamed loaders keep the +// legacy "run on every solve" behavior for backward compatibility. +static std::string _rxParLoaderNames[RX_MAX_PAR_LOADERS]; +static int _rxNParLoaders = 0; +static std::string _rxActiveParLoader; // set per-solve; empty = no named loader + +static void rxRegisterParLoaderImpl(const char* name, t_rxParLoader cb) { + if (cb == NULL) return; + for (int i = 0; i < _rxNParLoaders; ++i) if (_rxParLoaders[i] == cb) return; + if (_rxNParLoaders < RX_MAX_PAR_LOADERS) { + _rxParLoaderNames[_rxNParLoaders] = (name == NULL) ? std::string() : std::string(name); + _rxParLoaders[_rxNParLoaders++] = cb; + } else { + // registration usually happens during a package's .onLoad(); fail loudly so a + // dropped injector is diagnosable instead of silently never running. + Rf_warning("rxode2: parameter-loader registry full (max %d); loader '%s' not registered", + RX_MAX_PAR_LOADERS, (name == NULL) ? "" : name); + } +} + +extern "C" void rxRegisterParLoader(t_rxParLoader cb) { + rxRegisterParLoaderImpl(NULL, cb); // unnamed -> legacy always-run +} + +// Register a loader under a ":" name so it dispatches only to +// models that flag that name. +extern "C" void rxRegisterParLoaderNamed(const char* name, t_rxParLoader cb) { + rxRegisterParLoaderImpl(name, cb); +} + +extern "C" void rxRemoveParLoader(t_rxParLoader cb) { + for (int i = 0; i < _rxNParLoaders; ++i) { + if (_rxParLoaders[i] == cb) { + for (int k = i; k < _rxNParLoaders - 1; ++k) { + _rxParLoaders[k] = _rxParLoaders[k + 1]; + _rxParLoaderNames[k] = _rxParLoaderNames[k + 1]; + } + _rxParLoaders[--_rxNParLoaders] = NULL; + _rxParLoaderNames[_rxNParLoaders].clear(); + return; + } + } +} + +// The active injector flag for the next solve: set from the model about to be +// solved (rxSetActiveParLoader), READ by rxCallParLoaders to dispatch to the +// matching named loader, and cleared afterward by the R on-exit hook +// .rxClearActiveParLoaderC() -- rxCallParLoaders does not clear it itself. +extern "C" SEXP _rxode2_rxSetActiveParLoader(SEXP nameSxp) { + _rxActiveParLoader = (TYPEOF(nameSxp) == STRSXP && Rf_length(nameSxp) >= 1) ? + std::string(CHAR(STRING_ELT(nameSxp, 0))) : std::string(); + return R_NilValue; +} +extern "C" SEXP _rxode2_rxClearActiveParLoader(void) { + _rxActiveParLoader.clear(); + return R_NilValue; +} + +// ---- dydt forcing hooks --------------------------------------------------- +// Packages register a callback the generated model invokes at the end of its RHS +// (dydt) evaluation, so they can ADD forcing to designated state derivatives -- +// e.g. the b_j = dR/dg * dg/dw_j term for NN-weight forward-sensitivity +// (variational) states, on top of the J*s_j part rxode2's sensitivity codegen +// already produces. Called inside the per-subject solve (possibly parallel), so +// callbacks must be thread-safe (read-only shared state, write only their own +// dydt slots). neq[0]=nstate, neq[1]=cSub. +typedef void (*t_rxDydtForce)(int *neq, double t, double *y, double *dydt); +#define RX_MAX_DYDT_FORCE 32 +static t_rxDydtForce _rxDydtForce[RX_MAX_DYDT_FORCE] = {NULL}; +static int _rxNDydtForce = 0; + +extern "C" void rxRegisterDydtForce(t_rxDydtForce cb) { + if (cb == NULL) return; + for (int i = 0; i < _rxNDydtForce; ++i) if (_rxDydtForce[i] == cb) return; + if (_rxNDydtForce < RX_MAX_DYDT_FORCE) _rxDydtForce[_rxNDydtForce++] = cb; +} + +extern "C" void rxRemoveDydtForce(t_rxDydtForce cb) { + for (int i = 0; i < _rxNDydtForce; ++i) { + if (_rxDydtForce[i] == cb) { + for (int k = i; k < _rxNDydtForce - 1; ++k) _rxDydtForce[k] = _rxDydtForce[k + 1]; + _rxDydtForce[--_rxNDydtForce] = NULL; + return; + } + } +} + +// Invoked by generated model dydt (resolved via R_GetCCallable at model load). +// Fast no-op when nothing is registered -- one branch per RHS evaluation. +extern "C" void rxCallDydtForce(int *neq, double t, double *y, double *dydt) { + for (int i = 0; i < _rxNDydtForce; ++i) _rxDydtForce[i](neq, t, y, dydt); +} + +// Parameters injected by the loaders on the most recent solve, captured by +// diffing the population parameter block (column 0) before/after the loaders. +// Lets the solved object save/restore externally-injected values (e.g. trained +// neural-network weights) so re-solving from the object reproduces them even in +// a session where the injecting package's buffer is gone. +static std::vector _rxInjIdx; +static std::vector _rxInjVal; + +// Forced parameters for the current solve: (0-based param index, value) pairs set +// from R (a model/ui's `forcedPars` slot resolved to solve-param indices) and +// injected into EVERY gpars column before the registered par-loaders run. This +// is the first-class, plugin-free forcing mechanism: a ui carries its forced +// values (e.g. trained NN weights) and every solve of that ui reproduces them, so +// a fit is self-contained and portable. Population-constant (same value in every +// column), single-threaded (set before the parallel per-subject solve). +static std::vector _rxForcedIdx; +static std::vector _rxForcedVal; + +extern "C" SEXP _rxode2_rxSetForcedPars(SEXP idx, SEXP val) { + int n = Rf_length(idx); + if (Rf_length(val) != n) Rf_error("forcedPars idx/val length mismatch"); + SEXP idxI = PROTECT(Rf_coerceVector(idx, INTSXP)); + SEXP valR = PROTECT(Rf_coerceVector(val, REALSXP)); + _rxForcedIdx.assign(INTEGER(idxI), INTEGER(idxI) + n); + _rxForcedVal.assign(REAL(valR), REAL(valR) + n); + UNPROTECT(2); + return R_NilValue; +} + +extern "C" SEXP _rxode2_rxClearForcedPars(void) { + _rxForcedIdx.clear(); + _rxForcedVal.clear(); + return R_NilValue; +} + +static inline void rxCallParLoaders(rx_solve* rx, int npars, int ncols) { + _rxInjIdx.clear(); + _rxInjVal.clear(); + if (_rxNParLoaders == 0 && _rxForcedIdx.empty()) return; + static std::vector pre; + pre.assign(&_globals.gpars[0], &_globals.gpars[0] + npars); // subject 0 block + // ui-driven forced parameters: write to every subject/sim column first, so a + // registered loader (if any) can still override them. + for (size_t f = 0; f < _rxForcedIdx.size(); ++f) { + int k = _rxForcedIdx[f]; + if (k < 0 || k >= npars) continue; + double v = _rxForcedVal[f]; + for (int c = 0; c < ncols; ++c) _globals.gpars[(size_t)c * npars + k] = v; + } + for (int i = 0; i < _rxNParLoaders; ++i) { + // an UNNAMED loader runs always (legacy); a NAMED loader runs only when the + // model being solved flags its name (_rxActiveParLoader), so a package's + // injector never touches an unrelated model's par_ptr. + if (!_rxParLoaderNames[i].empty() && _rxParLoaderNames[i] != _rxActiveParLoader) continue; + _rxParLoaders[i](rx, &_globals.gpars[0], npars, ncols); + } + for (int k = 0; k < npars; ++k) { + if (_globals.gpars[k] != pre[k]) { + _rxInjIdx.push_back(k); + _rxInjVal.push_back(_globals.gpars[k]); + } + } +} + +// Injected params from the last solve as list(index0 = <0-based par indices>, +// value = ); the R layer maps indices to parameter names. +extern "C" SEXP _rxode2_rxGetInjectedPars(void) { + int n = (int) _rxInjIdx.size(); + SEXP idx = PROTECT(Rf_allocVector(INTSXP, n)); + SEXP val = PROTECT(Rf_allocVector(REALSXP, n)); + for (int i = 0; i < n; ++i) { + INTEGER(idx)[i] = _rxInjIdx[i]; + REAL(val)[i] = _rxInjVal[i]; + } + SEXP ret = PROTECT(Rf_allocVector(VECSXP, 2)); + SET_VECTOR_ELT(ret, 0, idx); + SET_VECTOR_ELT(ret, 1, val); + UNPROTECT(3); + return ret; +} + +// Test-only par loaders (tests/testthat/test-par-loader.R): confirm that +// multiple registered loaders are applied in series -- A writes a sentinel to +// parameter 0, B to parameter 1. +extern "C" void rxTestParLoaderA(rx_solve* rx, double* gpars, int npars, int ncols) { + (void) rx; + if (npars < 1) return; + for (int c = 0; c < ncols; ++c) gpars[(size_t) c * npars + 0] = 111.0; +} +extern "C" void rxTestParLoaderB(rx_solve* rx, double* gpars, int npars, int ncols) { + (void) rx; + if (npars < 2) return; + for (int c = 0; c < ncols; ++c) gpars[(size_t) c * npars + 1] = 222.0; +} +extern "C" SEXP _rxode2_rxRegisterTestParLoaders(SEXP nSEXP) { + int n = Rf_asInteger(nSEXP); + rxRegisterParLoader(rxTestParLoaderA); + if (n >= 2) rxRegisterParLoader(rxTestParLoaderB); + return R_NilValue; +} +extern "C" SEXP _rxode2_rxRemoveTestParLoaders(void) { + rxRemoveParLoader(rxTestParLoaderA); + rxRemoveParLoader(rxTestParLoaderB); + return R_NilValue; +} + static inline void rxSolve_resample(const RObject &obj, const List &rxControl, const Nullable &specParams, @@ -4715,6 +4929,11 @@ static inline void rxSolve_normalizeParms(const RObject &obj, const List &rxCont { gparsCovSetup(rx->npars, rxSolveDat->nPopPar, rx->nsub*rx->nsim, ev1, rx); rxSolve_assignGpars(rxSolveDat); + { + int _ncols = (int)(rx->nsub*rx->nsim); + if (rxSolveDat->nPopPar > _ncols) _ncols = rxSolveDat->nPopPar; + rxCallParLoaders(rx, rx->npars, _ncols); + } rxSolve_resample(obj, rxControl, specParams, extraArgs, pars, ev1, inits, rxSolveDat); curSolve=0; diff --git a/tests/testthat/test-forced-pars.R b/tests/testthat/test-forced-pars.R new file mode 100644 index 000000000..109d90930 --- /dev/null +++ b/tests/testthat/test-forced-pars.R @@ -0,0 +1,87 @@ +## Per-ui forced parameters: a named-numeric slot carried on the ui whose values +## override params/data/inits on EVERY solve (injected into every gpars column at +## solve setup, the par-loader injection point). This is the plugin-free forcing +## mechanism that lets a model/fit carry externally-owned values (e.g. trained +## neural-network weights) and stay self-contained. + +## slow-decay model: cl = exp(tcl)*WT with tcl small, so cp stays well away from 0 +## across the observed window (a robust, non-saturating response to forcing WT). +.forcedMod <- function() { + ini({ tcl <- -2 }) + model({ + cl <- exp(tcl) * WT + d/dt(depot) <- -cl * depot + cp <- depot + }) +} + +.forcedEv <- function() { + ev <- et(amt = 100) + ev <- et(ev, seq(1, 8, by = 1)) + ev$WT <- 1 + ev +} + +test_that("rxForcedPars getter/setter round-trips and validates", { + ui <- rxode2(.forcedMod) + expect_null(rxForcedPars(ui)) + rxForcedPars(ui) <- c(WT = 2) + expect_equal(rxForcedPars(ui), c(WT = 2)) + ## unnamed -> error + expect_error(`rxForcedPars<-`(ui, 3), "fully-named") + ## clear + rxForcedPars(ui) <- NULL + expect_null(rxForcedPars(ui)) +}) + +test_that("forcedPars overrides the data covariate on rxSolve; unset -> unchanged", { + ui <- rxode2(.forcedMod) + ev <- .forcedEv() + + s0 <- rxSolve(ui, ev) + + ## force WT = 2 -> cl doubles -> faster decay -> strictly lower cp at every t>0 + rxForcedPars(ui) <- c(WT = 2) + s1 <- rxSolve(ui, ev) + expect_true(all(s1$cp < s0$cp)) + + ## the injected value is captured on the solved object too + expect_equal(unname(rxInjectedPars(s1)["WT"]), 2) + + ## unset -> identical to the un-forced solve + rxForcedPars(ui) <- NULL + s2 <- rxSolve(ui, ev) + expect_equal(s2$cp, s0$cp, tolerance = 1e-8) +}) + +test_that("forcedPars is hidden from the printed model and survives piping", { + ui <- rxode2(.forcedMod) + rxForcedPars(ui) <- c(WT = 2) + + ## hidden: stored on the ui env, not the printed meta block + .txt <- paste(capture.output(print(ui)), collapse = "\n") + expect_false(grepl("forcedPars", .txt)) + + ## survives a model-block change (a "significant" pipe) via the sticky mechanism + ui2 <- ui |> model(cp <- depot * 1) + expect_equal(rxForcedPars(ui2), c(WT = 2)) + ## still hidden after piping + .txt2 <- paste(capture.output(print(ui2)), collapse = "\n") + expect_false(grepl("forcedPars", .txt2)) + + ## and still applied on the piped model's solve + ev <- .forcedEv() + s0 <- rxSolve(rxode2(.forcedMod), ev) + s2 <- rxSolve(ui2, ev) + expect_true(all(s2$cp < s0$cp)) +}) + +test_that("forcedPars names that are not model params are ignored", { + ui <- rxode2(.forcedMod) + ev <- .forcedEv() + s0 <- rxSolve(ui, ev) + ## a bogus name alongside a real one: bogus ignored, WT applied + rxForcedPars(ui) <- c(nonParam = 99, WT = 2) + s1 <- rxSolve(ui, ev) + expect_true(all(s1$cp < s0$cp)) +}) diff --git a/tests/testthat/test-occ.R b/tests/testthat/test-occ.R index 58b5b184b..c48c6aec7 100644 --- a/tests/testthat/test-occ.R +++ b/tests/testthat/test-occ.R @@ -383,6 +383,11 @@ rxTest({ expect_equal(curEval[curEval$parameter == "iov.cl", "curEval"], "exp") + # the post-UDF eta refresh must not fold IOV etas into the id-level eta + # list; that makes `theta + eta + iov` parse as 2 population etas and fail + expect_false("iov.cl" %in% f$eta) + expect_true("iov.cl" %in% f$level) + }) test_that("iov simulation with id()/occ() nesting (issue #323)", { diff --git a/tests/testthat/test-par-loader.R b/tests/testthat/test-par-loader.R new file mode 100644 index 000000000..9a0344309 --- /dev/null +++ b/tests/testthat/test-par-loader.R @@ -0,0 +1,86 @@ +## External parameter-block loader hooks (rxRegisterParLoader): a package can +## register callbacks that rxode2 invokes once per solve, after gpars is filled +## and before integration, to overwrite reserved par_ptr slots. Verifies that +## MULTIPLE registered loaders are applied in series (test loaders A and B write +## sentinels 111 and 222 to parameters 0 and 1). + +rxTest({ + + test_that("multiple par-loaders are applied in series", { + + .m <- rxode2({ + param(a, b) + oa <- a + ob <- b + d/dt(x) <- 0 + }) + ## a, b are the first two parameters -> par_ptr indices 0 and 1 + .pars <- rxModelVars(.m)$params + expect_equal(.pars[1], "a") + expect_equal(.pars[2], "b") + + .ev <- et(amt = 0) |> et(0, 1, by = 1) + + ## baseline: supplied parameter values pass through unchanged + .s0 <- rxSolve(.m, .ev, params = c(a = 1, b = 2), returnType = "data.frame") + expect_equal(.s0$oa[1], 1) + expect_equal(.s0$ob[1], 2) + + ## register two loaders -> both overwrite their slot, in series + .Call("_rxode2_rxRegisterTestParLoaders", 2L, PACKAGE = "rxode2") + on.exit(.Call("_rxode2_rxRemoveTestParLoaders", PACKAGE = "rxode2"), add = TRUE) + + .s <- rxSolve(.m, .ev, params = c(a = 1, b = 2), returnType = "data.frame") + expect_equal(.s$oa[1], 111) # loader A wrote parameter 0 + expect_equal(.s$ob[1], 222) # loader B wrote parameter 1 (second in series) + + ## a single loader overwrites only its slot + .Call("_rxode2_rxRemoveTestParLoaders", PACKAGE = "rxode2") + .Call("_rxode2_rxRegisterTestParLoaders", 1L, PACKAGE = "rxode2") + .s1 <- rxSolve(.m, .ev, params = c(a = 1, b = 2), returnType = "data.frame") + expect_equal(.s1$oa[1], 111) # loader A only + expect_equal(.s1$ob[1], 2) + + ## removing the loaders restores pass-through behavior + .Call("_rxode2_rxRemoveTestParLoaders", PACKAGE = "rxode2") + .s2 <- rxSolve(.m, .ev, params = c(a = 1, b = 2), returnType = "data.frame") + expect_equal(.s2$oa[1], 1) + expect_equal(.s2$ob[1], 2) + }) + + test_that("injected parameters are saved on the object and restored on re-solve", { + + .m <- rxode2({ + param(a, b) + oa <- a + ob <- b + d/dt(x) <- 0 + }) + .ev <- et(amt = 0) |> et(0, 1, by = 1) + + ## solve with two loaders injecting 111 -> a, 222 -> b + .Call("_rxode2_rxRegisterTestParLoaders", 2L, PACKAGE = "rxode2") + .obj <- rxSolve(.m, .ev, params = c(a = 1, b = 2)) + + ## the injected values are saved on the solved object + .inj <- rxInjectedPars(.obj) + expect_equal(.inj[["a"]], 111) + expect_equal(.inj[["b"]], 222) + + ## remove the loaders: a plain solve passes the supplied params through + .Call("_rxode2_rxRemoveTestParLoaders", PACKAGE = "rxode2") + .plain <- rxSolve(.m, .ev, params = c(a = 1, b = 2), returnType = "data.frame") + expect_equal(.plain$oa[1], 1) + expect_equal(.plain$ob[1], 2) + + ## re-solving from the saved object restores the injected values, even though + ## no loader is registered anymore + .re <- rxSolve(.obj, .ev, returnType = "data.frame") + expect_equal(.re$oa[1], 111) + expect_equal(.re$ob[1], 222) + + ## a model with no injection reports nothing + expect_null(rxInjectedPars(.plain)) + }) + +}) diff --git a/vignettes/articles/rxode2-solve-hooks.Rmd b/vignettes/articles/rxode2-solve-hooks.Rmd new file mode 100644 index 000000000..92bfb7f23 --- /dev/null +++ b/vignettes/articles/rxode2-solve-hooks.Rmd @@ -0,0 +1,197 @@ +--- +title: "Solve-time hooks for package developers" +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + eval = FALSE +) +``` + +## Overview + +Most `rxode2` models are fully specified by their parameters, data and initial +conditions. Some downstream tools, though, need to change *what a solve sees* from +outside the model text -- for example to inject a block of externally-owned +parameters (neural-network weights, a lookup table) on every solve, or to add a +forcing term to a state derivative. `rxode2` exposes a small set of hooks for this. +They fall into two groups: + +* **R-level:** `rxForcedPars()` -- override parameter values on a model, carried + with the model itself. +* **C-level:** the *par-loader* and *dydt forcing* hooks, plus the + **function-pointer table** that lets a downstream package call `rxode2`'s C + entry points. These are for package authors writing compiled code. + +This article is a companion to +[Providing Custom C/C++ Functions to rxode2](rxode2-pkg-exported-funs.html): that +one adds new *functions* callable inside a model; this one changes the *parameters +and derivatives* a solve uses. + +## `rxForcedPars()`: forced parameters on a model + +`rxForcedPars()` sets a named vector of parameter values that override `params`/data +and the model initial estimates on **every** solve of a ui. The values are stored +on the ui (hidden from the printed model, registered as a `sticky` item), so they +travel with it -- through model piping and into any `nlmixr2` fit built from it. +That keeps a model self-contained: a fit that carries, say, trained weights +re-solves, predicts and simulates with those weights and no external state. + +```{r} +ui <- function() { + model({ + d/dt(center) <- -(cl/v) * center + cp <- center/v + }) +} +ui <- rxode2(ui) +rxForcedPars(ui) <- c(cl = 1.2, v = 8) # force cl, v on every solve +rxForcedPars(ui) # getter +rxForcedPars(ui) <- NULL # clear +``` + +Names that are not model parameters are ignored at solve time. Forced parameters +are written into every subject/simulation column at solve setup -- the *same +injection point* the par-loader hook (below) uses, so a registered loader can still +override them if needed. + +## The par-loader hook (C) + +For a parameter block that is not a fixed vector but is *computed* -- owned by a +plugin and refreshed each solve -- register a **par-loader**. It runs once per +solve, single-threaded, right after `rxode2` lays out the global parameter matrix +(`gpars`) and before the parallel integration, so the injected values are in place +for the (thread-safe) solve. + +```c +/* signature (inst/include/rxode2.h) */ +typedef void (*t_rxParLoader)(rx_solve *rx, double *gpars, int npars, int ncols); +void rxRegisterParLoader(t_rxParLoader cb); +void rxRemoveParLoader(t_rxParLoader cb); +``` + +`gpars` is the flat `ncols x npars` parameter matrix (column `c`'s `par_ptr` is +`&gpars[c*npars]`). A loader writes its block into every column: + +```c +static void myLoader(rx_solve *rx, double *gpars, int npars, int ncols) { + for (int c = 0; c < ncols; ++c) { + double *par_ptr = gpars + (size_t)c * npars; + for (int k = 0; k < BLOCK_N; ++k) par_ptr[base + k] = myBlock[k]; /* base = block start */ + } +} +``` + +Register it in `.onLoad()` (via a `.Call` wrapper) and `rxRemoveParLoader()` in +`.onUnload()`. The block must occupy real `par_ptr` slots: declare its parameters +with `param(...)` (or as covariates) so `rxode2` reserves contiguous positions, +resolve `base` from the model's solve parameter order once, and cache it. +`rxInjectedPars()` reports what a loader changed on the last solve. + +### Named loaders: dispatch to the right model only + +A loader registered with `rxRegisterParLoader()` runs on **every** solve, so it must +guard itself, and it can still corrupt an unrelated model if its guards are loose. +Prefer a **named** loader instead: + +```c +rxRegisterParLoaderNamed("myPkg:myLoader", myLoader); /* runs only when flagged */ +``` + +A named loader runs *only* while a model flags its name. Flag the injector on the +model that needs it with `rxParLoader()` (a sticky ui item, like `rxForcedPars()`): + +```{r} +rxParLoader(ui) <- "myPkg:myLoader" # only this loader fires for ui's solves +``` + +`rxSolve()` sets that flag active for the model's solve and clears it afterward, so +the loader never touches any other model. For solves that bypass the ui path +(estimation internals, plain `rxode2()` models), a package can set the active flag +directly around its own solve batch and clear it after. This is how `nlmixr2nn` +keeps its neural-network weight loader from leaking into unrelated fits. + +## The dydt forcing hook (C) + +To add a term to a state derivative that the model text cannot express (e.g. a +plugin's contribution computed in C), register a **dydt-force** callback. The +generated model calls it at the end of `dydt`, so the added forcing is integrated +like any other RHS term: + +```c +/* signature (inst/include/rxode2.h) */ +typedef void (*t_rxDydtForce)(int *neq, double t, double *y, double *dydt); +void rxRegisterDydtForce(t_rxDydtForce cb); +void rxRemoveDydtForce(t_rxDydtForce cb); +``` + +```c +static void myForce(int *neq, double t, double *y, double *dydt) { + dydt[stateIdx] += myForcingTerm(t, y); /* on top of the model's own RHS */ +} +``` + +A model with no registered forcing is unaffected (the call is a null-pointer +no-op), so this is safe to leave installed. + +## `rxRegisterUiPrep()`: rehydrating state before a ui solve + +The C hooks above rely on state a package registers *in the running session*. +That state is lost when a model or fit is saved to disk and reloaded in a fresh +session -- the serializable ui survives, but the C-side registry behind it does +not. `rxRegisterUiPrep(name, fn)` closes that gap: `fn(ui)` is called with the +ui at the very start of every ui solve (before parameter loaders run), so a +package can rebuild transient C state from serializable ui slots. + +```r +# in a downstream package's .onLoad(): +rxode2::rxRegisterUiPrep("myPkg:rehydrate", function(ui) { + meta <- ui$myMeta # a slot the package stamped on the ui at fit time + if (is.null(meta)) return(invisible()) # no-op for models it does not own + # ... resolve indices BY NAME from the current parameter layout and re-register + # the C-side state so the reloaded ui solves correctly ... +}) +``` + +Guidelines for a prep hook: + +* **Be cheap and a no-op for unrelated models** -- it runs on *every* ui solve. +* **Resolve positions by name, not by a stored index.** A ui saved by one build + and reloaded by another may have a different parameter column layout; matching + a saved *name* against the current layout is robust, a cached integer is not. +* Errors are downgraded to a warning so a buggy hook cannot break other solves. + +This is how `nlmixr2nn` makes a saved neural-network fit reload: the trained +weights ride on the ui in `rxForcedPars()` and the network *shapes* in a sticky +`nnMeta` slot; the prep hook re-registers the shapes (base resolved by name) so +the persisted weight values land in a network that can stride them. Pair a +`rxRegisterUiPrep()` with `rxRemoveUiPrep(name)` in `.onUnload()`. + +## The function-pointer table + +The hooks above are C callables *inside* `rxode2`. Because CRAN discourages +`R_RegisterCCallable` ABI coupling across packages, `rxode2` shares its C entry +points through a positional **function-pointer table** rather than by name: + +* `rxode2::.rxode2ptrs()` returns the table (an external-pointer list) built by + `_rxode2_rxode2Ptr()` in `src/init.c`. +* A downstream package installs it into its own globals with `iniRxodePtrs()` / + `iniRxodePtrs0()` from the header `inst/include/rxode2ptr.h`, calling + `.Call("_myPkg_iniRxodePtrs", rxode2::.rxode2ptrs())` in `.onLoad()`. + +Each table slot is a fixed index; new entries are appended (so the layout stays +backward compatible). `rxRegisterParLoader` / `rxRegisterDydtForce` are exposed +this way, so a plugin can register its hooks from its own compiled code. See the +`rxode2` source (`src/init.c`, `inst/include/rxode2ptr.h`) for the current table, +and the `nlmixr2nn` package for a worked consumer that uses the par-loader to inject +neural-network weights. + +## See also + +* [Providing Custom C/C++ Functions to rxode2](rxode2-pkg-exported-funs.html) -- + adding model-callable functions. +* `nlmixr2est`: *"Extending nlmixr2est"* -- creating estimation methods, + intercepting estimations, and contributing to the objective, which build on + these solve hooks.