01_simulations_basic.Rmd

---
title: "Stochastic Block Model Prior with Ordering Constraints for Gaussian Graphical Models"
author:
    - Alessandro Colombi (Supervisor)^[a.colombi10@campus.unimib.it]
    - Teo Bucci^[teo.bucci@mail.polimi.it]
    - Filippo Cipriani^[filippo.cipriani@mail.polimi.it]
    - Filippo Pagella^[filippo.pagella@mail.polimi.it]
    - Flavia Petruso^[flavia.petruso@mail.polimi.it]
    - Andrea Puricelli^[andrea3.puricelli@mail.polimi.it]
    - Giulio Venturini^[giulio.venturini@mail.polimi.it]
output:
    pdf_document:
        toc: true
        toc_depth: 3
        number_section: true
    html_document:
        toc: true
        toc_float: true
        number_sections: true
#date: "2023-01-17"
editor_options:
    chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.align='center')
```

# Simulations

## Load the necessary packages

```{r, include = FALSE}
shhh = function(lib_name){ # It's a library, so shhh!
    suppressWarnings(suppressMessages(require(lib_name, character.only = TRUE)))
}
shhh("tidyverse")
shhh("ACutils")
shhh("mvtnorm")
shhh("salso")
shhh("FGM")
shhh("gmp")
shhh("mcclust")
shhh("mcclust.ext")
shhh("logr")
shhh("tidygraph")
shhh("ggraph")
shhh("igraph")
```

```{r, include = FALSE}
paths = c(
    "src/utility_functions.R",
    "src/bulky_functions.R",
    "src/data_generation.R"
)

for(p in paths){
    path = file.path(p)
    if(file.exists(path)){
        source(path)
    } else {
        cat("File",path,"was not found in directory, please check.")
    }
}
```


## Generate data

```{r}
# Define true clustering
rho_true = c(8,4,8,5)

# Set seed for data generation
seed_data_gen = 22111996

# Define number of observations
n = 500

# Define variance of the Beta
beta_sig2 = 1/16
```

```{r}
z_true = rho_to_z(rho_true)
r_true = z_to_r(z_true)
p = length(z_true)

# Generate data
sim = Generate_BlockDiagonal(n = n, z_true = z_true)
# sim = Generate_Block(
#     n=n,
#     z_true=z_true,
#     p_block_diag = 1,
#     p_block_extra_diag = 0,
#     p_inside_block = 0.95,
#     p_outside_block = 0.1,
#     elem_out = 5,
#     min_eigenval_correction = 3,
#     seed = 1
# )
```

```{r}
ACutils::ACheatmap(
    sim$Graph,
    use_x11_device = F,
    horizontal = F,
    main = "Graph",
    center_value = NULL,
    col.upper = "black",
    col.center = "grey50",
    col.lower = "white"
)

ACutils::ACheatmap(
    sim$Prec,
    use_x11_device = F,
    horizontal = F,
    main = "Precision matrix",
    center_value = NULL,
    col.upper = "black",
    col.center = "grey50",
    col.lower = "white"
)

ACutils::ACheatmap(
    sim$Cov,
    use_x11_device = F,
    horizontal = F,
    main = "Covariance",
    center_value = NULL,
    col.upper = "black",
    col.center = "grey50",
    col.lower = "white"
)
```

`sim$Prec` is the true precision matrix, `sim$Graph` is the true graph (adjacency matrix).

```{r}
graph_density = sum(sim$Graph) / (p*(p-1))
graph_density
```

## Set options for the simulation

```{r}
options = set_options(sigma_prior_0=0.5,
                      sigma_prior_parameters=list("a"=1,"b"=1,"c"=1,"d"=1),
                      theta_prior_0=1,
                      theta_prior_parameters=list("c"=1,"d"=1),
                      rho0=p, # start with one group
                      weights_a0=rep(1,p-1),
                      weights_d0=rep(1,p-1),
                      alpha_target=0.234,
                      beta_mu=graph_density,
                      beta_sig2=beta_sig2,
                      d=3,
                      alpha_add=0.5,
                      adaptation_step=1/(p*1000),
                      update_sigma_prior=TRUE,
                      update_theta_prior=TRUE,
                      update_weights=TRUE,
                      update_partition=TRUE,
                      update_graph=TRUE,
                      perform_shuffle=TRUE)
