tests/testthat/test-theta-eta.R

In nlmixrdevelopment/RxODE: Facilities for Simulating from ODE-Based Models

rxodeTest(
  {
    context("Test RxODE THETA/ETA support")
    library(digest)

    rigid <- RxODE({
      y1(0) <- 1
      y2(0) <- 0
      y3(0) <- 0.9
      a1 <- theta[1]
      a2 <- theta[2]
      a3 <- theta[3] + eta[1]
      d / dt(y1) <- a1 * y2 * y3
      d / dt(y2) <- a2 * y1 * y3
      d / dt(y3) <- a3 * y1 * y2
    })

    et <- eventTable()
    et$add.sampling(seq(0, 20, by = 0.01))

    out <- solve(rigid, et, theta = c(-2, 1.25, -0.5), eta = 0)

    test_that("Test rigid body example", {
      expect_equal(
        round(as.data.frame(out[c(1:15, seq(2001 - 15, 2001)), ]), 3),
        structure(list(time = c(0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 19.85, 19.86, 19.87, 19.88, 19.89, 19.9, 19.91, 19.92, 19.93, 19.94, 19.95, 19.96, 19.97, 19.98, 19.99, 20), a1 = c(-2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2), a2 = c(1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25), a3 = c(-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5), y1 = c(1, 1, 1, 0.999, 0.998, 0.997, 0.996, 0.995, 0.994, 0.992, 0.99, 0.988, 0.985, 0.983, 0.98, 0.749, 0.74, 0.731, 0.722, 0.713, 0.704, 0.694, 0.685, 0.675, 0.666, 0.656, 0.646, 0.636, 0.626, 0.616, 0.606), y2 = c(0, 0.011, 0.022, 0.034, 0.045, 0.056, 0.067, 0.079, 0.09, 0.101, 0.112, 0.123, 0.134, 0.145, 0.156, 0.524, 0.532, 0.539, 0.547, 0.554, 0.562, 0.569, 0.576, 0.583, 0.59, 0.597, 0.603, 0.61, 0.616, 0.623, 0.629), y3 = c(0.9, 0.9, 0.9, 0.9, 0.9, 0.899, 0.899, 0.899, 0.898, 0.898, 0.897, 0.897, 0.896, 0.895, 0.895, 0.837, 0.835, 0.833, 0.831, 0.829, 0.827, 0.825, 0.823, 0.821, 0.819, 0.817, 0.815, 0.813, 0.811, 0.809, 0.807)), row.names = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 1986L, 1987L, 1988L, 1989L, 1990L, 1991L, 1992L, 1993L, 1994L, 1995L, 1996L, 1997L, 1998L, 1999L, 2000L, 2001L), class = "data.frame")
      )
    })

    out <- solve(rigid, et, theta = c(-2, 1.25, -0.5), eta = 1)

    test_that("Test rigid body example", {
      expect_equal(
        round(as.data.frame(out[c(1:15, seq(2001 - 15, 2001)), ]), 3),
        structure(list(time = c(
          0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06,
          0.07, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 19.85, 19.86,
          19.87, 19.88, 19.89, 19.9, 19.91, 19.92, 19.93, 19.94, 19.95,
          19.96, 19.97, 19.98, 19.99, 20
        ), a1 = c(
          -2, -2, -2, -2, -2, -2,
          -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2, -2,
          -2, -2, -2, -2, -2, -2, -2, -2, -2
        ), a2 = c(
          1.25, 1.25, 1.25,
          1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25,
          1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25, 1.25,
          1.25, 1.25, 1.25, 1.25, 1.25, 1.25
        ), a3 = c(
          0.5, 0.5, 0.5, 0.5,
          0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
          0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
          0.5
        ), y1 = c(
          1, 1, 1, 0.999, 0.998, 0.997, 0.996, 0.995, 0.994,
          0.992, 0.99, 0.988, 0.985, 0.983, 0.98, 0.413, 0.427, 0.441,
          0.456, 0.47, 0.484, 0.497, 0.511, 0.525, 0.538, 0.551, 0.564,
          0.577, 0.59, 0.602, 0.614
        ), y2 = c(
          0, 0.011, 0.022, 0.034, 0.045,
          0.056, 0.067, 0.079, 0.09, 0.101, 0.112, 0.123, 0.135, 0.146,
          0.157, -0.72, -0.715, -0.709, -0.704, -0.698, -0.692, -0.686,
          -0.68, -0.673, -0.666, -0.66, -0.653, -0.646, -0.639, -0.631,
          -0.624
        ), y3 = c(
          0.9, 0.9, 0.9, 0.9, 0.9, 0.901, 0.901, 0.901,
          0.902, 0.902, 0.903, 0.903, 0.904, 0.905, 0.905, 1.009, 1.007,
          1.006, 1.004, 1.002, 1.001, 0.999, 0.997, 0.996, 0.994, 0.992,
          0.99, 0.988, 0.986, 0.985, 0.983
        )), row.names = c(
          1L, 2L, 3L,
          4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 1986L,
          1987L, 1988L, 1989L, 1990L, 1991L, 1992L, 1993L, 1994L, 1995L,
          1996L, 1997L, 1998L, 1999L, 2000L, 2001L
        ), class = "data.frame")
      )
    })

    theta <- RxODE({
      a <- cos(theta)
    })

    eta <- RxODE({
      a <- cos(eta)
    })

    test_that("theta/eta only parsing works", {
      expect_equal(class(theta), "RxODE")
      expect_equal(class(eta), "RxODE")
    })
  },
  silent = TRUE,
  test = "lvl2"
)