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


data("beltran2016HCMV24")
data("beltran2016MOCK24")
load("data/beltran24hypers.rda")

#par(mfrow = c(1,1))
#plot2D(beltran2016HCMV24, method = "t-SNE")
#plot2D(beltran2016MOCK24, method = "t-SNE")

idx <- intersect(rownames(beltran2016MOCK24), rownames(beltran2016HCMV24))

HCMV24_sub <- beltran2016HCMV24[idx, ]
MOCK24_sub <- beltran2016MOCK24[idx, ]


beltrandyn24 <- c(MOCK24_sub, HCMV24_sub) 
beltrandyn24 <- lapply(beltrandyn24, function(x) normalise(x, method = "center.median"))

K <- length(getMarkerClasses(MOCK24_sub))

dir_prior = diag(rep(1, K)) + matrix(0.01, nrow = K, ncol = K)
prior_pred_dir(object = MOCK24_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)
#hypres <- lapply(beltrandyn24, function(x) fitGPmaternPC(x, fcol = "markers", nu = 2, materncov = TRUE, hyppar = pc_prior)) 
#save(hypres, file = "beltran24hypers.rda")

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

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