```

## Running the simulation

Create output directory if needed

```{r}
dir.create(file.path("output", "data"), showWarnings = FALSE, recursive = TRUE)
dir.create(file.path("output", "log"), showWarnings = FALSE, recursive = TRUE)
```

```{r}
unique_ID = uuid::UUIDgenerate(use.time = TRUE, n = 1, output = c("string"))
unique_ID = dittodb::hash(unique_ID, n = 8)
cat("This simulation has been assigned ID:", unique_ID)

filename_data = paste("output/data/simulation_", unique_ID, ".rds", sep = "")
filename_log = paste("output/log/simulation_", unique_ID, ".log", sep = "")
```

Run the simulation

```{r}
#log_open(file_name = filename_log, show_notes=FALSE, logdir = FALSE)
res <- Gibbs_sampler(
    data = sim$data,
    niter = 1000,
    nburn = 2,
    thin = 1,
    options = options,
    seed = 123456,
    print = TRUE
)
#log_close()
```

Before saving, also append true data
```{r}
res$true_rho = rho_true
res$true_precision = sim$Prec
res$true_graph = sim$Graph
# remove self loops
res$true_graph[col(res$true_graph)==row(res$true_graph)] = 0
```

```{r}
# save an object to a file
saveRDS(res, file = filename_data)
```

# Posterior analysis

Restore the object

```{r}
#unique_ID = 'c7726e0e'
filename_data = paste("output/data/simulation_", unique_ID, ".rds", sep = "")
res = readRDS(file = filename_data)
```

## Partition

Recomputing the partition in other forms and the number of groups

```{r}
rho_true = res$true_rho
r_true = rho_to_r(rho_true)
z_true = rho_to_z(rho_true)
p = length(z_true)
num_clusters_true = length(rho_true)
rho = res$rho
r = do.call(rbind, lapply(res$rho, rho_to_r))
z = do.call(rbind, lapply(res$rho, rho_to_z))
num_clusters = do.call(rbind, lapply(res$rho, length))
num_clusters = as.vector(num_clusters)
```

### Acceptance frequency

```{r}
mean(res$accepted)
```


### Barplot of changepoints

```{r}
bar_heights = colSums(r)
cp_true = which(r_true==1)
color <- ifelse(seq_along(bar_heights) %in% c(cp_true), "red", "gray")

barplot(
    bar_heights,
    names = seq_along(bar_heights),
    border = "NA",
    space = 0,
    yaxt = "n",
    main="Changepoint frequency distribution",
    #col = color,
    cex.names=.6,
    las=2
)

abline(v=cp_true-0.5, col="red", lwd=2)
legend("topright", legend=c("True"), col=c("red"),
    bty = "n",
    lty = 1,
    cex = 0.6)
```

### Evolution of the number of clusters

```{r}
plot(
    x = seq_along(num_clusters),
    y = num_clusters,
    type = "n",
    xlab = "Iterations",
    ylab = "Number of groups",
    main = "Number of groups - Traceplot"
)
lines(x = seq_along(num_clusters), y = num_clusters)
abline(h = length(z_to_rho(z_true)),
       col = "red",
       lwd = 4)
legend("topleft", legend=c("True"), col=c("red"),
    lty = 1,
    cex = 1)
```

```{r}
barplot(
    prop.table(table(num_clusters)),
    xlab = "Number of groups",
    ylab = "Relative Frequency",
    main = paste(
        "Number of groups - Relative Frequency\n",
        "Last:",
        tail(num_clusters, n = 1),
        "- Mean:",
        round(mean(num_clusters), 2),
        "- True:",
        num_clusters_true
    )
)
```

### Evolution of the Rand Index

```{r}
# computing rand index for each iteration
rand_index = apply(z, 1, mcclust::arandi, z_true)

# plotting the traceplot of the index
plot(
    x = seq_along(rand_index),
    y = rand_index,
    type = "n",
    xlab = "Iterations",
    ylab = "Rand Index",
    main = paste(
        "Rand Index - Traceplot\n",
        "Last:",
        round(tail(rand_index, n=1), 3),
        "- Mean:",
        round(mean(rand_index), 2)
    )
)
lines(x = seq_along(rand_index), y = rand_index)
abline(h = 1, col = "red", lwd = 4)
```

### Retrieving best partition using VI on visited ones (order is guaranteed here)

Here we are satisfied with finding the optimal partition only in the set of those visited, not in all the possible ones. I expect it could work even worse. But at least it guarantees to find an admissible one.
I would say that it is the implementation of formula (13) of the Corradin-Danese paper (https://doi.org/10.1016/j.ijar.2021.12.019).

```{r}
library("Rcpp")
library("RcppArmadillo")
sourceCpp("src/wade.cpp")

