From b06c06e5b6fe9b3b2b662e7d86ff768212a876e0 Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Thu, 17 Oct 2024 16:00:39 -0500 Subject: [PATCH 1/5] Assign global poped environment to assess what model is loaded --- R/RcppExports.R | 4 ++-- R/poped.R | 13 +++++++------ src/RcppExports.cpp | 7 ++++--- src/init.c | 4 ++-- src/poped.cpp | 8 ++++++-- 5 files changed, 21 insertions(+), 15 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 97659764..29584c78 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -145,8 +145,8 @@ popedFree <- function() { .Call(`_babelmixr2_popedFree`) } -popedSetup <- function(e, full) { - .Call(`_babelmixr2_popedSetup`, e, full) +popedSetup <- function(e, eglobal, full) { + .Call(`_babelmixr2_popedSetup`, e, eglobal, full) } popedSolveIdME <- function(theta, umt, mt, ms, nend, id, totn) { diff --git a/R/poped.R b/R/poped.R index fd06dc0a..742e9c7d 100644 --- a/R/poped.R +++ b/R/poped.R @@ -17,13 +17,13 @@ #' This should not typically be called directly #' #' @param e environment with setup information for popEd +#' @param eglobal global environment for poped info #' @param full setup the full model #' @return nothing, called for side effects -#' @export #' @keywords internal #' @author Matthew L. Fidler -.popedSetup <- function(e, full=FALSE) { - invisible(.Call(`_babelmixr2_popedSetup`, e, full)) +.popedSetup <- function(e, eglobal, full=FALSE) { + invisible(.Call(`_babelmixr2_popedSetup`, e, eglobal, full)) } #' Solve poped problem for appropriate times with single/multiple endpoint models #' @@ -289,6 +289,7 @@ rxUiGet.popedFgFun <- function(x, ...) { attr(rxUiGet.popedFgFun, "desc") <- "PopED parameter model (fg_fun)" .poped <- new.env(parent=emptyenv()) +.poped$curLib <- NULL .poped$s <- NULL .poped$modelNumber <- 1L .poped$curNumber <- -1L @@ -418,7 +419,7 @@ attr(rxUiGet.popedFfFun, "desc") <- "PopED parameter model (ff_fun)" } if (.poped$setup != 1L) { rxode2::rxSolveFree() - .popedSetup(popedDb$babelmixr2, FALSE) + .popedSetup(popedDb$babelmixr2, .poped, FALSE) .poped$curNumber <- popedDb$babelmixr2$modelNumber .poped$setup <- 1L .poped$fullXt <- NULL @@ -488,7 +489,7 @@ attr(rxUiGet.popedFfFun, "desc") <- "PopED parameter model (ff_fun)" }))) .et <- rxode2::etTrans(.dat, .e$modelF) .e$dataF <- .et - .popedSetup(.e, TRUE) + .popedSetup(.e, .poped, TRUE) .poped$fullXt <- length(xt) .poped$curNumber <- popedDb$babelmixr2$modelNumber .poped$setup <- 2L @@ -554,7 +555,7 @@ attr(rxUiGet.popedFfFun, "desc") <- "PopED parameter model (ff_fun)" }))) .et <- rxode2::etTrans(.dat, .e$modelF) .e$dataF <- .et - .popedSetup(.e, TRUE) + .popedSetup(.e, .poped, TRUE) .poped$fullXt <- length(xt) .poped$curNumber <- popedDb$babelmixr2$modelNumber .poped$setup <- 2L diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index abe37bda..1fff9d19 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -94,14 +94,15 @@ BEGIN_RCPP END_RCPP } // popedSetup -RObject popedSetup(Environment e, bool full); -RcppExport SEXP _babelmixr2_popedSetup(SEXP eSEXP, SEXP fullSEXP) { +RObject popedSetup(Environment e, Environment eglobal, bool full); +RcppExport SEXP _babelmixr2_popedSetup(SEXP eSEXP, SEXP eglobalSEXP, SEXP fullSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< Environment >::type e(eSEXP); + Rcpp::traits::input_parameter< Environment >::type eglobal(eglobalSEXP); Rcpp::traits::input_parameter< bool >::type full(fullSEXP); - rcpp_result_gen = Rcpp::wrap(popedSetup(e, full)); + rcpp_result_gen = Rcpp::wrap(popedSetup(e, eglobal, full)); return rcpp_result_gen; END_RCPP } diff --git a/src/init.c b/src/init.c index f26569b5..4601685e 100644 --- a/src/init.c +++ b/src/init.c @@ -13,7 +13,7 @@ SEXP _babelmixr2_transDv(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP _babelmixr2_iniRxodePtrs(SEXP in); SEXP _babelmixr2_popedFree(void); -SEXP _babelmixr2_popedSetup(SEXP eSEXP, SEXP fullSEXP); +SEXP _babelmixr2_popedSetup(SEXP eSEXP, SEXP e2SEXP, SEXP fullSEXP); SEXP _babelmixr2_popedSolveIdME(SEXP thetaSEXP, SEXP umtSEXP, SEXP mtSEXP, SEXP msSEXP, SEXP nendSEXP, SEXP idSEXP, SEXP totnSEXP); SEXP _babelmixr2_popedSolveIdME2(SEXP thetaSEXP, SEXP umtSEXP, SEXP mtSEXP, SEXP msSEXP, SEXP nendSEXP, SEXP idSEXP, SEXP totnSEXP); SEXP _babelmixr2_popedGetMultipleEndpointModelingTimes(SEXP, SEXP, SEXP); @@ -29,7 +29,7 @@ static const R_CallMethodDef CallEntries[] = { {"_babelmixr2_popedGetMultipleEndpointModelingTimes", (DL_FUNC) &_babelmixr2_popedGetMultipleEndpointModelingTimes, 3}, {"_babelmixr2_popedFree", (DL_FUNC) &_babelmixr2_popedFree, 0}, - {"_babelmixr2_popedSetup", (DL_FUNC) &_babelmixr2_popedSetup, 2}, + {"_babelmixr2_popedSetup", (DL_FUNC) &_babelmixr2_popedSetup, 3}, {"_babelmixr2_popedSolveIdME", (DL_FUNC) &_babelmixr2_popedSolveIdME, 7}, {"_babelmixr2_popedSolveIdME2", (DL_FUNC) &_babelmixr2_popedSolveIdME2, 7}, {"_babelmixr2_iniRxodePtrs", (DL_FUNC) &_babelmixr2_iniRxodePtrs, 1}, diff --git a/src/poped.cpp b/src/poped.cpp index b3116229..0d37a4c1 100644 --- a/src/poped.cpp +++ b/src/poped.cpp @@ -23,6 +23,7 @@ timeIndexer globalTimeIndexer; + //' @title Get Multiple Endpoint Modeling Times //' //' @description @@ -247,6 +248,9 @@ using namespace Rcpp; #define popedOde(id) ind_solve(rx, id, rxInner.dydt_liblsoda, rxInner.dydt_lsoda_dum, rxInner.jdum_lsoda, rxInner.dydt, rxInner.update_inis, rxInner.global_jt) +Environment _popedE; +Environment _popedEglobal; + struct rxSolveF { // // std::string estStr; @@ -338,12 +342,12 @@ RObject popedFree() { return R_NilValue; } -Environment _popedE; //[[Rcpp::export]] -RObject popedSetup(Environment e, bool full) { +RObject popedSetup(Environment e, Environment eglobal, bool full) { popedFree(); _popedE=e; + _popedEglobal=eglobal; List control = e["control"]; List rxControl = as(e["rxControl"]); From 6058184a0c75dce82d8f66119a367dedd66a9ff2 Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Thu, 17 Oct 2024 19:10:18 -0500 Subject: [PATCH 2/5] Always assign model variables before running ode solving --- src/poped.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/poped.cpp b/src/poped.cpp index 0d37a4c1..9ca203fe 100644 --- a/src/poped.cpp +++ b/src/poped.cpp @@ -368,7 +368,9 @@ RObject popedSetup(Environment e, Environment eglobal, bool full) { e["paramCache"]=p2; e["lid"] = NA_INTEGER; List mvp = rxode2::rxModelVars_(model); - rxUpdateFuns(as(mvp["trans"]), &rxInner); + CharacterVector trans = mvp["trans"]; + _popedEglobal["curTrans"] = trans; + rxUpdateFuns(as(trans), &rxInner); // initial value of parameters CharacterVector pars = mvp[RxMv_params]; @@ -451,6 +453,7 @@ static inline bool solveCached(NumericVector &theta, int &id) { } void popedSolveFidMat(arma::mat &matMT, NumericVector &theta, int id, int nrow, int nend) { + rxUpdateFuns(_popedEglobal["curTrans"], &rxInner); // arma::vec ret(retD, nobs, false, true); rx_solving_options_ind *ind = updateParamRetInd(theta, id); rx_solving_options *op = getSolvingOptions(rx); @@ -459,12 +462,12 @@ void popedSolveFidMat(arma::mat &matMT, NumericVector &theta, int id, int nrow, int kk, k=0; double curT; bool isMT = false; + double *lhs = getIndLhs(ind); for (int j = 0; j < getIndNallTimes(ind); ++j) { setIndIdx(ind, j); kk = getIndIx(ind, j); curT = getTime(kk, ind); isMT = getIndEvid(ind, kk) >= 10 && getIndEvid(ind, kk) <= 99; - double *lhs = getIndLhs(ind); if (isDose(getIndEvid(ind, kk))) { rxInner.calc_lhs(id, curT, getOpIndSolve(op, ind, j), lhs); continue; @@ -500,6 +503,7 @@ Rcpp::DataFrame popedSolveIdME(NumericVector &theta, NumericVector &mt, IntegerVector &ms, int nend, int id, int totn) { if (solveCached(theta, id)) return(as(_popedE["s"])); + rxUpdateFuns(_popedEglobal["curTrans"], &rxInner); NumericVector t(totn); arma::vec f(totn); arma::vec w(totn); @@ -512,7 +516,6 @@ Rcpp::DataFrame popedSolveIdME(NumericVector &theta, std::fill(curLV.begin(), curLV.end(), 0); we[i] = curLV; } - popedSolveFidMat(matMT, theta, id, nrow, nend); // this gets the information from: // - the model time From 3f100abf749dfcef753a64ca9b5b1985d9d6b704 Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Thu, 17 Oct 2024 19:52:41 -0500 Subject: [PATCH 3/5] Reset time index whenever the model is setu --- R/poped.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/poped.R b/R/poped.R index 742e9c7d..5072a334 100644 --- a/R/poped.R +++ b/R/poped.R @@ -23,6 +23,7 @@ #' @keywords internal #' @author Matthew L. Fidler .popedSetup <- function(e, eglobal, full=FALSE) { + popedMultipleEndpointResetTimeIndex() invisible(.Call(`_babelmixr2_popedSetup`, e, eglobal, full)) } #' Solve poped problem for appropriate times with single/multiple endpoint models @@ -337,10 +338,10 @@ rxUiGet.popedFfFun <- function(x, ...) { # be in the right order if (.lu <= .(.poped$maxn)) { # only check for time reset if it is specified in the model + .popedRxRunSetup(poped.db) .p <- babelmixr2::popedMultipleEndpointParam(p, .u, model_switch, .(.poped$maxn), poped.db$babelmixr2$optTime) - .popedRxRunSetup(poped.db) .ret <- .popedSolveIdME(.p, .u, .xt, model_switch, .(length(.predDf$cond)), .id-1, .totn) } else if (.lu > .(.poped$maxn)) { From 62aa16f94dfb36c9e9db1f65dd8bbf436c74c636 Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Thu, 17 Oct 2024 20:10:58 -0500 Subject: [PATCH 4/5] Add test --- tests/testthat/test-poped.R | 122 ++++++++++++++++++++++++++++++++++++ 1 file changed, 122 insertions(+) diff --git a/tests/testthat/test-poped.R b/tests/testthat/test-poped.R index ce53ad11..07d88e01 100644 --- a/tests/testthat/test-poped.R +++ b/tests/testthat/test-poped.R @@ -316,6 +316,128 @@ if (requireNamespace("PopED", quietly=TRUE)) { }) + test_that("shrinkage", { + + f <- function() { + ini({ + tV <- 72.8 + tKa <- 0.25 + tCl <- 3.75 + tF <- fix(0.9) + eta.v ~ 0.09 + eta.ka ~ 0.09 + eta.cl ~0.25^2 + prop.sd <- fix(sqrt(0.04)) + add.sd <- fix(sqrt(5e-6)) + }) + model({ + V<-tV*exp(eta.v) + KA<-tKa*exp(eta.ka) + CL<-tCl*exp(eta.cl) + Favail <- tF + N <- floor(time/TAU)+1 + y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) * + (exp(-CL/V * (time - (N - 1) * TAU)) * + (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - + exp(-KA * (time - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) + y ~ prop(prop.sd) + add(add.sd) + }) + } + + # minxt, maxxt + e <- et(list(c(0, 10), + c(0, 10), + c(0, 10), + c(240, 248), + c(240, 248))) %>% + as.data.frame() + + #xt + e$time <- c(1,2,8,240,245) + + + babel.db <- nlmixr2(f, e, "poped", + popedControl(groupsize=20, + bUseGrouped_xt=TRUE, + a=list(c(DOSE=20,TAU=24), + c(DOSE=40, TAU=24)), + maxa=c(DOSE=200,TAU=24), + mina=c(DOSE=0,TAU=24))) + + expect_error(shrinkage(babel.db), NA) + + }) + + test_that("mixing 2 poped databases at the same time", { + library(PopED) + + # Define the model + f <- function() { + ini({ + tV <- 72.8 + tKa <- 0.25 + tCl <- 3.75 + tF <- fix(0.9) + pedCL <- 0.8 + + eta.v ~ 0.09 + eta.ka ~ 0.09 + eta.cl ~0.25^2 + + prop.sd <- fix(sqrt(0.04)) + add.sd <- fix(sqrt(5e-6)) + + }) + model({ + V<-tV*exp(eta.v) + KA<-tKa*exp(eta.ka) * (pedCL**isPediatric) # add covariate for pediatrics + CL<-tCl*exp(eta.cl) + Favail <- tF + + N <- floor(t/TAU)+1 + y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) * + (exp(-CL/V * (t - (N - 1) * TAU)) * + (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - + exp(-KA * (t - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) + + y ~ prop(prop.sd) + add(add.sd) + }) + } + + e <- et(c( 1,8,10,240,245)) + + babel.db <- nlmixr2(f, e, "poped", + popedControl(m = 2, + groupsize=20, + bUseGrouped_xt=TRUE, + a=list(c(DOSE=20,TAU=24,isPediatric = 0), + c(DOSE=40, TAU=24,isPediatric = 0)))) + + + # Note, to be able to use the adults FIM to combine with the pediatrics, + # both have to have the parameter "pedCL" defined and set notfixed_bpop to 1. + + ## Define pediatric model/design (isPediatric = 1) + ## One arm, 4 time points only + + e.ped <- et(c( 1,2,6,240)) + + babel.db.ped <- nlmixr2(f, e.ped, "poped", + popedControl(m = 1, + groupsize=6, + bUseGrouped_xt=TRUE, + a=list(c(DOSE=40,TAU=24,isPediatric = 1)))) + + + ## Create plot of model of adult data without variability + vdiffr::expect_doppelganger("poped-adult-first", + plot_model_prediction(babel.db, model_num_points = 300)) + + vdiffr::expect_doppelganger("poped-ped-next", + plot_model_prediction(babel.db.ped, + model_num_points = 300)) + }) + ## The tests run interactively runs OK ## However, the capture output seems to be interfere with the tests (from PopED) # So... these are commented out for now. From 033c09e721c6819fa8157bf8f2810919e284965e Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Thu, 17 Oct 2024 20:58:42 -0500 Subject: [PATCH 5/5] Update news --- NEWS.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/NEWS.md b/NEWS.md index 8fa3f358..9d0d5460 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,14 @@ * Fix bug where if/else clauses in the model could cause the model to not predict the values correctly. +* Fix bug so that `shrinkage()` calculation works + +* Fix bug so that you can mix 2 different `PopED` data bases in an + analysis without crashing R. While this didn't occur with every + database clash, it more frequently occurred when you interleaved + `PopED` code between two different `PopED` databases, like in issue + #131. + # babelmixr2 0.1.4 * Added experimental `PopED` integration