02c_simulation_grid_visualization_notebook.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, cache=FALSE)
```

```{r, include = FALSE}
suppressWarnings(suppressPackageStartupMessages(library(tidyverse)))
suppressWarnings(suppressPackageStartupMessages(library(ACutils)))
suppressWarnings(suppressPackageStartupMessages(library(mvtnorm)))
suppressWarnings(suppressPackageStartupMessages(library(salso)))
suppressWarnings(suppressPackageStartupMessages(library(FGM)))
suppressWarnings(suppressPackageStartupMessages(library(gmp)))
suppressWarnings(suppressPackageStartupMessages(library(mcclust)))
suppressWarnings(suppressPackageStartupMessages(library(mcclust.ext)))
suppressWarnings(suppressPackageStartupMessages(library(logr)))
suppressWarnings(suppressPackageStartupMessages(library(tidygraph)))
suppressWarnings(suppressPackageStartupMessages(library(ggraph)))
suppressWarnings(suppressPackageStartupMessages(library(igraph)))
suppressWarnings(suppressPackageStartupMessages(library(pbapply)))
suppressWarnings(suppressPackageStartupMessages(library(latex2exp)))
suppressWarnings(suppressPackageStartupMessages(library(knitr)))
suppressWarnings(suppressPackageStartupMessages(library(kableExtra)))

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 {
        stop("File",path,"was not found in directory, please check.")
    }
}
```

```{r, echo = FALSE}
posterior_analysis <- function(i,group_id, recap_table, grid){
    if(i %in% (c(0,cumsum(c(7,10,2,4,7,10,1)))+1)){
        cat(paste("## Group",group_id,"\n\n"))
        group_id = group_id+1
    }
    cat(paste("### Simulation ID: ", grid[i,]$simulation_id),"\n")
    
    # print i-th recap table row
    
    print(knitr::kable(recap_table[i, ], booktabs = TRUE) %>%
        kable_styling(latex_options = c("hold_position", "scale_down")))
    
    # extract simulation ID
    simulation_id = grid[i,]$simulation_id
    
    # read current simulation
    sim <-
        read_rds(file.path(
            "output",
            "data",
            paste("simulation_", simulation_id, ".rds", sep = "")
        ))
    
    # compute other partition forms
    rho_true = sim$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 = sim$rho
    r = do.call(rbind, lapply(sim$rho, rho_to_r))
    z = do.call(rbind, lapply(sim$rho, rho_to_z))
    num_clusters = do.call(rbind, lapply(sim$rho, length))
    num_clusters = as.vector(num_clusters)
    
    # graph related quantities
    last_plinks = tail(sim$G, n=1)[[1]]
    bfdr_select = BFDR_selection(last_plinks, tol = seq(0.1, 1, by = 0.001))
    G_est = bfdr_select$best_truncated_graph # estimated adjacency
    
    # --------------------------------------------------------------------------
    
    # kl distance
    kl_dist = do.call(rbind, lapply(sim$K, function(k) {
        ACutils::KL_dist(sim$true_precision, k)
    }))
    
    # computing rand index for each iteration
    rand_index = apply(z, 1, mcclust::arandi, z_true)
    
    # 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,]
    
    if(options$graph_and_K){
        par(mfrow=c(1,4))
        ACutils::ACheatmap(
            sim$true_graph,
            use_x11_device = F,
            horizontal = F,
            main = "Graph",
            center_value = NULL,
            col.upper = "black",
            col.center = "grey50",
            col.lower = "white"
        )
        
        ACutils::ACheatmap(
            G_est,
            use_x11_device = F,
            horizontal = F,
            main = "\nEstimated Graph",
            center_value = NULL,
            col.upper = "black",
            col.center = "grey50",
            col.lower = "white"
        )
        
        ACutils::ACheatmap(
            sim$true_precision,
            use_x11_device = F,
            horizontal = F,
            main = "Precision matrix",
            center_value = NULL,
            col.upper = "black",
            col.center = "grey50",
            col.lower = "white"
        )
        
        ACutils::ACheatmap(
            tail(sim$K,n=1)[[1]],
            use_x11_device = F,
            horizontal = F,
            main = "\nEstimated Precision matrix",
            center_value = NULL,
            col.upper = "black",
            col.center = "grey50",
            col.lower = "white"
        )
        
        par(mfrow=c(1,1))
    }
    
    if(options$changepoint_kl){
        par(mfrow=c(1,3))
        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 = "\nChangepoint\n frequency distribution",
            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)
        
        # --------------------------------------------
        
        plot(
            x = seq_along(rand_index),
            y = rand_index,
            type = "n",
            xlab = "Iterations",
            ylab = "Rand Index",
            main = paste(
                "\nRand 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)
        
        # --------------------------------------------
        
        plot(
            x = seq_along(kl_dist),
            y = kl_dist,
            type = "n",
            xlab = "Iterations",
            ylab = "K-L distance",
            main = paste(
                "\nKullback-Leibler distance\nLast:",
                round(tail(kl_dist, n=1), 3)
                )
        )
        lines(x = seq_along(kl_dist), y = kl_dist)
        
        par(mfrow=c(1,1))
    }
    
    par(mar=c(5,4,4,2) + 0.1)
    
    if(options$theta_sigma_numgroupsfreq){
        par(mfrow=c(1,3))
        plot(
            x = seq_along(sim$sigma),
            y = sim$sigma,
            type = "n",
            xlab = "Iterations",
            ylab = TeX(r'($\sigma$ prior)'),
            main = TeX(r'($\sigma$ prior - Traceplot)')
        )
        lines(x = seq_along(sim$sigma), y = sim$sigma, lwd = 0.3)

        # --------------------------------------------
        
        plot(
            x = seq_along(sim$theta),
            y = sim$theta,
            type = "n",
            xlab = "Iterations",
            ylab = TeX(r'($\theta$ prior)'),
            main = TeX(r'($\theta$ prior - Traceplot)')
        )
        lines(x = seq_along(sim$theta), y = sim$theta, lwd = 0.3)
        
        # --------------------------------------------
        
        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
            )
        )
        
        par(mfrow=c(1,1))
    }
    
    if(options$numgroups_traceplot){
        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)
    }
    
    if(options$plot_graph){
        par(mfrow=c(1,1))
        
        # create graph for visualization
        g1 <- graph.adjacency(sim$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)
    }
    
    cat("\n\n\\pagebreak\n")
    #writeLines("ValueForV")
    return(group_id)
}
```

\newpage

# Recap table

```{r, echo = FALSE, include=FALSE, warning = FALSE}
# Read table
filename_data = "output/simulation_table.rds"
grid = readRDS(file = filename_data)