# compute VI
sim_matrix <- salso::psm(z)
dists <- VI_LB(z, psm_mat = sim_matrix)

# select best partition (among the visited ones)
best_partition_index = which.min(dists)
rho_est = rho[[best_partition_index]]
z_est = z[best_partition_index,]

# VI loss
dists[best_partition_index]

# select best partition
unname(z_est)

# compute Rand Index
mcclust::arandi(z_est, z_true)
```

## Graph

Extract last plinks

```{r}
last_plinks = tail(res$G, n=1)[[1]]
```

Criterion 1 to select the threshold (should not work very well) and assign final graph

```{r}
threshold = 0.5
G_est <- matrix(0,p,p)
G_est[which(last_plinks>threshold)] = 1
```

Criterion 2 to select the threshold

```{r}
bfdr_select = BFDR_selection(last_plinks, tol = seq(0.1, 1, by = 0.001))
```

Inspect the threshold and assign final graph

```{r}
bfdr_select$best_treshold
G_est = bfdr_select$best_truncated_graph
```

### Standardized Hamming distance

```{r}
SHD = sum(abs(sim$true_graph - G_est)) / (p^2 - p)
SHD
```

### Plot estimated matrices


```{r}
ACutils::ACheatmap(
    last_plinks,
    use_x11_device = F,
    horizontal = F,
    main = "Estimated plinks matrix",
    center_value = NULL,
    col.upper = "black",
    col.center = "grey50",
    col.lower = "white"
)

ACutils::ACheatmap(
    G_est,
    use_x11_device = F,
    horizontal = F,
    main = "Estimated Graph",
    center_value = NULL,
    col.upper = "black",
    col.center = "grey50",
    col.lower = "white"
)

ACutils::ACheatmap(
    tail(res$K,n=1)[[1]],
    use_x11_device = F,
    horizontal = F,
    main = "Estimated Precision matrix",
    center_value = NULL,
    col.upper = "black",
    col.center = "grey50",
    col.lower = "white"
)
```

### Evolution of the Kullback-Leibler

```{r}
kl_dist = do.call(rbind, lapply(res$K, function(k) {
    ACutils::KL_dist(res$true_precision, k)
}))

last = round(tail(kl_dist, n=1), 3)
plot(
    x = seq_along(kl_dist),
    y = kl_dist,
    type = "n",
    xlab = "Iterations",
    ylab = "K-L distance",
    main = paste("Kullback-Leibler distance\nLast value:", last)
)
lines(x = seq_along(kl_dist), y = kl_dist)
```

### Graph visualization

Plot against the original

```{r}
g1 <- graph.adjacency(res$true_graph)
edges1 <- get.edgelist(g1)
edges1 <- cbind(edges1, rep("true", nrow(edges1)))
g2 <- graph.adjacency(G_est)
edges2 <- get.edgelist(g2)
edges2 <- cbind(edges2, rep("estimated", nrow(edges2)))
edges <- as.data.frame(rbind(edges1, edges2))
names(edges) = c("from", "to", "graph")
nodes = data.frame(
    vertices = 1:p,
    clust_true = as.factor(z_true),
    clust_est = as.factor(z_est)
)
# nodes
g = graph_from_data_frame(edges, directed = FALSE, nodes)
lay = create_layout(g, layout = "linear", circular = TRUE)

output_plot <- ggraph(lay) +
    geom_edge_arc(edge_colour = "grey") +
    geom_node_point(aes(color = clust_true, shape = clust_est), size = 2) +
    geom_node_text(aes(label = name), repel = TRUE) +
    facet_edges( ~ graph)
print(output_plot)
```



```{r, eval=FALSE, include=FALSE}

# ================================================================
# =           DELETE FROM THIS TO THE END OF THE FILE?           =
# ================================================================

# read edges and nodes
g1 <- graph.adjacency(res$true_graph)
edges1 <- get.edgelist(g1)
edges1 <- cbind(edges1,rep("true",nrow(edges1)))

g2 <- graph.adjacency(final_graph)
edges2 <- get.edgelist(g2)
edges2 <- cbind(edges2,rep("estimated",nrow(edges2)))

