decay <- sd_model(c(x = 100), params = list(k = 0.1),
                  flows = function(time, stocks, params) {
                    -params$k * stocks[["x"]]
                  })

test_that("model validation works", {
  expect_s3_class(decay, "sd_model")
  expect_error(sd_model(c(1), function(time, stocks, params) 0))
  expect_error(sd_model(c(x = 1), function(t, s, p) 0), "arguments")
  expect_error(sd_model(c(x = 1, y = 2),
                        function(time, stocks, params) 0), "2 flow")
  expect_output(print(decay), "Stocks")
})

test_that("Euler and RK4 approximate exponential decay", {
  e <- sd_simulate(decay, stop = 10, dt = 0.01)
  r <- sd_simulate(decay, stop = 10, dt = 0.5, method = "rk4")
  exact <- 100 * exp(-1)
  expect_equal(e$x[nrow(e)], exact, tolerance = 1e-3)
  expect_equal(r$x[nrow(r)], exact, tolerance = 1e-5)
  expect_equal(nrow(e), 11)
  expect_equal(e$time, 0:10)
})

test_that("simulation options are checked", {
  expect_error(sd_simulate(decay, stop = 0))
  expect_error(sd_simulate(decay, dt = 0.3, save_every = 1), "multiple")
  s <- sd_simulate(decay, stop = 4, params = list(k = 0))
  expect_true(all(s$x == 100))
  expect_s3_class(summary(s), "data.frame")
})

test_that("sd_update and helpers", {
  m <- sd_update(decay, params = list(k = 0.2), stocks = c(x = 5))
  expect_equal(m$params$k, 0.2)
  expect_equal(m$stocks[["x"]], 5)
  expect_error(sd_update(decay, stocks = c(z = 1)), "Unknown")
  f <- sd_lookup(c(0, 1), c(0, 10))
  expect_equal(f(c(-1, 0.5, 2)), c(0, 5, 10))
  expect_equal(sd_step(0:3, 2, 2), c(0, 0, 2, 2))
  expect_equal(sd_pulse(0:4, 6, 1, 2), c(0, 3, 3, 0, 0))
  expect_equal(sd_ramp(0:4, 1, 1, 3), c(0, 0, 1, 2, 2))
})