require(mcgraph)
require(ggplot2)
require(gridExtra)

source("../Comparative_analysis_of_Monte_Carlo_data/mcg.g2d.wrapper.R")
source("../Comparative_analysis_of_Monte_Carlo_data/huge.wrapper.R")
source("cluster.network.R")
res.raw <- readRDS("../Comparative_analysis_of_Monte_Carlo_data/results/results_2022-22-01-15:22:53")
res <- res.raw$data

###### PARAMETERS ######
iter <- c(15, 30, 60, 120)
prop <- c(0.05, 0.1, 0.2)
noise <- c(0.1, 1, 10)
ns <- c(100, 200, 1000)

# Thresholds
thresholds <- seq(0, 0.1, length.out=11)

# Plotting parameters
pal <- hcl.colors(10, palette="Spectral")[c(3, 9, 2)]
prop_dat <- data.frame()
noise_dat <- data.frame()
n_dat <- data.frame()

for (j in iter) {
    # prop
    for (i in prop) {
        G <- mcg.g2d.wrapper(nodes=clu_raw$nodes, graph="mcg.cluster", edges=clu_raw$edges,
            cluster=clu_raw$cluster, seed=clu_raw$seed, iter=j, prop=i)
        roc <- mcg.roc(d=G$data, G=G$theta, thresholds=thresholds, fun=mcg.lvs)
        prop_dat <- rbind(prop_dat, data.frame(sen = roc$PRC$sen, pre = roc$PRC$pre, prop = factor(rep(i, length(roc$PRC$sen))), iter = factor(rep(j, length(roc$PRC$sen)))))
    }
    # noise
    for (i in noise) {
        G <- mcg.g2d.wrapper(nodes=clu_raw$nodes, graph="mcg.cluster", edges=clu_raw$edges,
            cluster=clu_raw$cluster, seed=clu_raw$seed, iter=j, noise=i)
        roc <- mcg.roc(d=G$data, G=G$theta, thresholds=thresholds, fun=mcg.lvs)
        noise_dat <- rbind(noise_dat, data.frame(sen = roc$PRC$sen, pre = roc$PRC$pre, noise = factor(rep(i, length(roc$PRC$sen))), iter = factor(rep(j, length(roc$PRC$sen)))))
    }
    # n
    for (i in ns) {
        G <- mcg.g2d.wrapper(n = i, nodes=clu_raw$nodes, graph="mcg.cluster", edges=clu_raw$edges,
            cluster=clu_raw$cluster, seed=clu_raw$seed, iter=j)
        roc <- mcg.roc(d=G$data, G=G$theta, thresholds=thresholds, fun=mcg.lvs)
        n_dat <- rbind(n_dat, data.frame(sen = roc$PRC$sen, pre = roc$PRC$pre, n = factor(rep(i, length(roc$PRC$sen))), iter = factor(rep(j, length(roc$PRC$sen)))))
    }
}

# Generate plots

plt1 <- ggplot(data = prop_dat, aes(pre, sen, colour = prop, shape = prop)) +
    geom_abline(intercept = 1, slope = -1) +
    geom_point() +
    geom_line() +
    facet_wrap(vars(iter), nrow = 4, labeller = "label_both") +
    scale_colour_manual(values = pal) +
    theme_bw(base_size = 14) +
    theme(legend.position = "bottom") +
    labs(tag = bquote(~bold("A"))) +
    xlab("Sensitivity") +
    ylab("Precision") +
    coord_fixed()

plt2 <- ggplot(data = noise_dat, aes(pre, sen, colour = noise, shape = noise)) +
    geom_abline(intercept = 1, slope = -1) +
    geom_point() +
    geom_line() +
    facet_wrap(vars(iter), nrow = 4, labeller = "label_both") +
    scale_colour_manual(values = pal) +
    theme_bw(base_size = 14) +
    theme(legend.position = "bottom") +
    labs(tag = bquote(~bold("B"))) +
    xlab("Sensitivity") +
    ylab("Precision") +
    coord_fixed()

plt3 <- ggplot(data = n_dat, aes(pre, sen, colour = n, shape = n)) +
    geom_abline(intercept = 1, slope = -1) +
    geom_point() +
    geom_line() +
    facet_wrap(vars(iter), nrow = 4, labeller = "label_both") +
    theme_bw(base_size = 14) +
    scale_colour_manual(values = pal) +
    theme(legend.position = "bottom") +
    labs(tag = bquote(~bold("C"))) +
    xlab("Sensitivity") +
    ylab("Precision") +
    coord_fixed()

varParam <- grid.arrange(
    plt1, plt2, plt3,
    layout_matrix = rbind(c(1, 2, 3), c(1, 2, 3))
)
setwd("../figure")
ggsave("Figure6.pdf", width=10, height=10, varParam)
setwd("../Analysis_of_parameters_of_the_Monte_Carlo_algorithm")
