tests/testthat/test-theta-eta.R

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

rxTest({
  rigid <- rxode2({
    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 <- rxode2({
    a <- cos(theta)
  })

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

  test_that("theta/eta only parsing works", {
    expect_s3_class(theta, "rxode2")
    expect_s3_class(eta, "rxode2")
  })
})