# install.packages("gips")
library("gips")

# install.packages(c("ggplot2", "HSAUR2", "DAAG", "dplyr", "tidyr"))
library("ggplot2")

#####
# 1.2 aspirin

theme_set(theme_bw())

# This will add labels and numbers inside tiles
plot_cosmetic_modifications <- function(gg_plot_object) {
  my_col_names <- c(
    "deaths after placebo", "deaths after Aspirin",
    "treated with placebo", "treated with Aspirin"
  )

  suppressMessages( # message from ggplot2
    out <- gg_plot_object +
      scale_x_continuous(
        labels = my_col_names,
        breaks = 1:4
      ) +
      scale_y_reverse(
        labels = my_col_names,
        breaks = 1:4
      ) +
      theme(
        title = element_text(face = "bold", size = 18),
        legend.text = element_text(face = "bold", size = 10),
        axis.text.y = element_text(face = "bold", size = 11),
        axis.text.x = element_text(face = "bold", size = 11)
      ) +
      scale_fill_gradient2(
        low = "#F0EA3E", mid = "#A41836", high = "#95E956",
        midpoint = 3172432
      )
  )

  out +
    geom_text(aes(label = round(covariance, -3)),
      fontface = "bold", size = 8
    ) +
    theme(legend.position = "none")
}


data("aspirin", package = "HSAUR2")

Z <- aspirin

# Renumber the columns for better readability:
Z[, c(2, 3)] <- Z[, c(3, 2)]
names(Z) <- names(Z)[c(1, 3, 2, 4)]
rownames(Z) <- NULL

head(Z, 4)

n <- nrow(Z) # 7
p <- ncol(Z) # 4

S <- cov(Z)

g <- gips(S, n)
my_aspirin_id_ggplot <- plot_cosmetic_modifications(plot(g, type = "heatmap")) +
  labs(
    title = "Standard MLE estimator",
    subtitle = "of the covariance matrix"
  )
# Figure 1:
my_aspirin_id_ggplot


g_MAP <- find_MAP(g,
  optimizer = "brute_force",
  save_all_perms = TRUE, return_probabilities = TRUE
)
g_MAP

get_probabilities_from_gips(g_MAP)

compare_posteriories_of_perms(g_MAP, "(34)")
compare_posteriories_of_perms(g_MAP, "(12)")
compare_posteriories_of_perms(g_MAP, "()")

S_projected <- project_matrix(S, g_MAP)

my_aspirin_perm_ggplot <- plot_cosmetic_modifications(plot(g_MAP, type = "heatmap"))
# Figure 2:
my_aspirin_perm_ggplot


#####
# 1.2 books
theme_set(theme_bw())

# This will add labels and numbers inside tiles
plot_cosmetic_modifications <- function(gg_plot_object) {
  my_col_names <- c("thick", "height", "breadth")

  suppressMessages( # message from ggplot2
    out <- gg_plot_object +
      scale_x_continuous(
        labels = my_col_names,
        breaks = 1:3
      ) +
      scale_y_reverse(
        labels = my_col_names,
        breaks = 1:3
      ) +
      theme(
        title = element_text(face = "bold", size = 18),
        legend.text = element_text(face = "bold", size = 10),
        axis.text.y = element_text(face = "bold", size = 17),
        axis.text.x = element_text(face = "bold", size = 17)
      ) +
      scale_fill_gradient2(
        low = "#F0EA3E", mid = "#A41836", high = "#95E956",
        midpoint = 1.239099
      )
  )

  out +
    geom_text(aes(label = round(covariance, 1)),
      fontface = "bold", size = 8
    ) +
    theme(legend.position = "none")
}

data("oddbooks", package = "DAAG")

head(oddbooks, 4)

Z <- oddbooks[, c(1, 2, 3)]

number_of_observations <- nrow(Z) # 12
p <- ncol(Z) # 3

Z$height <- Z$height / sqrt(2) # A-series paper size have a height / breadth = sqrt(2)

S <- cov(Z)
g <- gips(S, number_of_observations)
my_books_id_ggplot <- plot_cosmetic_modifications(plot(g, type = "heatmap")) +
  labs(
    title = "Standard MLE estimator",
    subtitle = "of the covariance matrix"
  )
# Figure 3:
my_books_id_ggplot


g_MAP <- find_MAP(g,
  optimizer = "brute_force", show_progress_bar = FALSE,
  return_probabilities = TRUE, save_all_perms = TRUE
)

get_probabilities_from_gips(g_MAP)
g_MAP

my_books_map_ggplot <- plot_cosmetic_modifications(plot(g_MAP, type = "heatmap"))
# Figure 4:
my_books_map_ggplot


#####
# 2 Methodological Background


p <- 5
n <- 10

