From ae998deaa5e5f019494103c953034beef0baacab Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Fri, 1 Nov 2024 15:16:31 -0500 Subject: [PATCH] Add more stability enhancements; can now compare 3 forms of rxode2 solve --- NEWS.md | 5 ++++ R/RcppExports.R | 4 +++ R/poped.R | 16 ++++++++++++ R/zzz.R | 1 + src/RcppExports.cpp | 10 ++++++++ src/init.c | 3 +++ src/poped.cpp | 47 +++++++++++++++++++++++++++++++++ vignettes/articles/PopED.Rmd | 50 +++++++++++++++++++++++++++++++----- 8 files changed, 129 insertions(+), 7 deletions(-) diff --git a/NEWS.md b/NEWS.md index a4e1af09..0b52177e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,10 @@ # babelmixr2 (development version) +* Check loaded `rxode2` information and compare to what the loaded + model information should be. This allows better checking of which + model is loaded and even more robust stability. It requires + `rxode2` > `3.0.2`. + # babelmixr2 0.1.5 * Fix bug where `PopED` could error with certain `dvid` values diff --git a/R/RcppExports.R b/R/RcppExports.R index 29584c78..929ca33b 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -145,6 +145,10 @@ popedFree <- function() { .Call(`_babelmixr2_popedFree`) } +popedGetLoadedInfo <- function() { + .Call(`_babelmixr2_popedGetLoadedInfo`) +} + popedSetup <- function(e, eglobal, full) { .Call(`_babelmixr2_popedSetup`, e, eglobal, full) } diff --git a/R/poped.R b/R/poped.R index 096380b9..4b2b96f2 100644 --- a/R/poped.R +++ b/R/poped.R @@ -6,6 +6,7 @@ #' @return nothing, called for side effects #' #' @export +#' #' @author Matthew L. Fidler #' @keywords internal .popedFree <- function() { @@ -418,6 +419,11 @@ attr(rxUiGet.popedFfFun, "desc") <- "PopED parameter model (ff_fun)" } else if (.poped$curNumber != popedDb$babelmixr2$modelNumber) { .poped$setup <- 0L } + if (.poped$setup != 0L && + !identical(.poped$loadInfo, + popedGetLoadedInfo())) { + .poped$setup <- 0L + } if (.poped$setup != 1L) { rxode2::rxSolveFree() .popedSetup(popedDb$babelmixr2, .poped, FALSE) @@ -447,6 +453,11 @@ attr(rxUiGet.popedFfFun, "desc") <- "PopED parameter model (ff_fun)" if (.poped$curNumber != popedDb$babelmixr2$modelNumber) { .poped$setup <- 0L } + if (.poped$setup != 0L && + !identical(.poped$loadInfo, + popedGetLoadedInfo())) { + .poped$setup <- 0L + } if (.poped$setup == 2L) { if (!identical(.poped$fullXt, length(xt))) { .poped$setup <- 0L @@ -518,6 +529,11 @@ attr(rxUiGet.popedFfFun, "desc") <- "PopED parameter model (ff_fun)" if (.poped$curNumber != popedDb$babelmixr2$modelNumber) { .poped$setup <- 0L } + if (.poped$setup != 0L && + !identical(.poped$loadInfo, + popedGetLoadedInfo())) { + .poped$setup <- 0L + } if (.poped$setup == 2L) { if (!identical(.poped$fullXt, length(xt))) { .poped$setup <- 0L diff --git a/R/zzz.R b/R/zzz.R index 0b3205de..56e94238 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -30,4 +30,5 @@ rxode2.api <- names(rxode2::.rxode2ptrs()) rxode2::.s3register("rxode2::rxUiDeparse", "pkncaControl") rxode2::.s3register("rxode2::rxUiDeparse", "popedControl") .iniRxode2Ptr() + .poped$loadInfo <- popedGetLoadedInfo() } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 1fff9d19..82cd196d 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -93,6 +93,16 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// popedGetLoadedInfo +Rcpp::IntegerVector popedGetLoadedInfo(); +RcppExport SEXP _babelmixr2_popedGetLoadedInfo() { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + rcpp_result_gen = Rcpp::wrap(popedGetLoadedInfo()); + return rcpp_result_gen; +END_RCPP +} // popedSetup RObject popedSetup(Environment e, Environment eglobal, bool full); RcppExport SEXP _babelmixr2_popedSetup(SEXP eSEXP, SEXP eglobalSEXP, SEXP fullSEXP) { diff --git a/src/init.c b/src/init.c index 4601685e..d5a93f57 100644 --- a/src/init.c +++ b/src/init.c @@ -21,7 +21,10 @@ SEXP _babelmixr2_popedMultipleEndpointResetTimeIndex(void); SEXP _babelmixr2_popedMultipleEndpointIndexDataFrame(SEXP); SEXP _babelmixr2_popedMultipleEndpointParam(SEXP, SEXP, SEXP, SEXP, SEXP); +SEXP _babelmixr2_popedGetLoadedInfo(void); + static const R_CallMethodDef CallEntries[] = { + {"_babelmixr2_popedGetLoadedInfo", (DL_FUNC) &_babelmixr2_popedGetLoadedInfo, 0}, {"_babelmixr2_popedMultipleEndpointParam", (DL_FUNC) &_babelmixr2_popedMultipleEndpointParam, 5}, {"_babelmixr2_popedMultipleEndpointIndexDataFrame", (DL_FUNC) &_babelmixr2_popedMultipleEndpointIndexDataFrame, 1}, diff --git a/src/poped.cpp b/src/poped.cpp index 96c6c540..7fcfc170 100644 --- a/src/poped.cpp +++ b/src/poped.cpp @@ -349,6 +349,50 @@ RObject popedFree() { return R_NilValue; } +//[[Rcpp::export]] +Rcpp::IntegerVector popedGetLoadedInfo() { + rx = getRxSolve_(); + if (rx == NULL) { + return Rcpp::IntegerVector::create(_["nsub"]=0, + _["nall"]=0, + _["nobs"]=0, + _["nobs2"]=0, + _["neq"]=0, + _["nlhs"]=0, + _["stiff"]=0, + _["npars"]=0); + } + Rcpp::IntegerVector ret(8); + Rcpp::CharacterVector retN(8); + retN[0] = "nsub"; + ret[0] = getRxNsub(rx); + + retN[1] = "nall"; + ret[1] = getRxNall(rx); + + retN[2] = "nobs"; + ret[2] = getRxNobs(rx); + + retN[3] = "nobs2"; + ret[3] = getRxNobs2(rx); + + rx_solving_options *op = getSolvingOptions(rx); + + retN[4] = "neq"; + ret[4] = getOpNeq(op); + + retN[5] = "nlhs"; + ret[5] = getOpNlhs(op); + + retN[6] = "stiff"; + ret[6] = getOpStiff(op); + + retN[7] = "npars"; + ret[7] = getRxNpars(rx); + ret.attr("names") = retN; + return ret; +} + //[[Rcpp::export]] RObject popedSetup(Environment e, Environment eglobal, bool full) { @@ -377,6 +421,7 @@ RObject popedSetup(Environment e, Environment eglobal, bool full) { List mvp = rxode2::rxModelVars_(model); CharacterVector trans = mvp["trans"]; _popedEglobal["curTrans"] = trans; + rxUpdateFn(as(trans)); // initial value of parameters @@ -406,9 +451,11 @@ RObject popedSetup(Environment e, Environment eglobal, bool full) { R_NilValue, // inits 1);//const int setupOnly = 0 rx = getRxSolve_(); + _popedEglobal["loadInfo"] = popedGetLoadedInfo(); return R_NilValue; } + void popedSolve(int &id) { rx_solving_options *op = getSolvingOptions(rx); rx_solving_options_ind *ind = getSolvingOptionsInd(rx, id); diff --git a/vignettes/articles/PopED.Rmd b/vignettes/articles/PopED.Rmd index b582a6e4..9afcc993 100644 --- a/vignettes/articles/PopED.Rmd +++ b/vignettes/articles/PopED.Rmd @@ -186,6 +186,7 @@ create the event table like a typical `rxode2` simulation, but it is used to specify the study design: ```{r babelmixr2-events} + library(babelmixr2) library(PopED) @@ -226,8 +227,7 @@ Other things you may have to include in your `PopED` model data frame are: the `PopED` database as `G_xt` - `id` becomes an ID for a design (which you can use as a covariate to - pool different designs or different). - + pool different designs or different regimens for optimal design). Once the design is setup, we need to specify a model. It is easy to specify the model using the `nlmixr2`/`rxode2` function/ui below: @@ -284,11 +284,9 @@ what is being added: This is the function that is run to generate the predictions: - ```{r} -# This can be retrieved directly from the database as: -# (Though it can also be viewed with f$popedFfFun) -poped_db_ode_babelmixr2$ff_fun +# The ff_fun can be retrieved from the ui with f$popedFfFun +f$popedFfFun ``` Some things to note in this function: @@ -378,13 +376,49 @@ f$popedD # PopED d f$popedNotfixedD # PopED notfixed_d f$popedCovd # PopED covd f$popedNotfixedCovd # PopED notfixed_covd -f$popedSigma # PopED sigma +f$popedSigma # PopED sigma (not variance is exported, not SD) f$popedNotfixedSigma # PopED notfixed_sigma ``` The rest of the parameters are generated in conjunction with the `popedControl()`. +## linear comparment models in `babelmixr2` + +You can also specify the models using the `linCmt()` solutions as +below: + +```{r lincmt} +f2 <- function() { + ini({ + tV <- 72.8 + tKA <- 0.25 + tCL <- 3.75 + Favail <- fix(0.9) + eta.ka ~ 0.09 + eta.cl ~ 0.25 ^ 2 + eta.v ~ 0.09 + prop.sd <- sqrt(0.04) + add.sd <- fix(sqrt(5e-6)) + }) + model({ + ka <- tKA * exp(eta.ka) + v <- tV * exp(eta.v) + cl <- tCL * exp(eta.cl) + cp <- linCmt() + f(depot) <- DOSE + cp ~ add(add.sd) + prop(prop.sd) + }) +} + +poped_db_analytic_babelmixr2 <- nlmixr(f, e, + popedControl(a=list(c(DOSE=20), + c(DOSE=40)), + maxa=c(DOSE=200), + mina=c(DOSE=0))) +``` + + ## Comparing method to the speed of other methods ```{r compare} @@ -393,9 +427,11 @@ library(microbenchmark) compare <- microbenchmark( evaluate_design(poped_db_analytic), + evaluate_design(poped_db_analytic_babelmixr2), evaluate_design(poped_db_ode_babelmixr2), evaluate_design(poped_db_ode_mrg), evaluate_design(poped_db_ode_pkpdsim), + evaluate_design(poped_db_ode_rx), times = 100L)