From efa096769e7e8e7c777a46f61c7a496580ce551f Mon Sep 17 00:00:00 2001 From: Bill Denney Date: Mon, 17 Oct 2022 13:42:00 -0400 Subject: [PATCH 1/4] Make table integration more robust --- R/resid.R | 89 ++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 58 insertions(+), 31 deletions(-) diff --git a/R/resid.R b/R/resid.R index 73ebfccf..50830334 100644 --- a/R/resid.R +++ b/R/resid.R @@ -84,45 +84,72 @@ nmObjGet.foceiThetaEtaParameters <- function(x, ...) { stop("cannot solve with `model` NULL", call.=FALSE) } keep <- unique(c(keep, "nlmixrRowNums")) - .res <- .foceiSolveWithId(model, pars, fit$dataSav, - returnType = returnType, - atol = fit$atol[1], rtol = fit$rtol[1], - maxsteps = fit$maxstepsOde, - hmin = fit$hmin, hmax = fit$hmax, hini = fit$hini, - maxordn = fit$maxordn, - maxords = fit$maxords, method = fit$methodOde, - keep=keep, addDosing=addDosing, subsetNonmem=subsetNonmem, addCov=addCov) - rxode2::rxSolveFree() - if (any(is.na(.res$rx_pred_)) && fit$methodOde == 2L) { - .res <- .foceiSolveWithId(model, pars, fit$dataSav, - returnType = returnType, - atol = fit$atol[1], rtol = fit$rtol[1], - maxsteps = fit$maxstepsOde * 2, - hmin = fit$hmin, hmax = fit$hmax / 2, hini = fit$hini, - maxordn = fit$maxordn, - maxords = fit$maxords, method = "lsoda", - keep=keep, addDosing=addDosing, subsetNonmem=subsetNonmem, addCov=addCov) - rxode2::rxSolveFree() - if (any(is.na(.res$rx_pred_))) { + # The numeric versions are at + # https://github.com/nlmixr2/rxode2/blob/7e27a7842ca0b5dd849ea75833bc7c34be729e31/R/rxsolve.R#L804, + # but keeping them in sync will be fragile. Only using the character + # versions. + currentOdeMethod <- fit$methodOde + allOdeMethods <- + setdiff( + eval(formals(rxode2::rxSolve)$method), + # ignore indLin for now + "indLin" + ) + # Fallback methods based on discussion in + # https://github.com/nlmixr2/nlmixr2est/issues/254 + if (currentOdeMethod %in% "dop853") { + allOdeMethods <- "liblsoda" + } else if (currentOdeMethod %in% c("liblsoda", "lsoda")) { + allOdeMethods <- "dop853" + } # otherwise, use all the methods + odeMethods <- + append( + list(currentOdeMethod), + as.list(setdiff(allOdeMethods, currentOdeMethod)) + ) + startHmax <- fit$hmax + startMaxStepsOde <- fit$maxstepsOde + failedMethods <- character() + isFirstFit <- TRUE + recalc <- TRUE + while (recalc & length(odeMethods) > 0) { + # Iterate through ODE methods + recalcN <- 0 + currentOdeMethod <- odeMethods[[1]] + odeMethods <- odeMethods[-1] + currentHmax <- fit$hmax + currentMaxStepsOde <- fit$maxStepsOde + while (recalc & recalcN < fit$control$stickyRecalcN) { + # Iterate down Hmax and up Max Steps .res <- .foceiSolveWithId(model, pars, fit$dataSav, returnType = returnType, atol = fit$atol[1], rtol = fit$rtol[1], - maxsteps = fit$maxstepsOde * 2, - hmin = fit$hmin, hmax = fit$hmax / 2, hini = fit$hini, + maxsteps = currentMaxStepsOde, + hmin = fit$hmin, hmax = currentHmax, hini = fit$hini, maxordn = fit$maxordn, - maxords = fit$maxords, method = "dop853", + maxords = fit$maxords, method = currentOdeMethod, keep=keep, addDosing=addDosing, subsetNonmem=subsetNonmem, addCov=addCov) rxode2::rxSolveFree() - if (any(is.na(.res$rx_pred_))) { - warning("Problems solving ", what, " liblsoda, lsoda and dop853") - } else { - warning("Problems solving ", what, " liblsoda and lsoda switched to dop853") - } - } else { - warning("Problems solving ", what, " liblsoda switched to lsoda") + recalcN <- recalcN + 1 + currentHmax <- currentHmax / 2 + currentMaxStepsOde <- currentMaxStepsOde * 2 + recalc <- any(is.na(.res$rx_pred_)) + } + if (recalc) { + failedMethods <- c(failedMethods, currentOdeMethod) } + if (isFirstFit) { + isFirstFit <- FALSE + .resFirst <- .res + } + } + if (recalc) { + .res <- .resFirst + warning("Problems solving ", what, " with ", paste(failedMethods, collapse = ", "), ", returning results from the first method") + } else if (length(failedMethods) > 0) { + warning("Problems solving ", what, " with ", paste(failedMethods, collapse = ", "), ", returning results from ", currentOdeMethod) } - return(.res) + .res } #' Create a ipred/pred list from the focei style model From ed45f00cc0296c990935e44527c291fcf052518e Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Tue, 18 Oct 2022 18:09:56 -0500 Subject: [PATCH 2/4] Add test and redo solving --- R/resid.R | 34 +++++++++++++++++++++++----------- tests/testthat/test-restart.R | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 56 insertions(+), 11 deletions(-) create mode 100644 tests/testthat/test-restart.R diff --git a/R/resid.R b/R/resid.R index 50830334..0c0e7add 100644 --- a/R/resid.R +++ b/R/resid.R @@ -89,6 +89,13 @@ nmObjGet.foceiThetaEtaParameters <- function(x, ...) { # but keeping them in sync will be fragile. Only using the character # versions. currentOdeMethod <- fit$methodOde + if (!inherits(currentOdeMethod, "character")) { + cur <- as.integer(currentOdeMethod)+1L + attr(cur, "levels") <- c("dop853", "lsoda", "liblsoda", "indLin") + attr(cur, "class") <- "factor" + currentOdeMethod <- as.character(cur) + print(currentOdeMethod) + } allOdeMethods <- setdiff( eval(formals(rxode2::rxSolve)$method), @@ -107,33 +114,38 @@ nmObjGet.foceiThetaEtaParameters <- function(x, ...) { list(currentOdeMethod), as.list(setdiff(allOdeMethods, currentOdeMethod)) ) - startHmax <- fit$hmax - startMaxStepsOde <- fit$maxstepsOde failedMethods <- character() isFirstFit <- TRUE recalc <- TRUE + maxAtolRtol <- fit$foceiControl$rxControl$maxAtolRtolFactor + recalcFactor <- fit$foceiControl$odeRecalcFactor while (recalc & length(odeMethods) > 0) { # Iterate through ODE methods recalcN <- 0 currentOdeMethod <- odeMethods[[1]] odeMethods <- odeMethods[-1] - currentHmax <- fit$hmax - currentMaxStepsOde <- fit$maxStepsOde - while (recalc & recalcN < fit$control$stickyRecalcN) { - # Iterate down Hmax and up Max Steps + .atol <- fit$atol[1] + .rtol <- fit$rtol[1] + while (recalc & recalcN < fit$foceiControl$stickyRecalcN) { + # Iterate up atol/rtol .res <- .foceiSolveWithId(model, pars, fit$dataSav, returnType = returnType, - atol = fit$atol[1], rtol = fit$rtol[1], - maxsteps = currentMaxStepsOde, - hmin = fit$hmin, hmax = currentHmax, hini = fit$hini, + atol = .atol, rtol = .rtol, + maxsteps = fit$maxstepsOde*ifelse(isFirstFit, 1, 2), + hmin = fit$hmin, hmax = fit$hmax/*ifelse(isFirstFit, 1, 2), hini = fit$hini, maxordn = fit$maxordn, maxords = fit$maxords, method = currentOdeMethod, keep=keep, addDosing=addDosing, subsetNonmem=subsetNonmem, addCov=addCov) rxode2::rxSolveFree() recalcN <- recalcN + 1 - currentHmax <- currentHmax / 2 - currentMaxStepsOde <- currentMaxStepsOde * 2 recalc <- any(is.na(.res$rx_pred_)) + if (recalc) { + .atol <- min(.atol*recalcFactor, maxAtolRtol) + .rtol <- min(.rtol*recalcFactor, maxAtolRtol) + if (.atol == maxAtolRtol && .rtol == maxAtolRtol) { + recalcN <- fit$foceiControl$stickyRecalcN + 1 + } + } } if (recalc) { failedMethods <- c(failedMethods, currentOdeMethod) diff --git a/tests/testthat/test-restart.R b/tests/testthat/test-restart.R new file mode 100644 index 00000000..5d90ae54 --- /dev/null +++ b/tests/testthat/test-restart.R @@ -0,0 +1,33 @@ +test_that("restart resid", { + + lobo <- function() { + ini({ + lkng <- log(0.02) + ltau <- log(c(1, 34.1, 500)) + lec50 <- log(c(1, 50, 2000)) + kmax <- 0.01 + propSd <- c(0, 0.3) + addSd <- c(0, 50) + }) + model({ + kng <- exp(lkng) + tau <- exp(ltau) + taulast <- tau + ec50 <- exp(lec50) + edrug <- kmax * cp/(ec50 + cp) + tumor(0) <- tumor0 + d/dt(transit1) <- (edrug - transit1)/tau + d/dt(transit2) <- (transit1 - transit2)/tau + d/dt(transit3) <- (transit2 - transit3)/tau + d/dt(transitlast) <- transit3/tau - transitlast/taulast + d/dt(tumor) <- kng * tumor - transitlast * tumor + tumor ~ prop(propSd) + add(addSd) + }) + } + + prepfit <- qs::qdeserialize(qs::base91_decode("un]\"BAAA@QRtHACAAAAAAAuWeBAABdk1kus^^d8Ah9}?=Z:alBMc4Iv(F\":C?hVBAMZRxwFfBB7IB.y6FTL)yFQA@v;KlgSH3Vn~rL/,{CP/ez~`.3>$Wj$rcy==//#}Pu?\"V(RgKtf1J3qQg/yI7*1]/WV=iegVmPs3?a\":kEMu~/*zmX#;E4`i@It`]Ouu[N]T8G3!4A.^j0M%xmxwMU%ET6#~xvX9rH!;S53gbLTWY]Fcri\"]7\"|Z^W{xobiiTc~DLN_;.Itj(INGKCupDYxEA^!GzfHO=aDW&(I)z}0*mZD\"^b.O!QdY2rVRD;~Z*HB(]G_Dpj]*0A`]+7VoGDF@,vk>jx}tFI>MVOnZojuABN9Bt\"O~V[n6U[kn|W74&xR7CL(Skn:CA)NP`||hQ%w/i+&c8$#KxsFdb4,qI\"Fl&lLg,?$eh&s{`QxtwPWi$GX<[*<0{to@[:NAy}a=O`wedEA*Abqhz2bL2sfII3ZRJR#5q~:FeBW%/F<]`(?Q:c(qc,DZ_d.&|J(NW~Q4kz;Us(7e+Z0YGMdvf.%XRgD]FA2D10sl^KxuPXvXSm+p}ndVY!3`o}Iq+M;i~mLmr1In0~ymm]K2x9g9Ij.UkBOTriq+93#eKT9uaMmvF7L^aDwL?sj>s|}[XMQdEx(yS|vIDfwqH:YXc(EfrzuplGn3`|X=ObNnD%;3(ST3tWr^D+vDG=cjKk!^:5ZfpXK5/dqF@dW6+*lb~@\"*H_t@3rRG9w|kG1!SRghm{sUOcDQ?.gYd?c;:IC~RF6lEfd4X;I3+4Y&8.x0P[h8]Ffo)$Y=7)|08Vid5@2btxd4I[0>E7~+CCX~|{Ve0qf^aQ1M4vsQN[uCto@(_b&&NA0eV*R}~Tp#Kz+AL~%_Wu@J&>C__{[#o3Lmxh$%R;264SAA)O;,:*:p.s:z+=pot[NkrFE@[}VoBF=Pv:;Cj#YR]M@a{68?;s&4Y4Zd[D`znL|kp1n)}0+Q`(/LcQM42JIA}ZFH+9>#dQLs>JX?#GT_M&#wui|N/HWs/WyJ/CSn*BaX!:]**?VL`zPu?xA|.+kn)M[xCjLHQ@8A6HX^T4hXT7F~JQw])fR/^Y>YBD([_QsK;[iYPcbW2Z\"BA9q~@t^rOgPg95OnA01xfA00E8k)eQtyH#MzpQ$4|et6W(eGPRmx0B?bGrM@DA31ha^UZWMe52Egm6nn3=lcs3[gq!vBvrKiTSf],3$wjk7^`BPE#{kPh.CcbM*h:9+XPRfGx1P^)X!utw{#;**|LZtNdAwx\"IC?m[]Wil9!m3Tx}{Xx@9;`=+ODLPp)w]zc\"D.rl7/d:ismtgiesu]/CF[9v7ICZ:U?wky`@wQn])q@|Dd#IfQbwAQF`7]{f;hJCK=Vv)N6M4x4>G.u%o7I5fCay=\"^,?M!600bO\"MBTL9~4}eDkY:YkGc6G9oMzo3<9[Ye`GFvPyXiw5$Z]*PByzRza9ik5LJ;}@p;JH}+.,1om5gxAkHdn%`#kbKMJuiDOa#*M!h.2dUje501>,MIW$207SwS$M~FCMaG:`&l?jK]e\"!>HG]Tw0Xbm7SOIjt.B%cH5Ewy_YgfI6(#$I]1Q4Jte1:^W%XdqR4C#OV#._I`G$k|#>NfxcZp!Y]YnX~yZ0nx),%jfE\"u7yS,:y&oJV`;4nMg0g[a]Tidjt|y[*~_8H;*X61T/1kv9x$f/3N?~>xL$mO^0Nk6W|U[wD")) + + + sim <- nlmixr2(lobo, data=prepfit, "focei") + +}) From bb41e9d71a4dfb9a77ea0d8ca685ded7ef9bac7b Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Tue, 18 Oct 2022 18:48:34 -0500 Subject: [PATCH 3/4] More redo solving settings --- R/resid.R | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/R/resid.R b/R/resid.R index 0c0e7add..d45d4cdc 100644 --- a/R/resid.R +++ b/R/resid.R @@ -94,7 +94,6 @@ nmObjGet.foceiThetaEtaParameters <- function(x, ...) { attr(cur, "levels") <- c("dop853", "lsoda", "liblsoda", "indLin") attr(cur, "class") <- "factor" currentOdeMethod <- as.character(cur) - print(currentOdeMethod) } allOdeMethods <- setdiff( @@ -126,19 +125,21 @@ nmObjGet.foceiThetaEtaParameters <- function(x, ...) { odeMethods <- odeMethods[-1] .atol <- fit$atol[1] .rtol <- fit$rtol[1] + ## message(currentOdeMethod) while (recalc & recalcN < fit$foceiControl$stickyRecalcN) { # Iterate up atol/rtol + ## message("\t", .atol, " ", .rtol) .res <- .foceiSolveWithId(model, pars, fit$dataSav, returnType = returnType, atol = .atol, rtol = .rtol, - maxsteps = fit$maxstepsOde*ifelse(isFirstFit, 1, 2), - hmin = fit$hmin, hmax = fit$hmax/*ifelse(isFirstFit, 1, 2), hini = fit$hini, - maxordn = fit$maxordn, - maxords = fit$maxords, method = currentOdeMethod, + maxsteps = fit$maxstepsOde, + hmin = fit$hmin, hmax = fit$hmax, hini = fit$hini, + maxordn = fit$maxordn, maxords = fit$maxords, + method = rxode2::odeMethodToInt(currentOdeMethod), keep=keep, addDosing=addDosing, subsetNonmem=subsetNonmem, addCov=addCov) rxode2::rxSolveFree() - recalcN <- recalcN + 1 recalc <- any(is.na(.res$rx_pred_)) + recalcN <- recalcN + 1 if (recalc) { .atol <- min(.atol*recalcFactor, maxAtolRtol) .rtol <- min(.rtol*recalcFactor, maxAtolRtol) From 4fd55657c847a324f0d94965e7dcbacd64002789 Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Tue, 18 Oct 2022 19:16:41 -0500 Subject: [PATCH 4/4] Test restart works with the non submitted CRAN version of rxode2 --- tests/testthat/test-restart.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-restart.R b/tests/testthat/test-restart.R index 5d90ae54..250c3fa3 100644 --- a/tests/testthat/test-restart.R +++ b/tests/testthat/test-restart.R @@ -28,6 +28,6 @@ test_that("restart resid", { prepfit <- qs::qdeserialize(qs::base91_decode("un]\"BAAA@QRtHACAAAAAAAuWeBAABdk1kus^^d8Ah9}?=Z:alBMc4Iv(F\":C?hVBAMZRxwFfBB7IB.y6FTL)yFQA@v;KlgSH3Vn~rL/,{CP/ez~`.3>$Wj$rcy==//#}Pu?\"V(RgKtf1J3qQg/yI7*1]/WV=iegVmPs3?a\":kEMu~/*zmX#;E4`i@It`]Ouu[N]T8G3!4A.^j0M%xmxwMU%ET6#~xvX9rH!;S53gbLTWY]Fcri\"]7\"|Z^W{xobiiTc~DLN_;.Itj(INGKCupDYxEA^!GzfHO=aDW&(I)z}0*mZD\"^b.O!QdY2rVRD;~Z*HB(]G_Dpj]*0A`]+7VoGDF@,vk>jx}tFI>MVOnZojuABN9Bt\"O~V[n6U[kn|W74&xR7CL(Skn:CA)NP`||hQ%w/i+&c8$#KxsFdb4,qI\"Fl&lLg,?$eh&s{`QxtwPWi$GX<[*<0{to@[:NAy}a=O`wedEA*Abqhz2bL2sfII3ZRJR#5q~:FeBW%/F<]`(?Q:c(qc,DZ_d.&|J(NW~Q4kz;Us(7e+Z0YGMdvf.%XRgD]FA2D10sl^KxuPXvXSm+p}ndVY!3`o}Iq+M;i~mLmr1In0~ymm]K2x9g9Ij.UkBOTriq+93#eKT9uaMmvF7L^aDwL?sj>s|}[XMQdEx(yS|vIDfwqH:YXc(EfrzuplGn3`|X=ObNnD%;3(ST3tWr^D+vDG=cjKk!^:5ZfpXK5/dqF@dW6+*lb~@\"*H_t@3rRG9w|kG1!SRghm{sUOcDQ?.gYd?c;:IC~RF6lEfd4X;I3+4Y&8.x0P[h8]Ffo)$Y=7)|08Vid5@2btxd4I[0>E7~+CCX~|{Ve0qf^aQ1M4vsQN[uCto@(_b&&NA0eV*R}~Tp#Kz+AL~%_Wu@J&>C__{[#o3Lmxh$%R;264SAA)O;,:*:p.s:z+=pot[NkrFE@[}VoBF=Pv:;Cj#YR]M@a{68?;s&4Y4Zd[D`znL|kp1n)}0+Q`(/LcQM42JIA}ZFH+9>#dQLs>JX?#GT_M&#wui|N/HWs/WyJ/CSn*BaX!:]**?VL`zPu?xA|.+kn)M[xCjLHQ@8A6HX^T4hXT7F~JQw])fR/^Y>YBD([_QsK;[iYPcbW2Z\"BA9q~@t^rOgPg95OnA01xfA00E8k)eQtyH#MzpQ$4|et6W(eGPRmx0B?bGrM@DA31ha^UZWMe52Egm6nn3=lcs3[gq!vBvrKiTSf],3$wjk7^`BPE#{kPh.CcbM*h:9+XPRfGx1P^)X!utw{#;**|LZtNdAwx\"IC?m[]Wil9!m3Tx}{Xx@9;`=+ODLPp)w]zc\"D.rl7/d:ismtgiesu]/CF[9v7ICZ:U?wky`@wQn])q@|Dd#IfQbwAQF`7]{f;hJCK=Vv)N6M4x4>G.u%o7I5fCay=\"^,?M!600bO\"MBTL9~4}eDkY:YkGc6G9oMzo3<9[Ye`GFvPyXiw5$Z]*PByzRza9ik5LJ;}@p;JH}+.,1om5gxAkHdn%`#kbKMJuiDOa#*M!h.2dUje501>,MIW$207SwS$M~FCMaG:`&l?jK]e\"!>HG]Tw0Xbm7SOIjt.B%cH5Ewy_YgfI6(#$I]1Q4Jte1:^W%XdqR4C#OV#._I`G$k|#>NfxcZp!Y]YnX~yZ0nx),%jfE\"u7yS,:y&oJV`;4nMg0g[a]Tidjt|y[*~_8H;*X61T/1kv9x$f/3N?~>xL$mO^0Nk6W|U[wD")) - sim <- nlmixr2(lobo, data=prepfit, "focei") + sim <- nlmixr2(lobo, data=prepfit, "focei", control=foceiControl(print=0)) })