library("Rcpp")
library("RcppArmadillo")
sourceCpp("src/wade.cpp")
    
# build the final table row by row
recap_table = data.frame()
for(i in seq_len(nrow(grid))){
    
    # extract simulation ID
    simulation_id = grid[i,]$simulation_id
    
    # read current simulation
    sim <-
        read_rds(file.path(
            "output",
            "data",
            paste("simulation_", simulation_id, ".rds", sep = "")
        ))
    
    # compute other partition forms
    rho_true = sim$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 = sim$rho
    r = do.call(rbind, lapply(sim$rho, rho_to_r))
    z = do.call(rbind, lapply(sim$rho, rho_to_z))
    num_clusters = do.call(rbind, lapply(sim$rho, length))
    num_clusters = as.vector(num_clusters)
    
    # acceptance rate
    accept <- mean(sim$accepted)
    
    # estimated partition and indexes
    sim_matrix <- salso::psm(z)
    dists <- VI_LB(z, psm_mat = sim_matrix)
    best_partition_index = which.min(dists)
    rho_est = rho[[best_partition_index]]
    z_est = z[best_partition_index,]
    randindex = mcclust::arandi(z_est,z_true)
    VI_loss = dists[best_partition_index]
    
    # KL distance from the true precision to the last precision
    kl_dist = ACutils::KL_dist(sim$true_precision, sim$K[[length(sim$K)]])
    
    # execution time
    time = sim$execution_time
    
    # standardized Hamming distance
    last_plinks = tail(sim$G, n=1)[[1]]
    bfdr_select = BFDR_selection(last_plinks, tol = seq(0.1, 1, by = 0.001))
    G_est = bfdr_select$best_truncated_graph
    SHD = sum(abs(sim$true_graph - G_est)) / (p^2 - p)
    
    # create new row and append it
    new_row = data.frame(
        sim_id = grid[i,]$simulation_id,
        n = as.character(grid[i,]$n),
        p = as.character(grid[i,]$p),
        data_gen = grid[i,]$type_data_gen,
        seed = as.character(grid[i,]$seed_data_gen),
        rho0 = paste(grid[i,]$rho0,collapse=','),
        beta_sig2 = as.character(grid[i,]$beta_sig2),
        rho_true = paste(grid[i,]$rho_true,collapse=','),
        rho_est = paste(rho_est,collapse=','),
        accept = round(accept,3),
        VI = round(VI_loss,3),
        RI = round(randindex,3),
        KL = round(kl_dist,3),
        SHD = round(SHD,3),
        time = round(time,3)
    )
    
    recap_table = rbind(recap_table, new_row)
}
recap_table[19,]$rho0 <- "singletons"