get_plotted_matrix <- function(my_matrix, my_title) {
  suppressMessages( # message from ggplot2
    gips:::pretty_plot_matrix(my_matrix, title = my_title) +
      geom_text(
        aes(label = round(covariance, 1)),
        fontface = "bold", size = 8
      ) +
      labs(fill = "covariance") +
      theme(
        legend.position = "none",
        plot.title = element_text(size = 25),
      ) +
      scale_fill_gradient2(
        low = "#F0EA3E", mid = "#A41836", high = "#95E956",
        midpoint = 8.660486
      )
  )
}


# Section 2.1

set.seed(2022)
sigma_matrix <- matrix(
  data = c(
    10, 08, 06, 06, 08,
    08, 10, 08, 06, 06,
    06, 08, 10, 08, 06,
    06, 06, 08, 10, 08,
    08, 06, 06, 08, 10
  ),
  nrow = p, byrow = TRUE
) # sigma_matrix is a matrix invariant under permutation (1,2,3,4,5)
Z <- MASS::mvrnorm(n, mu = rep(0, p), Sigma = sigma_matrix)
S <- cov(Z)


# Plots

# full symmetry
off_diagonal_element <- mean(S[row(S) != col(S)])
on_diagonal_element <- mean(S[row(S) == col(S)])
S_full <- matrix(rep(off_diagonal_element, p * p), nrow = p) + diag(on_diagonal_element - off_diagonal_element, p)
my_full_symmetry_ggplot <- get_plotted_matrix(
  S_full,
  bquote(atop("Example of a matrix" ~ bold("with full"), "permutation symmetry           "))
)
# Figure 5 top left:
my_full_symmetry_ggplot

# long perm symmetry
S_long <- project_matrix(S, "(12345)")

my_full_long_perm_ggplot <- get_plotted_matrix(
  S_long,
  bquote(atop("Example of a matrix" ~ bold("with long"), "permutation symmetry             "))
)
# Figure 5 top right:
my_full_long_perm_ggplot

# short perm symmetry
S_short <- project_matrix(S, "(123)")

my_short_perm_ggplot <- get_plotted_matrix(
  S_short,
  bquote(atop("Example of a matrix" ~ bold("with short"), "permutation symmetry              "))
)
# Figure 5 bottom left:
my_short_perm_ggplot

# no symmetry
my_no_perm_ggplot <- get_plotted_matrix(
  S,
  bquote(atop("Example of a matrix" ~ bold("without"), "permutation symmetry          "))
)
# Figure 5 bottom right:
my_no_perm_ggplot


# Section 2.2
g <- gips(S, n, perm = "(12345)", was_mean_estimated = FALSE)
summary(g)$n_parameters


# Section 2.3
g <- gips(S, n, perm = "(12345)", was_mean_estimated = FALSE)
summary(g)$n0

g <- gips(S, n, perm = "()", was_mean_estimated = FALSE)
summary(g)$n0

S_projected <- project_matrix(S, "(12345)")


# Section 2.4
g <- gips(S, n, perm = "(12345)", was_mean_estimated = FALSE)
exp(log_posteriori_of_gips(g))

g2 <- gips(S, n, perm = "(123)", was_mean_estimated = FALSE)
exp(log_posteriori_of_gips(g2))

compare_posteriories_of_perms(g, "(123)")


# Section 2.5
g <- gips(S, n, was_mean_estimated = FALSE)

set.seed(2022)
g_MAP_MH_25 <- find_MAP(g, max_iter = 25, optimizer = "MH") # message is anticipated
g_MAP_MH_25

g_MAP_BF <- find_MAP(g, optimizer = "BF")
g_MAP_BF

compare_posteriories_of_perms(g_MAP_BF, g_MAP_MH_25)

g_MAP_BF_with_probs <- find_MAP(g,
  optimizer = "BF",
  save_all_perms = TRUE, return_probabilities = TRUE
)
head(get_probabilities_from_gips(g_MAP_BF_with_probs), 10)

g_MAP_MH_20000 <- find_MAP(g,
  optimizer = "MH", max_iter = 20000,
  save_all_perms = TRUE, return_probabilities = TRUE
)
head(get_probabilities_from_gips(g_MAP_MH_20000), 10)


#####
# Appendix A

p <- 4
n <- 50

set.seed(2022)
Z <- matrix(rnorm(n * p), ncol = p)
S <- cov(Z)

g <- gips(S, n)
g_MAP <- find_MAP(g,
  optimizer = "BF", show_progress_bar = FALSE,
  return_probabilities = TRUE, save_all_perms = TRUE
)
get_probabilities_from_gips(g_MAP)


#####
# Section 3

# This part is commented out because it is computationally intensive.
# We used the cluster to run the simulations in "v112i07-sec3.R".
# source("v112i07-sec3.R")

sessionInfo()