### Warning: some of the analyses below may take a long time to complete ###

library("dmbc")
library("hexbin")
library("ggplot2")
library("bayesplot")

# Section 4.2
data("simdiss", package = "dmbc")

plot(simdiss, colors = c("white", "darkgray"), cex.font = 0.75)

p <- 2
G <- 3
control <- dmbc_control(nsim = 10000, burnin = 20000, thin = 10, nchains = 2,
  threads = 2, parallel = "snow", seed = 601, verbose = TRUE)
prior <- dmbc_prior(eta = list(a = rep(1, G), b = rep(2, G)))
res <- dmbc(simdiss, p = p, G = G, control = control, prior = prior,
  post_all = TRUE)

# Section 4.4
summary(res, include.burnin = FALSE, summary.Z = FALSE)

est <- dmbc_get_configuration(res)
summary(est)

clusters(est)

set.seed(101)  # needed to replicate the label positioning
color_scheme_set(rep("black", 6))
graph <- plot(est, size = 2, size_lbl = 3, label_objects = TRUE)
graph + panel_bg(fill = "white", color = NA)

color_scheme_set("brightblue")
plot(res, what = "areas_ridges", regex_pars = "lambda", include.burnin = TRUE)

color_scheme_set("mix-blue-red")
plot(res, what = "trace", regex_pars = "eta", include.burnin = TRUE,
  transformation = "log", facet_args = list(labeller = ggplot2::label_parsed))

color_scheme_set("viridisA")
plot(res, what = "pairs", regex_pars = "alpha", diag_fun = "hist",
  off_diag_fun = "hex", include.burnin = TRUE)

# Section 4.5
pmax <- 3
Gmax <- 4
control <- list(burnin = 7000, nsim = 3000, thin = 10, z.prop = 1.5,
  alpha.prop = 0.75, random.start = TRUE, verbose = TRUE, store.burnin = TRUE,
  seed = 301)
ic <- dmbc_IC(data = simdiss, pmax = pmax, Gmax = Gmax, control = control)

# summarize results
print(ic)  # or simply ic

# plot results
color_scheme_set("darkgray")
plot(ic, size = c(4, 1.5)) +
  theme(panel.grid.major = element_line(colour = "gray",
    size = 0.05), panel.grid.minor = element_line(colour = "gray", size = 0.05),
    panel.border = element_rect(colour = "black", fill = NA),
    panel.background = element_blank())

# update the dmbc_ic object with an additional value of p
# (note: the following plot is not reported in the paper)
new.pmax <- pmax + 1
new.ic <- update(ic, pmax = new.pmax, Gmax = Gmax)
plot(new.ic, size = c(4, 1.5)) + 
  theme(panel.grid.major = element_line(colour = "gray",
    size = 0.05), panel.grid.minor = element_line(colour = "gray", size = 0.05),
    panel.border = element_rect(colour = "black", fill = NA),
    panel.background = element_blank())

# Section 5 - Empirical Application
data("animals", package = "dmbc")

plot(animals, colors = c("white", "darkgray"), cex.font = 0.75)

pmax <- 5
Gmax <- 5
prm.prop <- list(z = 2.5, alpha = 1)
control <- list(burnin = 10000, nsim = 10000, z.prop = prm.prop$z,
  alpha.prop = prm.prop$alpha, verbose = TRUE, seed = 201)
animals_ic <- dmbc_IC(animals, pmax = pmax, Gmax = Gmax, control = control)
animals_ic
# color_scheme_set('darkgray') plot(animals_ic)

# Solution: p = 1 -- G = 4
prm.prop <- list(z = 3, alpha = 1)
control <- list(burnin = 1e+05, nsim = 50000, thin = 10, z.prop = prm.prop$z,
  alpha.prop = prm.prop$alpha, verbose = TRUE, seed = 1406, nchains = 2,
  threads = 2, parallel = "snow")
animals_1_4 <- dmbc(animals, p = 1, G = 4, control = control,
  prior = dmbc_prior())

color_scheme_set("teal")

plot(animals_1_4, what = "trace", regex_pars = "eta", transformation = "log")
plot(animals_1_4, what = "trace", regex_pars = "sigma2", transformation = "log")
plot(animals_1_4, what = "trace", regex_pars = "lambda")
plot(animals_1_4, what = "trace", regex_pars = "alpha")

# cl1 <- clusters(dmbc_get_configuration(animals_1_4, chain = 1)) cl2 <-
# clusters(dmbc_get_configuration(animals_1_4, chain = 2)) table(cl1, cl2)

animals_lbl <- attr(animals@diss[[1]], "Labels")
z <- dmbc_get_configuration(animals_1_4, labels = animals_lbl)

color_scheme_set(rep("black", 6))
z_p <- plot(z, label_objects = TRUE, nudge_x = 1.5)
z_p + panel_bg(fill = "white", color = NA)

drop(z@Z.est)
clusters(z)