# tests/testthat/test_generator.R In CharlotteJana/pdmpsim: Simulate Piecewise Deterministic Markov Processes

```context("generator")

test_that("method generator is correct for a pdmp with one discrete variable", {
skip_if_not_installed("Deriv")

simplePdmp <- pdmpModel(
descr = "a simple pdmp",
init = c(f = 0, d = 0),
times = c(from = 0, to = 10, by = 0.01),
discStates = list(d = c(-1, 0, 1)),
dynfunc = function(t, x, parms) c(x["d"], 0),
ratefunc = function(t, x, parms) c(1+x["d"], 1-x["d"]),
jumpfunc = function(t, x, parms, jtype){
c(0, switch(jtype, x["d"]-1, x["d"]+1))
}
)

fvals <- seq(from = 0, to = 3, by = 0.5)

g <- function(d, f) d*f
Qg <- function(d, fvals)
sapply(fvals, function(f)
generator(simplePdmp)(g)(t = 10, x = c("d" = d, "f" = f)))
test1 <- cbind(Qg(-1, fvals),
Qg(0, fvals),
Qg(1, fvals))
test2 <- cbind(1+2*fvals,
rep(0, length(fvals)),
1-2*fvals)
expect_equal(test1, test2)

h <- function(d, f) exp(f) + d
Qh <- function(d, fvals)
sapply(fvals, function(f)
generator(simplePdmp)(h)(t = 10, x = c("d" = d, "f" = f)))
test1 <- cbind(Qh(-1, fvals),
Qh(0, fvals),
Qh(1, fvals))
test2 <- cbind(-3*exp(fvals)+4,
-2*exp(fvals)+2 ,
-exp(fvals))
expect_equal(test1, test2)
})

test_that("method generator is correct for a pdmp with 2 discrete variables", {
skip_if_not_installed("Deriv")

discPdmp <- pdmpModel(
descr = "a pdmp with 2 discrete variables and 2 jumptypes",
init = c(f = 0, d1 = 0, d2 = 0),
discStates = list(d1 = 0:1, d2 = 0:1),
times = c(from = 0, to = 10, by = 0.01),
dynfunc = function(t, x, parms) c(x["d1"]-x["d2"], 0, 0),
ratefunc = function(t, x, parms) c(1,1),
jumpfunc = function(t, x, parms, jtype){
c(0, switch(jtype, c((x["d1"]-1)^2, x["d2"]),
c(x["d1"], (x["d2"]-1)^2)))
}
)

# The simulations of this model are equivalent ot the simulations of the
# model used in the test before. Therefore the generators should be the same.

fvals <- seq(from = 0, to = 3, by = 0.5)

g <- function(d1, d2, f) (d1-d2)*f
Qg <- function(d1, d2, fvals)
sapply(fvals, function(f)
generator(discPdmp)(g)(t = 10, x = c("d1" = d1, "d2" = d2, "f" = f)))
test1 <- cbind(Qg(0, 1, fvals),
Qg(0, 0, fvals),
Qg(1, 1, fvals),
Qg(1, 0, fvals))
test2 <- cbind(1+2*fvals,
rep(0, length(fvals)),
rep(0, length(fvals)),
1-2*fvals)
expect_equal(test1, test2)

h <- function(d1, d2, f) exp(f) + d1 - d2
Qh <- function(d1, d2, fvals)
sapply(fvals, function(f)
generator(discPdmp)(h)(t = 10, x = c("d1" = d1, "d2" = d2, "f" = f)))
test1 <- cbind(Qh(0, 1, fvals),
Qh(0, 0, fvals),
Qh(1, 1, fvals),
Qh(1, 0, fvals))
test2 <- cbind(-3*exp(fvals)+4,
-2*exp(fvals)+2 ,
-2*exp(fvals)+2,
-exp(fvals))
expect_equal(test1, test2)
})

test_that("method generator is correct for the toggle switch model", {
skip_if_not_installed("Deriv")

toggleSwitch <- pdmpModel(
descr = "toggle switch with two promotors",
parms = list(bA = 0.5, bB = 1,   aA = 2, aB = 4,
k01A = 0.5, k10A = 2, k01B = 1/3, k10B = 3),
init = c(fA = 0.5, fB = 0.5, dA = 1, dB = 1),
discStates = list(dA = 0:1, dB = 0:1),
dynfunc = function(t, x, parms) {
df <- with(as.list(c(x, parms)),
c(-bA*fA + aA*dA, -bB*fB + aB*dB))
return(c(df, 0, 0))
},
ratefunc = function(t, x, parms) {
return(with(as.list(c(x, parms)),
c(switch(dB+1, k01A, k10A*fB),
switch(dA+1, k01B, k10B*fA))))
},
jumpfunc = function(t, x, parms, jtype) {
return(with(as.list(c(x, parms)),
c(fA, fB, switch(jtype, c(dA, 1-dB), c(1-dA, dB)))))
},
times = c(from = 0, to = 10, by = 0.01)
)

fvals <- expand.grid(1:10, -1:3)
fAs <- fvals[,1]
fBs <- fvals[,2]

g <- function(dA, dB, fA, fB) fA * fB + dA # the order of the arguments
# is different to the order in
# init(toggleSwitch)
Qg <- function(dA, dB, fvals) {
sapply(seq_len(nrow(fvals)), function(i)
generator(toggleSwitch)(g)(t = 10, x = c("dA" = dA, "dB" = dB,
"fA" = fvals[i,1],
"fB" = fvals[i,2])))
}
test1 <- cbind(Qg(0, 0, fvals),
Qg(0, 1, fvals),
Qg(1, 0, fvals),
Qg(1, 1, fvals))
test2 <- with(as.list(parms(toggleSwitch)), -bA*fAs*fBs - bB*fBs*fAs +
cbind(k01B, k01B + aB*fAs, aA*fBs -k10B*fAs,
aA*fBs - k10B*fAs + aB*fAs))
dimnames(test2) <- NULL
expect_equal(test1, test2)
})
```
CharlotteJana/pdmpsim documentation built on July 2, 2019, 5:37 a.m.