# rename with latex format
#names(recap_table)[names(recap_table) == "sim_id"] <- "\\texttt{sim\\_id}"
#names(recap_table)[names(recap_table) == "data_gen"] <- "\\texttt{data\\_gen}"
#names(recap_table)[names(recap_table) == "seed"] <- "\\texttt{seed}"
#names(recap_table)[names(recap_table) == "time"] <- "\\texttt{time}"
#names(recap_table)[names(recap_table) == "accept"] <- "\\texttt{accept}"
#names(recap_table)[names(recap_table) == "beta_sig2"] <- "\\texttt{beta\\_sig2}"
#names(recap_table)[names(recap_table) == "rho0"] <- "$\\rho_0$"
#names(recap_table)[names(recap_table) == "rho_true"] <- "$\\rho_{\\text{true}}$"
#names(recap_table)[names(recap_table) == "rho_est"] <- "$\\rho_{\\text{est}}$"
#names(recap_table)[names(recap_table) == "n"] <- "$n$"
#names(recap_table)[names(recap_table) == "p"] <- "$p$"
```

```{r, echo=FALSE, results='asis'}
# export for LaTeX
#writeLines(knitr::kable(recap_table, "latex", booktabs=TRUE,escape=FALSE) %>% 
knitr::kable(recap_table, booktabs=TRUE) %>% 
    add_header_above(c(" ", "Data" = 7, "Analysis" = 7)) %>%
    pack_rows("Group 1: varying seed", 1, 7) %>%
    pack_rows("Group 2: varying beta_sig2", 8, 17) %>%
    pack_rows("Group 3: varying initial partition", 18, 19) %>%
    pack_rows("Group 4: varying group numerosities", 20, 23) %>%
    pack_rows("Group 5: varying n", 24, 30) %>%
    pack_rows("Group 6: varying p", 31, 40) %>%
    pack_rows("Group 7: using noised block structure", 41, 41) %>%
    kable_styling(latex_options=c("striped", "hold_position","scale_down"),
                  stripe_index = c(1:7,8:17,18:19,20:23,24:30,31:40,41))
#)
```

\newpage

# Simulations

```{r, echo = FALSE, fig.align='center', results = 'asis', fig.width=10, fig.height=2.8, fig.fullwidth=TRUE, warning=FALSE}

options <- list(graph_and_K = TRUE,
                plot_graph = TRUE,
                changepoint_kl = TRUE,
                numgroups_traceplot = FALSE,
                theta_sigma_numgroupsfreq = TRUE)
group_id = 1
for(i in 1:nrow(grid)){
    group_id <- posterior_analysis(i,group_id, recap_table, grid)
}
```