02e_simulation_grid_kl_comparison.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
        #keep_tex: yes
    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(
    dev = "pdf",
    echo = FALSE,
    cache = FALSE,
    fig.path = "output/",
    fig.align = 'center'
)
```

```{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(viridis)))
suppressWarnings(suppressPackageStartupMessages(library(knitr)))
suppressWarnings(suppressPackageStartupMessages(library(kableExtra)))
```
```{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.")
    }
}
```

\newpage

# Kullback-Leibler distances comparison

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

filename_data = file.path("output_noburnin", "simulation_table.rds")
grid = readRDS(file = filename_data)

kl_matrix = NULL

for(i in 1:nrow(grid)){
    simulation_id = grid[i,]$simulation_id
    sim <-
        read_rds(file.path(
            "output_noburnin",
            "data",
            paste("simulation_", simulation_id, ".rds", sep = "")
        ))
    
    kl_dist = do.call(rbind, lapply(sim$K, function(k) {
        ACutils::KL_dist(sim$true_precision, k)
    }))
    
    kl_matrix <- cbind(kl_matrix, kl_dist)
}

compare_kl_distances = function(kl_matrix, subset, grid, param, round=FALSE){
    
    labels= sapply(subset,function(x){grid[x,param]})
    
    if(round){
        labels = lapply(grid[subset, param], function(x){round(x,2)})
    }

    # Save current graphical parameters
    opar <- par(no.readonly = TRUE)
    
    # Change the margins of the plot (the fourth is the right margin)
    par(mar = c(5, 5, 4, 8))
    par(mfrow = c(1, 1))
    
    plot(
        x = seq_len(dim(kl_matrix[,subset])[1]),
        #y = kl_dist,
        type = "n",
        xlab = "Iterations",
        ylab = "K-L distance",
        main = "Kullback-Leibler distance",
        ylim = range(kl_matrix[,subset]),
        log = 'y'
    )
    
    vect_n_palette = seq_along(subset)
    color = viridis(length(subset))
    
    legend(
        "topright",
        title = param,
        legend = labels,
        fill = color,
        inset = c(-0.3, 0),
        xpd = TRUE
    )
    
    matlines(kl_matrix[,subset], col = color, lty = 1, lwd=2)
    grid()

}
```

## Varying $B_{\sigma^2}$

```{r kl_dist_comparison_beta_sig2, results = 'asis', out.width="50%", fig.height = 4.5}
compare_kl_distances(kl_matrix, 8:17, grid, "beta_sig2")
```

## Varying $p$

```{r kl_dist_comparison_p, results = 'asis', out.width="50%", fig.height = 4.5}
compare_kl_distances(kl_matrix, 31:40, grid, "p")
```

## Varying $n$

```{r kl_dist_comparison_n, results = 'asis', out.width="50%", fig.height = 4.5}
compare_kl_distances(kl_matrix, 24:30, grid, "n")
```