edges <- as.data.frame(rbind(edges1,edges2))
names(edges) = c("from","to","graph")



nodes = data.frame(vertices=1:p,true_clust=as.factor(z_true), estim_clust=as.factor(final_partition_VI))
nodes

g = graph_from_data_frame(edges, directed = FALSE, nodes)

lay = create_layout(g, layout = "fr")

# plot with ggraph
ggraph(lay) + 
  geom_edge_link() + 
  geom_node_point() +
  theme_graph()

# unset graph theme 
# unset_graph_style()

# add node names
ggraph(lay) + 
  geom_edge_link(edge_colour = "grey") + 
  geom_node_point(aes(color=true_clust, shape = estim_clust), size = 4) +
  geom_node_text(aes(label = name), repel=TRUE)

ggraph(lay) + 
  geom_edge_link(edge_colour = "grey") + 
  geom_node_point(aes(color=true_clust, shape = estim_clust), size = 4) +
  geom_node_text(aes(label = name), repel=TRUE) +
  facet_edges(~graph)

library(magick)
library(gganimate)
anim <- 
  ggraph(lay) + 
  geom_edge_link(edge_colour = "grey") + 
  geom_node_point(aes(color=true_clust, shape = estim_clust), size = 4) +
  transition_states(graph,
                    transition_length = 0.5,
                    state_length = 0.5)


anim

```

```{r, eval=FALSE, include=FALSE}
p <- 25
estim_clust=as.factor(final_partition_VI)
nodes = data.frame(name=1:p,true_clust=as.factor(z_true), estim_clust = as.factor(z[1,]))

temp_g <- matrix(0,p,p)
temp_g[which(res$G[[1]]>threshold)] = 1
rownames(temp_g) <- 1:p
colnames(temp_g) <- 1:p
g1 <- graph.adjacency(temp_g)
g1 <- set_vertex_attr(g1, "name", index = V(g1), nodes$name)
g1 <- set_vertex_attr(g1, "tc", index = V(g1), nodes$true_clust)
g1 <- set_vertex_attr(g1, "ec", index = V(g1), as.factor(z[1,]))
graphs <- list(g1)
# read edges and nodes
for(i in 2:400){
    temp_g <- matrix(0,p,p)
    temp_g[which(res$G[[i]]>threshold)] = 1
    rownames(temp_g) <- 1:p
    colnames(temp_g) <- 1:p
    g1 <- graph.adjacency(temp_g)
    g1 <- set_vertex_attr(g1, "name", index = V(g1), nodes$name)
    g1 <- set_vertex_attr(g1, "tc", index = V(g1), nodes$true_clust)
    g1 <- set_vertex_attr(g1, "ec", index = V(g1), as.factor(z[i,]))
    temp_list <- list(g1)
    graphs <- append(graphs,temp_list)
}

#graphs
#x11(width= 16, height=9)
xy <- layout_as_dynamic(graphs, alpha = 0.85)
for(j in 1:length(xy)){
    for(i in 1:nrow(xy[[j]])){
        if(xy[[j]][i,1]<0){ xy[[j]][i,1] = -sqrt(-xy[[j]][i,1])}
        if(xy[[j]][i,2]>0){ xy[[j]][i,2] =  sqrt( xy[[j]][i,2])}
    }
}
xy
#pList <- vector("list", length(graphs))
#for (i in 1:length(graphs)) {
#  pList[[i]] <- ggraph(graphs[[i]], layout = "manual", x = xy[[i]][, 1], y = xy[[i]][, 2]) +
#    geom_edge_link0(edge_width = 0.6, edge_colour = "grey66") +
#    geom_node_point(shape = 21, aes(fill = as.factor(z_true)), size = 6) +
#    geom_node_text(label = 1:p, repel = FALSE, color = "white", size = 4) +
#    #scale_fill_manual(
#    #  values = c("forestgreen", "grey25", "firebrick"),
#    #  guide = ifelse(i != 2, "none", "legend"),
#    #  name = "smoking",
#    #  labels = c("never", "occasionally", "regularly")
#    #) +
#    theme_graph() +
#    theme(legend.position = "bottom") +
#    labs(title = paste0("Iter ", i))
#}
#Reduce("+", pList)


