Skip to content

Commit

Permalink
Merge pull request #671 from nlmixr2/670-a-different-model-with-est-=…
Browse files Browse the repository at this point in the history 
…-rxsolve-crashes-r

670 a different model with est = rxsolve crashes r
  • Loading branch information
@mattfidler
mattfidler authored Apr 7, 2024
2 parents 316b236 + e7698f0 commit 9042f6a
Show file tree
Hide file tree
Showing 3 changed files with 211 additions and 2 deletions.
6 changes: 5 additions & 1 deletion NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,11 @@

## Bug fixes

- Make sure that the object is a uncompressed rxode2 ui for solving with `rxSolve` (See #661)
- Make sure that the object is a uncompressed rxode2 ui for solving with `rxSolve` (See #661)

- Fix #670 by using the last simulated observation residual when there
are trailing doses.


## New features

Expand Down
2 changes: 1 addition & 1 deletion src/rxode2_df.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -360,7 +360,7 @@ extern "C" SEXP rxode2_df(int doDose0, int doTBS) {
if ((doDose && evid!= 9) || (evid0 == 0 && isObs(evid)) || (evid0 == 1 && evid==0)){
// Only increment if this is an observation or of this a
// simulation that requests dosing information too.
kk++;
kk=min2(kk+1, errNrow-1);
}
}
if (nlhs){
Expand Down
205 changes: 205 additions & 0 deletions tests/testthat/test-sim-err-670.R
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,205 @@
test_that("Simulate Error, issue #670", {

m <- function() {
ini({
lfdose_max <- -0.1
lfdose_50 <- 1
lktr <- 0.1
lka <- -1
lcl <- 1.1
lvc <- 1.1
lvp <- 1.1
lq <- 1.1
cl_covar <- -1
propSd <- c(0, 0.1)
addSd <- c(0, 1)
iiv_ka ~ 0.1
iiv_cl ~ 0.1
iiv_vc ~ 0.1
})
model({
fdose_max <- exp(lfdose_max)
fdose_50 <- exp(lfdose_50)
fdose <- 1 + fdose_max * DOSE/(fdose_50 + DOSE)
fdepot <- fdose
ktr <- exp(lktr)
ka <- exp(lka + iiv_ka)
cl <- exp(lcl + iiv_cl + COVAR * cl_covar)
vc <- exp(lvc + iiv_vc)
vp <- exp(lvp)
q <- exp(lq)
kel <- cl/vc
k12 <- q/vc
k21 <- q/vp
d/dt(depot) <- -ktr * depot
d/dt(transit) <- -ka * transit + ktr * fdepot * depot
d/dt(central) <- ka * transit - kel * central - k12 *
central + k21 * peripheral1
d/dt(peripheral1) <- k12 * central - k21 * peripheral1
cc <- 1000 * central/vc
cc ~ prop(propSd) + add(addSd)
})
}

m <- m()

d <-
structure(list(ID = c("A131", "A131", "A131", "A131", "A131",
"A131", "A131", "A131", "A131", "A131", "A131", "A131", "A131",
"A131", "A131", "A131", "A131", "A131", "A131", "A131", "A131",
"A131", "A131", "A131", "A131", "A131", "A131", "A131", "A131",
"A131", "A131", "A131", "A131", "A131", "A131", "A131", "A131",
"A131", "A131", "A131", "A131", "A131", "A131", "A131", "A131",
"A131", "A131", "A131", "A131", "A131", "A131", "A131", "A131",
"A131", "A131", "A131", "A131", "A131", "A131", "A131", "A131",
"A131", "A131", "A131", "A131", "A131", "A131", "A131", "A131",
"A131", "A131", "A131", "A131", "A131", "A131", "A131", "A131",
"A131", "A131", "A131", "A131", "A131", "A131", "A131", "A131",
"A131", "A131", "A131", "A131", "A131", "A131", "A131", "A131",
"A131", "A131", "A131", "A131", "A132", "A132", "A132", "A132",
"A132", "A132", "A132", "A132", "A132", "A132", "A132", "A132",
"A132", "A132", "A132", "A132", "A132", "A132", "A132", "A132",
"A132", "A132", "A132", "A132", "A132", "A132", "A132", "A132",
"A132", "A132", "A132", "A132", "A132", "A132", "A132", "A132",
"A132", "A132", "A132", "A132", "A133", "A133", "A133", "A133",
"A133", "A133", "A133", "A133", "A133", "A133", "A133", "A133",
"A133", "A133", "A133", "A133", "A133", "A133", "A133", "A133",
"A133", "A133", "A133", "A133", "A133", "A133", "A133", "A133",
"A133", "A133", "A133", "A133", "A133", "A133", "A133", "A133",
"A133", "A133", "A133", "A134", "A134", "A134", "A134", "A134",
"A134", "A134", "A134", "A134", "A134", "A134", "A135", "A135",
"A135", "A135", "A135", "A135", "A135", "A135", "A135", "A135",
"A135", "A135", "A135", "A135", "A135", "A135", "A135", "A135",
"A135", "A135", "A135", "A135", "A135", "A135", "A135", "A135",
"A135", "A135", "A135", "A135", "A135"),
TIME = c(0, 0, 1.3,
1.883333333, 2.766666667, 4.85, 6.883333333, 7.816666667, 23.81666667,
47.81666667, 71.81666667, 95.81666667, 166.2, 166.2, 190.2, 214.2,
238.2, 262.2, 332.3333333, 332.3333333, 332.95, 333.5333333,
334.5333333, 336.5333333, 338.5333333, 339.45, 355.45, 379.45,
403.45, 427.45, 504.3166667, 504.3166667, 528.3166667, 552.3166667,
576.3166667, 672.3166667, 696.3166667, 720.3166667, 744.3166667,
840.3166667, 864.3166667, 888.3166667, 909.5333333, 909.5333333,
910.3666667, 910.8166667, 911.7, 913.45, 915.4, 917, 1009.783333,
1009.783333, 1012.3, 1034.3, 1058.3, 1082.3, 1178.3, 1202.3,
1226.3, 1250.3, 2016.1, 2016.1, 2040.1, 2064.1, 2088.1, 2184.1,
2208.1, 2232.1, 2256.1, 2352.1, 2376.1, 2400.1, 2424.1, 2522.1,
2522.1, 2522.45, 2544.45, 2568.45, 2592.45, 2688.45, 2712.45,
2736.45, 2760.45, 2856.45, 2880.45, 2904.45, 2928.45, 3025.8,
3025.8, 3028.266667, 3050.266667, 3074.266667, 3098.266667, 3194.266667,
3218.266667, 3242.266667, 3266.266667, 0, 0, 0.85, 24.35, 48.35,
72.35, 96.35, 336.35, 360.35, 384.35, 408.35, 432.35, 500.7,
500.7, 503.0166667, 525.0166667, 549.0166667, 573.0166667, 597.0166667,
666.1333333, 666.1333333, 666.9833333, 667.4833333, 668.4833333,
669.9833333, 672.0333333, 673.7, 689.7, 713.7, 737.7, 761.7,
833.7, 881.7, 905.7, 1001.7, 1025.7, 1049.7, 1082.816667, 1082.816667,
1085.7, 0, 0, 1.616666667, 2.033333333, 2.866666667, 4.866666667,
6.283333333, 24.28333333, 170.2833333, 170.2833333, 194.2833333,
218.2833333, 242.2833333, 261.6833333, 333.1666667, 357.1666667,
381.1666667, 405.1666667, 429.1666667, 505.35, 505.35, 508.3333333,
530.3333333, 554.3333333, 578.3333333, 595.4, 595.4, 596.9333333,
666.9333333, 690.9333333, 714.9333333, 738.9333333, NA, NA, NA,
NA, NA, NA, NA, 0, 0, 1.133333333, 1.716666667, 2.633333333,
4.416666667, 6.216666667, 24.21666667, 48.21666667, 72.21666667,
96.21666667, 0, 0, 0.9666666667, 1.466666667, 2.266666667, 4.35,
6.016666667, 24.01666667, 48.01666667, 168.0166667, 192.0166667,
216.0166667, 240.0166667, 264.0166667, 336.0166667, 360.0166667,
384.0166667, 408.0166667, 432.0166667, 503.3, 503.3, 505.9, 527.9,
551.9, 575.9, 599.9, 1175.566667, 1175.566667, 1178, 1200, 1224),
CMT = c("central", "depot", "central", "central", "central",
"central", "central", "central", "depot", "depot", "depot", "depot",
"central", "depot", "depot", "depot", "depot", "depot", "central",
"depot", "central", "central", "central", "central", "central",
"central", "depot", "depot", "depot", "depot", "central", "depot",
"depot", "depot", "depot", "depot", "depot", "depot", "depot",
"depot", "depot", "depot", "central", "depot", "central", "central",
"central", "central", "central", "central", "central", "depot",
"central", "depot", "depot", "depot", "depot", "depot", "depot",
"depot", "central", "depot", "depot", "depot", "depot", "depot",
"depot", "depot", "depot", "depot", "depot", "depot", "depot",
"central", "depot", "central", "depot", "depot", "depot", "depot",
"depot", "depot", "depot", "depot", "depot", "depot", "depot",
"central", "depot", "central", "depot", "depot", "depot", "depot",
"depot", "depot", "depot", "central", "depot", "central", "depot",
"depot", "depot", "depot", "depot", "depot", "depot", "depot",
"depot", "central", "depot", "central", "depot", "depot", "depot",
"depot", "central", "depot", "central", "central", "central",
"central", "central", "central", "depot", "depot", "depot", "depot",
"depot", "depot", "depot", "depot", "depot", "depot", "central",
"depot", "central", "central", "depot", "central", "central",
"central", "central", "central", "depot", "central", "depot",
"depot", "depot", "depot", "depot", "depot", "depot", "depot",
"depot", "depot", "central", "depot", "central", "depot", "depot",
"depot", "central", "depot", "central", "depot", "depot", "depot",
"depot", "central", "central", "central", "central", "central",
"central", "central", "central", "depot", "central", "central",
"central", "central", "central", "depot", "depot", "depot", "depot",
"central", "depot", "central", "central", "central", "central",
"central", "depot", "depot", "depot", "depot", "depot", "depot",
"depot", "depot", "depot", "depot", "depot", "depot", "central",
"depot", "central", "depot", "depot", "depot", "depot", "central",
"depot", "central", "depot", "depot"),
AMT = c(NA, 400, NA, NA,
NA, NA, NA, NA, 400, 400, 400, 400, NA, 400, 400, 400, 400, 400,
NA, 400, NA, NA, NA, NA, NA, NA, 400, 400, 400, 400, NA, 400,
400, 400, 400, 400, 400, 400, 400, 400, 400, 400, NA, 400, NA,
NA, NA, NA, NA, NA, NA, 400, NA, 400, 400, 400, 400, 400, 400,
400, NA, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400,
400, NA, 400, NA, 400, 400, 400, 400, 400, 400, 400, 400, 400,
400, 400, NA, 400, NA, 400, 400, 400, 400, 400, 400, 400, NA,
500, NA, 500, 500, 500, 500, 450, 450, 450, 450, 450, NA, 450,
NA, 450, 450, 450, 450, NA, 450, NA, NA, NA, NA, NA, NA, 450,
450, 450, 450, 450, 450, 450, 450, 450, 450, NA, 450, NA, NA,
350, NA, NA, NA, NA, NA, 350, NA, 300, 300, 300, 300, 300, 300,
300, 300, 300, 300, NA, 300, NA, 300, 300, 300, NA, 300, NA,
300, 300, 300, 300, NA, NA, NA, NA, NA, NA, NA, NA, 350, NA,
NA, NA, NA, NA, 350, 350, 350, 350, NA, 400, NA, NA, NA, NA,
NA, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400,
NA, 400, NA, 400, 400, 400, 400, NA, 400, NA, 400, 400),
DOSE = c(0,
400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400,
400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400,
400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400,
400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400,
400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400,
400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400,
400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400,
400, 400, 400, 400, 400, 0, 500, 500, 500, 500, 500, 500, 450,
450, 450, 450, 450, 450, 450, 450, 450, 450, 450, 450, 450, 450,
450, 450, 450, 450, 450, 450, 450, 450, 450, 450, 450, 450, 450,
450, 450, 450, 450, 450, 450, 0, 350, 350, 350, 350, 350, 350,
350, 350, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300,
300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300,
300, 300, 300, 300, 300, 300, 0, 350, 350, 350, 350, 350, 350,
350, 350, 350, 350, 0, 400, 400, 400, 400, 400, 400, 400, 400,
400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400, 400,
400, 400, 400, 400, 400, 400, 400, 400, 400),
EVID = c(0, 1,
0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0,
0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1,
0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1,
1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0,
0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1,
1, 0, 1, 0, 1, 1),
COVAR = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1,
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, 1, 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,
1, 1, 1, 1, 1, 1, 1, 1, 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)),
row.names = c(NA, -218L),
class = "data.frame")

expect_error(rxSolve(object = m, events = d), NA)

})

0 comments on commit 9042f6a

Please sign in to comment.