## Analysis of dynamic HCMV infection
source("bandleFUN.R")
require(pRoloc)
require(pRolocdata)


data("beltran2016HCMV48")
data("beltran2016MOCK48")
load("data/beltran48hypers.rda")

idx <- intersect(rownames(beltran2016MOCK48), rownames(beltran2016HCMV48))

HCMV48_sub <- beltran2016HCMV48[idx, ]
MOCK48_sub <- beltran2016MOCK48[idx, ]


beltrandyn48 <- c(MOCK48_sub, HCMV48_sub) 
beltrandyn48 <- lapply(beltrandyn48, function(x) normalise(x, method = "center.median"))

K <- length(getMarkerClasses(MOCK48_sub))

dir_prior = diag(rep(1, K)) + matrix(0.01, nrow = K, ncol = K)
prior_pred_dir(object = MOCK48_sub, dir_prior = dir_prior)

par(mfrow = c(3,4))
pc_prior <- matrix(NA, ncol = 3, K)
pc_prior[seq.int(1:K), ] <- matrix(rep(c(10, 60, 250), each = K), ncol = 3)

par(mfrow = c(4,4))
bandleres_beltran48 <- bandle(object_cond1 = c(beltrandyn48[[1]]),
                              object_cond2 = c(beltrandyn48[[2]]), params = hypres,
                              numIter = 20000, burnin = 10000, thin = 10, pg = FALSE, hypLearn = "LFBGS", 
                              maternCov = TRUE, nu = 2, pc_prior = pc_prior, numChains = 9, dir_prior = dir_prior)

save(bandleres_beltran48, file = "bandleres_beltran48.rda")