nodes_lst <- lapply(1:length(graphs), function(i) {
  cbind(igraph::as_data_frame(graphs[[i]], "vertices"),
    x = xy[[i]][, 1], y = xy[[i]][, 2], frame = i
  )
})


edges_lst <- lapply(1:length(graphs), function(i) cbind(igraph::as_data_frame(graphs[[i]], "edges"), frame = i))
edges_lst <- lapply(1:length(graphs), function(i) {
  edges_lst[[i]]$x <- nodes_lst[[i]]$x[match(edges_lst[[i]]$from, nodes_lst[[i]]$name)]
  edges_lst[[i]]$y <- nodes_lst[[i]]$y[match(edges_lst[[i]]$from, nodes_lst[[i]]$name)]
  edges_lst[[i]]$xend <- nodes_lst[[i]]$x[match(edges_lst[[i]]$to, nodes_lst[[i]]$name)]
  edges_lst[[i]]$yend <- nodes_lst[[i]]$y[match(edges_lst[[i]]$to, nodes_lst[[i]]$name)]
  edges_lst[[i]]$id <- paste0(edges_lst[[i]]$from, "-", edges_lst[[i]]$to)
  edges_lst[[i]]$status <- TRUE
  edges_lst[[i]]
})

all_edges <- do.call("rbind", lapply(graphs, get.edgelist))
all_edges <- all_edges[!duplicated(all_edges), ]
all_edges <- cbind(all_edges, paste0(all_edges[, 1], "-", all_edges[, 2]))

edges_lst <- lapply(1:length(graphs), function(i) {
  idx <- which(!all_edges[, 3] %in% edges_lst[[i]]$id)
  if (length(idx != 0)) {
    tmp <- data.frame(from = all_edges[idx, 1], to = all_edges[idx, 2], id = all_edges[idx, 3])
    tmp$x <- nodes_lst[[i]]$x[match(tmp$from, nodes_lst[[i]]$name)]
    tmp$y <- nodes_lst[[i]]$y[match(tmp$from, nodes_lst[[i]]$name)]
    tmp$xend <- nodes_lst[[i]]$x[match(tmp$to, nodes_lst[[i]]$name)]
    tmp$yend <- nodes_lst[[i]]$y[match(tmp$to, nodes_lst[[i]]$name)]
    tmp$frame <- i
    tmp$status <- FALSE
    edges_lst[[i]] <- rbind(edges_lst[[i]], tmp)
  }
  edges_lst[[i]]
})


edges_df <- do.call("rbind", edges_lst)
nodes_df <- do.call("rbind", nodes_lst)

p2 <- ggplot() +
  geom_segment(
    data = edges_df,
    aes(x = x, xend = xend, y = y, yend = yend, group = id, alpha = status),
    colour = "grey",
    show.legend = FALSE
  ) +
  geom_point(
    data = nodes_df, aes(x, y, group=name, fill = as.factor(ec)),
    shape = 21, size = 5, show.legend = FALSE
  ) +
  #geom_jitter()+
  #scale_fill_manual(values = c("forestgreen", "grey25", "firebrick")) +
  scale_alpha_manual(values = c(0, 1)) +
  ease_aes("quadratic-in") +
  enter_fade() + 
  exit_shrink() +
  transition_states(frame, transition_length = 2, wrap = FALSE) +
  labs(title = "Iter {closest_state}") +
  theme_void()
#animate(p2,nframes = 50*25,fps=50)
anim_save("prova.gif", p2, nframes=50*70, fps=50)
```




```{r, eval=FALSE, include=FALSE}
g = graph_from_data_frame(edges, directed = FALSE, nodes)
lay = create_layout(g, layout = "fr")

# add node names
ggraph(lay) + 
  geom_edge_link(edge_colour = "grey") + 
  geom_node_point(aes(color=true_clust, shape = estim_clust), size = 4) +
  geom_node_text(aes(label = name), repel=TRUE)

ggraph(lay) + 
  geom_edge_link(edge_colour = "grey") + 
  geom_node_point(aes(color=true_clust, shape = estim_clust), size = 4) +
  geom_node_text(aes(label = name), repel=TRUE) +
  facet_edges(~graph)


anim <- 
  ggraph(lay) + 
  geom_edge_link(edge_colour = "grey") + 
  geom_node_point(aes(color=true_clust, shape = estim_clust), size = 4) +
  transition_states(graph,
                    transition_length = 0.1,
                    state_length = 0.1)
```
ì