Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

581 rxode2 covariate sync #582

Merged
merged 5 commits into from
Sep 12, 2023
Merged

581 rxode2 covariate sync #582

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
c972f6b
add failing covarite sync test
mattfidler Sep 12, 2023
8bdde06
Unify same time with par_solve
mattfidler Sep 12, 2023
fcc6462
drop < eps
mattfidler Sep 12, 2023
14c2e9f
Fix sorted update of par_ptr
mattfidler Sep 12, 2023
c1fc2d1
Update news
mattfidler Sep 12, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions NEWS.md
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading