tests/testthat/test-dfdy.R

In nlmixr2/rxode2: Facilities for Simulating from ODE-Based Models

rxTest({
  test_that("Specified jacobian is captured", {

    Vtpol2 <- rxode2("
d/dt(y)  = dy
d/dt(dy) = mu*(1-y^2)*dy - y
## Jacobian
df(y)/dy(dy)  = 1
df(dy)/dy(y)  = -2*dy*mu*y - 1
df(dy)/dy(dy) = mu*(1-y^2)
## Initial conditions
y(0) = 2
dy(0) = 0
## mu
mu = 1 ## nonstiff; 10 moderately stiff; 1000 stiff
")

    expect_true(any(rxDfdy(Vtpol2) == "df(y)/dy(dy)"))
    expect_true(any(rxDfdy(Vtpol2) == "df(dy)/dy(y)"))
    expect_true(any(rxDfdy(Vtpol2) == "df(dy)/dy(dy)"))
  })

  ## fixme multiple jacobian definitions raise an error.

  tmp <- rxode2("
d/dt(y)  = dy
d/dt(dy) = mu*(1-y^2)*dy - y
## Jacobian
df(y)/dy(dy)  = 1
df(dy)/dy(y)  = -2*dy*mu*y - 1
df(dy)/dy(dy) = mu*(1-y^2)
df(dy)/dy(dy) = mu*(1-y^2)
## Initial conditions
y(0) = 2
dy(0) = 0
## mu
mu = 1 ## nonstiff; 10 moderately stiff; 1000 stiff
")

  test_that("Doubled jacobain will compile correctly", {
    expect_s3_class(tmp, "rxode2")
  })

  test_that("Jacobian and sensitivity not specified.", {
    norm <- rxode2("
d/dt(y)  = dy
d/dt(dy) = mu*(1-y^2)*dy - y
## Initial conditions
y(0) = 2
dy(0) = 0
## mu
mu = 1 ## nonstiff; 10 moderately stiff; 1000 stiff
")

    expect_false(norm$calcJac)
    expect_false(norm$calcSens)
    rxDelete(norm)
  })

  test_that("Jacobian specified but sensitivity not specified.", {
    jac <- suppressMessages(rxode2("
d/dt(y)  = dy
d/dt(dy) = mu*(1-y^2)*dy - y
## Initial conditions
y(0) = 2
dy(0) = 0
## mu
mu = 1 ## nonstiff; 10 moderately stiff; 1000 stiff
", calcJac = TRUE))

    expect_true(jac$calcJac)
    expect_false(jac$calcSens)
    rxDelete(jac)
  })

  test_that("Sensitivity specified.", {
    sens <- suppressMessages(rxode2("
d/dt(y)  = dy
d/dt(dy) = mu*(1-y^2)*dy - y
## Initial conditions
y(0) = 2
dy(0) = 0
## mu
mu = 1 ## nonstiff; 10 moderately stiff; 1000 stiff
", calcSens = TRUE))

    expect_false(sens$calcJac)
    expect_true(sens$calcSens)
    rxDelete(sens)
  })

  test_that("Jac/Sens can be calculated from 'normal' model", {
    norm <- rxode2("
d/dt(y)  = dy
d/dt(dy) = mu*(1-y^2)*dy - y
## Initial conditions
y(0) = 2
dy(0) = 0
## mu
mu = 1 ## nonstiff; 10 moderately stiff; 1000 stiff
")
    jac <- norm
    jac <- suppressMessages(rxode2(jac, calcJac = TRUE))
    expect_true(jac$calcJac)
    expect_false(jac$calcSens)
    sens <- suppressMessages(rxode2(jac, calcSens = TRUE))
    expect_false(sens$calcJac)
    expect_true(sens$calcSens)
    full <- suppressMessages(rxode2(jac, calcSens = TRUE, calcJac = TRUE))
    expect_false(sens$calcJac)
    expect_true(sens$calcSens)
    rxDelete(jac)
    rxDelete(sens)
    rxDelete(norm)
    rxDelete(full)
  })

  test_that("Jacobian and sensitivity specified.", {
    sens <- suppressMessages(rxode2("
d/dt(y)  = dy
d/dt(dy) = mu*(1-y^2)*dy - y
## Initial conditions
y(0) = 2
dy(0) = 0
## mu
mu = 1 ## nonstiff; 10 moderately stiff; 1000 stiff
", calcSens = TRUE))

    norm <- suppressMessages(rxode2(sens, calcSens = FALSE))
    expect_false(norm$calcJac)
    expect_false(norm$calcSens)

    jac <- suppressMessages(rxode2(sens, calcJac = TRUE))
    expect_true(jac$calcJac)
    expect_false(jac$calcSens)
    expect_false(sens$calcJac)
    expect_true(sens$calcSens)
  })

  test_that("Conditional Sensitivities", {

    transit.if <- rxode2({
      ## Table 3 from Savic 2007
      cl <- 17.2 # (L/hr)
      vc <- 45.1 # L
      ka1 <- 0.38 # 1/hr
      ka2 <- 0.2 # 1/hr
      mtt <- 0.37 # hr
      bio <- 1
      n <- 20.1
      k <- cl / vc
      if (pop1 == 1) {
        ka <- ka1
      } else {
        ka <- ka2
      }
      d/dt(depot) <- transit(n, mtt, bio) - ka * depot
      d/dt(cen) <- ka * depot - k * cen
    })

    expect_false(transit.if$calcJac)
    expect_false(transit.if$calcSens)

    jac <- suppressMessages(rxode2(transit.if, calcJac = TRUE))
    expect_true(jac$calcJac)
    expect_false(jac$calcSens)

    sens <- suppressMessages(rxode2(transit.if, calcSens = TRUE))
    expect_false(sens$calcJac)
    expect_true(sens$calcSens)

    full <- suppressMessages(rxode2(transit.if, calcSens = TRUE, calcJac = TRUE))
    expect_true(full$calcJac)
    expect_true(full$calcSens)

  })

  test_that("Transit Sensitivities", {

    mod <- rxode2("
## Table 3 from Savic 2007
cl = 17.2 # (L/hr)
vc = 45.1 # L
ka = 0.38 # 1/hr
mtt = 0.37 # hr
bio=1
n = 20.1
k = cl/vc
ktr = (n+1)/mtt
## note that lgammafn is the same as lgamma in R.
d/dt(depot) = exp(log(bio*podo(depot))+log(ktr)+n*log(ktr*tad(depot))-ktr*tad(depot)-lgammafn(n+1))-ka*depot
d/dt(cen) = ka*depot-k*cen
")

    mod <- suppressMessages(rxode2(mod, calcSens = TRUE))

    et <- eventTable()
    et$add.sampling(seq(0, 10, length.out = 200))
    et$add.dosing(20, start.time = 0, evid=7)

    transit <- suppressWarnings({
      rxSolve(mod, et)
    })

    ## Used the log(0) protection since the depot_mtt sensitivity
    ## includes log(t*...) and t=0
    ##
    ## These results are now system dependent since it uses
    ## log(sqrt(DOUBLE_EPS)) and DOUBLE_EPS varies by platform, so
    ## just make sure the results are not NA.
    expect_true(all(!is.na(transit[["_sens_depot_mtt"]])))
    expect_equal(transit$depot.mtt, transit[["_sens_depot_mtt"]])
    expect_equal(transit$depot_mtt, transit[["_sens_depot_mtt"]])
    expect_equal(transit[["depot_mtt"]], transit[["_sens_depot_mtt"]])
    expect_equal(transit[["depot.mtt"]], transit[["_sens_depot_mtt"]])

    mod <- rxode2({
      ## Table 3 from Savic 2007
      cl <- 17.2 # (L/hr)
      vc <- 45.1 # L
      tvka <- 0.38 # 1/hr
      tvmtt <- 0.37 # hr
      ka <- tvka * exp(eta_ka)
      mtt <- tvmtt * exp(eta_mtt)
      bio <- 1
      n <- 20.1
      k <- cl / vc
      ktr <- (n + 1) / mtt
      ## note that lgammafn is the same as lgamma in R.
      d/dt(depot) <- exp(log(bio * podo(depot)) + log(ktr) + n * log(ktr * tad(depot)) - ktr * tad(depot) -
                           lgammafn(n + 1)) - ka * depot
      d/dt(cen) <- ka * depot - k * cen
    })

    transit <- suppressMessages(rxode2(mod, calcSens = c("eta_ka", "eta_mtt")))
    expect_true(all(!is.na(transit[["_sens_depot_mtt"]])))

    ## tmp <- rxode2(mod, calcSens=list(eta=c("eta_ka", "eta_mtt"), theta=c("cl", "vc")));
  })
})