From c972f6b8a9baadff062df814dfb22bcca22e4b15 Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Tue, 12 Sep 2023 08:33:59 -0500 Subject: [PATCH 1/5] add failing covarite sync test --- tests/testthat/test-cov.R | 140 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 140 insertions(+) diff --git a/tests/testthat/test-cov.R b/tests/testthat/test-cov.R index b6a9d04fd..e87926a66 100644 --- a/tests/testthat/test-cov.R +++ b/tests/testthat/test-cov.R @@ -1,4 +1,144 @@ rxTest({ + + test_that("sync covariates", { + + dat <- structure(list(ID = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, + 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, + 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, + 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, + 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, + 3L, 3L, 3L, 3L, 3L, 3L, 3L), + TIME = c(0, 0, 0, 0.5, 0.5, 1, 1, + 1.5, 1.5, 2, 2, 2.5, 2.5, 3, 3, 3.5, 3.5, 4, 4, 4.5, 4.5, 5, + 5, 5.5, 5.5, 6, 6, 0, 0, 0, 0.5, 0.5, 1, 1, 1.5, 1.5, 2, 2, 2.5, + 2.5, 3, 3, 3.5, 3.5, 4, 4, 4.5, 4.5, 5, 5, 5.5, 5.5, 6, 6, 0, + 0, 0, 0.5, 0.5, 1, 1, 1.5, 1.5, 2, 2, 2.5, 2.5, 3, 3, 3.5, 3.5, + 4, 4, 4.5, 4.5, 5, 5, 5.5, 5.5, 6, 6), + EVID = c(80101L, 60101L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 80101L, 60101L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 80101L, 60101L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, + 0L, 0L, 0L), + AMT = c(10, 10, NA, NA, NA, NA, NA, NA, NA, NA, + NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, + NA, 10, 10, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, + NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 10, 10, NA, NA, + NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, + NA, NA, NA, NA, NA, NA, NA), + II = c(0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0), + DV = c(NA, NA, 0, 3.8, 0.9, 5.8, + 2.5, 3.2, 3.2, 1.8, 2.9, 1, 2.2, 0.5, 1.6, 0.3, 1.1, 0.2, 0.8, + 0.1, 0.5, 0, 0.3, 0, 0.2, 0, 0.1, NA, NA, 0, 3.4, 0.7, 7.3, 2.4, + 2.8, 2.5, 1.3, 2.7, 1, 1.6, 0.6, 1.1, 0.3, 0.8, 0.2, 0.8, 0.1, + 0.5, 0, 0.3, 0, 0.2, 0, 0.2, NA, NA, 0, 2.8, 1.1, 6.7, 2.7, 3.9, + 3.6, 1.5, 3.1, 0.7, 2.4, 0.5, 1.5, 0.3, 1.2, 0.1, 0.5, 0.1, 0.6, + 0.1, 0.4, 0, 0.2, 0, 0.2), + CMT = c(1L, 1L, 4L, 3L, 4L, 3L, 4L, + 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, + 3L, 4L, 3L, 4L, 1L, 1L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, + 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 1L, + 1L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, + 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L, 3L, 4L), + nlmixrRowNums = + c(1, + NA, 14, 2, 15, 3, 16, 4, 17, 5, 18, 6, 19, 7, 20, 8, 21, 9, 22, + 10, 23, 11, 24, 12, 25, 13, 26, 27, NA, 40, 28, 41, 29, 42, 30, + 43, 31, 44, 32, 45, 33, 46, 34, 47, 35, 48, 36, 49, 37, 50, 38, + 51, 39, 52, 53, NA, 66, 54, 67, 55, 68, 56, 69, 57, 70, 58, 71, + 59, 72, 60, 73, 61, 74, 62, 75, 63, 76, 64, 77, 65, 78)), + class = "data.frame", row.names = c(NA, -81L)) + + par <- structure(c(0, 0, 0, 0.693147180559945, 0.693147180559945, 0.693147180559945, + 2.30258509299405, 2.30258509299405, 2.30258509299405, 2.30258509299405, + 2.30258509299405, 2.30258509299405, 0, 0, 0, 1, 1, 1, 1, 1, 1, + -1.86782087447842, 0, 0, -2.44813462994314, 0, 0, 0.0831313226164289, + 0, 0, -2.28843604434224, 0, 0, 1.97816823869187e-05, 0, 0), + dim = c(3L, 12L), + dimnames = list(NULL, + c("THETA[1]", "THETA[2]", "THETA[3]", + "THETA[4]", "THETA[5]", "THETA[6]", "THETA[7]", "ETA[1]", "ETA[2]", + "ETA[3]", "ETA[4]", "ETA[5]"))) + + + rx <- rxode2({ + param(THETA[1], THETA[2], THETA[3], THETA[4], THETA[5], THETA[6], + THETA[7], ETA[1], ETA[2], ETA[3], ETA[4], ETA[5]) + cmt(center) + cmt(meta) + rx_expr_5 ~ ETA[3] + THETA[1] + rx_expr_6 ~ ETA[1] + THETA[2] + rx_expr_7 ~ ETA[2] + THETA[4] + rx_expr_10 ~ exp(rx_expr_5) + d/dt(center) = -rx_expr_10 * center - exp(rx_expr_6 - (rx_expr_7)) * + center + rx_expr_8 ~ ETA[5] + THETA[5] + rx_expr_11 ~ exp(rx_expr_8) + dur(center) = rx_expr_11 + rx_expr_9 ~ ETA[4] + THETA[3] + d/dt(meta) = rx_expr_10 * center - exp(rx_expr_9 - (rx_expr_7)) * + meta + center(0) = 0 + meta(0) = 0 + rx_expr_0 ~ CMT == 4 + rx_expr_1 ~ CMT == 3 + rx_expr_2 ~ 1 - (rx_expr_0) + rx_yj_ ~ 2 * (rx_expr_2) * (rx_expr_1) + 2 * (rx_expr_0) + rx_expr_3 ~ (rx_expr_0) + rx_expr_4 ~ (rx_expr_2) + rx_expr_12 ~ rx_expr_4 * (rx_expr_1) + rx_lambda_ ~ rx_expr_12 + rx_expr_3 + rx_hi_ ~ rx_expr_12 + rx_expr_3 + rx_low_ ~ 0 + rx_expr_13 ~ exp(-(rx_expr_7)) + rx_expr_14 ~ rx_expr_13 * meta + rx_expr_15 ~ rx_expr_13 * center + rx_expr_17 ~ rx_expr_14 * (rx_expr_0) + rx_expr_18 ~ rx_expr_15 * (rx_expr_2) + rx_expr_19 ~ rx_expr_18 * (rx_expr_1) + rx_pred_ = (rx_expr_0) * (rx_expr_17 + rx_expr_19) + rx_expr_18 * + Rx_pow_di((rx_expr_1), 2) + rx_r_ = (rx_expr_0) * Rx_pow_di(((rx_expr_17 + rx_expr_19) * + THETA[7]), 2) + (rx_expr_2) * Rx_pow_di((rx_expr_15 * + THETA[6] * (rx_expr_1)), 2) * (rx_expr_1) + tkm = THETA[1] + tcl = THETA[2] + tclm = THETA[3] + tv = THETA[4] + tdur0 = THETA[5] + prop.err = THETA[6] + prop.err2 = THETA[7] + eta.cl = ETA[1] + eta.v = ETA[2] + eta.km = ETA[3] + eta.clm = ETA[4] + eta.dur0 = ETA[5] + km = rx_expr_10 + cl = exp(rx_expr_6) + clm = exp(rx_expr_9) + v = exp(rx_expr_7) + dur0 = rx_expr_11 + cp = rx_expr_15 + cm = rx_expr_14 + tad = tad() + dosenum = dosenum() + cmt(cp) + cmt(cm) + dvid(3, 4) + }) + + s <- rxSolve(rx, dat, par, returnType="data.frame", + keep=c("nlmixrRowNums", "DV"), subsetNonmem=TRUE, addCov=TRUE) + + expect_equal(sort(unique(s[s$time==0.5,"CMT"])), 3:4) + + }) + skip_if_not_installed("units") for (meth in c("liblsoda", "lsoda")) { ## Dop is very close but doesn't match precisely. From 8bdde065c3a7d04d23ecbc88c053c867d2d81b5e Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Tue, 12 Sep 2023 08:42:31 -0500 Subject: [PATCH 2/5] Unify same time with par_solve --- src/approx.c | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/approx.c b/src/approx.c index 90bad93a5..3a53ce9f2 100644 --- a/src/approx.c +++ b/src/approx.c @@ -14,6 +14,8 @@ #define _as_zero(a) (fabs(a) < sqrt(DBL_EPSILON) ? 0.0 : a) #define _as_dbleps(a) (fabs(a) < sqrt(DBL_EPSILON) ? ((a) < 0 ? -sqrt(DBL_EPSILON) : sqrt(DBL_EPSILON)) : a) +#define isSameTimeOp(xout, xp) (op->stiff == 0 ? isSameTimeDop(xout, xp) : isSameTime(xout, xp)) + #ifdef ENABLE_NLS #include #define _(String) dgettext ("rxode2", String) @@ -208,9 +210,9 @@ void _update_par_ptr(double t, unsigned int id, rx_solve *rx, int idxIn) { } } // Pick best match - if (v == getTime(j, ind)) { + if (isSameTimeOp(v, getTime(j, ind))) { idx = j; - } else if (v == getTime(i, ind)) { + } else if (isSameTimeOp(v, getTime(i, ind))) { idx = i; } else if (op->is_locf == 2) { // nocb @@ -229,7 +231,6 @@ void _update_par_ptr(double t, unsigned int id, rx_solve *rx, int idxIn) { ind->_update_par_ptr_in = 1; if (ISNA(t)) { // functional lag, rate, duration, mtime - // Update all covariate parameters int k, idxSample; int ncov = op->ncov; @@ -252,7 +253,8 @@ void _update_par_ptr(double t, unsigned int id, rx_solve *rx, int idxIn) { ind->par_ptr[op->par_cov[k]-1] = getValue(idxSample, y, indSample, op); if (idx == 0){ ind->cacheME=0; - } else if (getValue(idxSample, y, indSample, op) != getValue(idxSample-1, y, indSample, op)) { + } else if (!isSameTimeOp(getValue(idxSample, y, indSample, op), + getValue(idxSample-1, y, indSample, op))) { ind->cacheME=0; } } @@ -279,13 +281,15 @@ void _update_par_ptr(double t, unsigned int id, rx_solve *rx, int idxIn) { double *par_ptr = ind->par_ptr; //double *all_times = indSample->all_times; double *y = indSample->cov_ptr + indSample->n_all_times*k; - if (idxSample == 0 && fabs(t- getAllTimes(indSample, idxSample)) < DBL_EPSILON) { + if (idxSample == 0 && + isSameTimeOp(t, getAllTimes(indSample, idxSample)) < DBL_EPSILON) { par_ptr[op->par_cov[k]-1] = y[0]; ind->cacheME=0; } else if (idxSample > 0 && idxSample < indSample->n_all_times && - fabs(t- getAllTimes(indSample, idxSample)) < DBL_EPSILON) { + isSameTimeOp(t, getAllTimes(indSample, idxSample))) { par_ptr[op->par_cov[k]-1] = getValue(idxSample, y, indSample, op); - if (getValue(idxSample, y, indSample, op) != getValue(idxSample-1, y, indSample, op)) { + if (!isSameTimeOp(getValue(idxSample, y, indSample, op), + getValue(idxSample-1, y, indSample, op))) { ind->cacheME=0; } } else { From fcc646208194c99a23675861e20deb8bfa1219ab Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Tue, 12 Sep 2023 09:26:59 -0500 Subject: [PATCH 3/5] drop < eps --- src/approx.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/approx.c b/src/approx.c index 3a53ce9f2..54cb55948 100644 --- a/src/approx.c +++ b/src/approx.c @@ -282,7 +282,7 @@ void _update_par_ptr(double t, unsigned int id, rx_solve *rx, int idxIn) { //double *all_times = indSample->all_times; double *y = indSample->cov_ptr + indSample->n_all_times*k; if (idxSample == 0 && - isSameTimeOp(t, getAllTimes(indSample, idxSample)) < DBL_EPSILON) { + isSameTimeOp(t, getAllTimes(indSample, idxSample))) { par_ptr[op->par_cov[k]-1] = y[0]; ind->cacheME=0; } else if (idxSample > 0 && idxSample < indSample->n_all_times && From 14c2e9f958cfd3aa115e9d6cb2b951cce313d096 Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Tue, 12 Sep 2023 09:41:00 -0500 Subject: [PATCH 4/5] Fix sorted update of par_ptr --- src/approx.c | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/approx.c b/src/approx.c index 54cb55948..0d9c3c993 100644 --- a/src/approx.c +++ b/src/approx.c @@ -191,28 +191,28 @@ void _update_par_ptr(double t, unsigned int id, rx_solve *rx, int idxIn) { idx = -1-ind->extraDoseTimeIdx[ind->extraDoseN[0]-1]; } // extra dose time, find the closest index - double v = getAllTimes(ind, idxIn); + double v = getTime(idxIn, ind); int i, j, ij, n = ind->n_all_times; i = 0; j = n - 1; - if (v < getTime(i, ind)) { + if (v < getTime(ind->ix[i], ind)) { idx = i; - } else if (v > getTime(j, ind)) { + } else if (v > getTime(ind->ix[j], ind)) { idx = j; } else { /* find the correct interval by bisection */ while(i < j - 1) { /* T(i) <= v <= T(j) */ ij = (i + j)/2; /* i+1 <= ij <= j-1 */ - if (v < getTime(ij, ind)) { + if (v < getTime(ind->ix[ij], ind)) { j = ij; } else { i = ij; } } // Pick best match - if (isSameTimeOp(v, getTime(j, ind))) { + if (isSameTimeOp(v, getTime(ind->ix[j], ind))) { idx = j; - } else if (isSameTimeOp(v, getTime(i, ind))) { + } else if (isSameTimeOp(v, getTime(ind->ix[i], ind))) { idx = i; } else if (op->is_locf == 2) { // nocb @@ -282,11 +282,11 @@ void _update_par_ptr(double t, unsigned int id, rx_solve *rx, int idxIn) { //double *all_times = indSample->all_times; double *y = indSample->cov_ptr + indSample->n_all_times*k; if (idxSample == 0 && - isSameTimeOp(t, getAllTimes(indSample, idxSample))) { + isSameTimeOp(t, getTime(ind->ix[idxSample], indSample))) { par_ptr[op->par_cov[k]-1] = y[0]; ind->cacheME=0; } else if (idxSample > 0 && idxSample < indSample->n_all_times && - isSameTimeOp(t, getAllTimes(indSample, idxSample))) { + isSameTimeOp(t, getTime(ind->ix[idxSample], indSample))) { par_ptr[op->par_cov[k]-1] = getValue(idxSample, y, indSample, op); if (!isSameTimeOp(getValue(idxSample, y, indSample, op), getValue(idxSample-1, y, indSample, op))) { From c1fc2d170aa8460c5d4b4cf8f34eb189656780b1 Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Tue, 12 Sep 2023 09:59:37 -0500 Subject: [PATCH 5/5] Update news --- NEWS.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/NEWS.md b/NEWS.md index d70bf061a..825e20cb9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -128,6 +128,9 @@ mu-referencing style to run the optimization. - Bug fix for `geom_amt()` when the `aes` transformation has `x` +- Bug fix for some covariate updates that may affect multiple compartment + models (like issue #581) + # rxode2 2.0.13 ## Bug fixes