Skip to content

Commit

Permalink
Merge pull request #582 from nlmixr2/581-rxode2-covariate-sync
Browse files Browse the repository at this point in the history 
581 rxode2 covariate sync
  • Loading branch information
@mattfidler
mattfidler authored Sep 12, 2023
2 parents 49c3593 + c1fc2d1 commit becaef1
Show file tree
Hide file tree
Showing 3 changed files with 158 additions and 11 deletions.
3 changes: 3 additions & 0 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
26 changes: 15 additions & 11 deletions src/approx.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <libintl.h>
#define _(String) dgettext ("rxode2", String)
Expand Down Expand Up @@ -189,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 (v == getTime(j, ind)) {
if (isSameTimeOp(v, getTime(ind->ix[j], ind))) {
idx = j;
} else if (v == getTime(i, ind)) {
} else if (isSameTimeOp(v, getTime(ind->ix[i], ind))) {
idx = i;
} else if (op->is_locf == 2) {
// nocb
Expand All @@ -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;
Expand All @@ -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;
}
}
Expand All @@ -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, 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 &&
fabs(t- getAllTimes(indSample, idxSample)) < DBL_EPSILON) {
isSameTimeOp(t, getTime(ind->ix[idxSample], indSample))) {
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 {
Expand Down
140 changes: 140 additions & 0 deletions tests/testthat/test-cov.R
View file
Original file line number Diff line number Diff line change
@@ -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.
Expand Down

0 comments on commit becaef1

Please sign in to comment.