From 2991f241562c5ac04da49d33febd2fba61e493e1 Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Sat, 6 Apr 2024 23:18:47 -0500 Subject: [PATCH 1/6] Add crash as failing test --- tests/testthat/test-sim-err-670.R | 206 ++++++++++++++++++++++++++++++ 1 file changed, 206 insertions(+) create mode 100644 tests/testthat/test-sim-err-670.R diff --git a/tests/testthat/test-sim-err-670.R b/tests/testthat/test-sim-err-670.R new file mode 100644 index 000000000..fcf211e76 --- /dev/null +++ b/tests/testthat/test-sim-err-670.R @@ -0,0 +1,206 @@ +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") + + crsh <- rxSolve(object = m, + events = d) + +}) From 40e51ed9dc9fa9b7283c44ba73a22c17c808e0ea Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Sat, 6 Apr 2024 23:19:01 -0500 Subject: [PATCH 2/6] Add hack fix --- src/rxData.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/rxData.cpp b/src/rxData.cpp index f807444b8..fd8f772e8 100644 --- a/src/rxData.cpp +++ b/src/rxData.cpp @@ -3292,6 +3292,9 @@ static inline void rxSolve_simulate(const RObject &obj, if (!cbindPar1) { params0 = as>(rxSolveDat->par1); } + curObs++; // sometimes there is trailing non-observation evids + // (though this is an edge case, probably related to NA + // times; See #670) List lst = rxSimThetaOmega(params0, omega, omegaDf, From cb3ee65cc536574630f8cc095fab63d786a08ca6 Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Sat, 6 Apr 2024 23:39:51 -0500 Subject: [PATCH 3/6] Update news --- NEWS.md | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 602728cf5..6054b6726 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,7 +2,14 @@ ## 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) + +- Ensure the residual simulation includes one simulation above the + number of observations; this allows trailing doses to access the + extra simulated residuals without accessing memory that has not been + allocated. This fixes the R crashed observed in #670. This did no + happen often and likely does not affect much code since the + `etTrans()` typically drops trailing doses. ## New features From 67c22b574a72b20170876cc864926330e634f9b1 Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Sat, 6 Apr 2024 23:42:16 -0500 Subject: [PATCH 4/6] Revert "Add hack fix" This reverts commit 40e51ed9dc9fa9b7283c44ba73a22c17c808e0ea. --- src/rxData.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/rxData.cpp b/src/rxData.cpp index fd8f772e8..f807444b8 100644 --- a/src/rxData.cpp +++ b/src/rxData.cpp @@ -3292,9 +3292,6 @@ static inline void rxSolve_simulate(const RObject &obj, if (!cbindPar1) { params0 = as>(rxSolveDat->par1); } - curObs++; // sometimes there is trailing non-observation evids - // (though this is an edge case, probably related to NA - // times; See #670) List lst = rxSimThetaOmega(params0, omega, omegaDf, From d08eee8d762f2d9f73668b5655b15d5a98c05f6a Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Sat, 6 Apr 2024 23:46:47 -0500 Subject: [PATCH 5/6] Add fix that does not change the random numbers --- NEWS.md | 9 +++------ src/rxode2_df.cpp | 2 +- tests/testthat/test-sim-err-670.R | 3 +-- 3 files changed, 5 insertions(+), 9 deletions(-) diff --git a/NEWS.md b/NEWS.md index 6054b6726..545036d89 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,12 +4,9 @@ - Make sure that the object is a uncompressed rxode2 ui for solving with `rxSolve` (See #661) -- Ensure the residual simulation includes one simulation above the - number of observations; this allows trailing doses to access the - extra simulated residuals without accessing memory that has not been - allocated. This fixes the R crashed observed in #670. This did no - happen often and likely does not affect much code since the - `etTrans()` typically drops trailing doses. +- Fix #670 by using the last simulated observation residual when there + are trailing doses. + ## New features diff --git a/src/rxode2_df.cpp b/src/rxode2_df.cpp index 7c5959a61..07f50a71a 100644 --- a/src/rxode2_df.cpp +++ b/src/rxode2_df.cpp @@ -355,7 +355,7 @@ extern "C" SEXP rxode2_df(int doDose0, int doTBS) { if (updateErr){ for (j=0; j < errNcol; j++){ // The error pointer is updated if needed - par_ptr[svar[j]] = errs[errNrow*j+kk]; + par_ptr[svar[j]] = errs[errNrow*j+min2(kk, errNrow-1)]; } if ((doDose && evid!= 9) || (evid0 == 0 && isObs(evid)) || (evid0 == 1 && evid==0)){ // Only increment if this is an observation or of this a diff --git a/tests/testthat/test-sim-err-670.R b/tests/testthat/test-sim-err-670.R index fcf211e76..9599a3d1d 100644 --- a/tests/testthat/test-sim-err-670.R +++ b/tests/testthat/test-sim-err-670.R @@ -200,7 +200,6 @@ test_that("Simulate Error, issue #670", { row.names = c(NA, -218L), class = "data.frame") - crsh <- rxSolve(object = m, - events = d) + expect_error(rxSolve(object = m, events = d), NA) }) From e7698f05dae298b1838a44625bdd21f61d556f2f Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Sun, 7 Apr 2024 00:12:44 -0500 Subject: [PATCH 6/6] Different fix that doesn't change random numbers --- src/rxode2_df.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/rxode2_df.cpp b/src/rxode2_df.cpp index 07f50a71a..6c6633a31 100644 --- a/src/rxode2_df.cpp +++ b/src/rxode2_df.cpp @@ -355,12 +355,12 @@ extern "C" SEXP rxode2_df(int doDose0, int doTBS) { if (updateErr){ for (j=0; j < errNcol; j++){ // The error pointer is updated if needed - par_ptr[svar[j]] = errs[errNrow*j+min2(kk, errNrow-1)]; + par_ptr[svar[j]] = errs[errNrow*j+kk]; } 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){