In adamsardar/stoneTrees: Solve Node Centric Steiner Tree Problems

library(ggplot2)
library(igraph)
library(purrr)
library(data.table)
library(stringr)
library(GGally)
library(foreach)

ScatterPlot <- function(dt,x,y){ 

                dt1 <- dt[,.(x = get(x), y = get(y)), by = eval(colnames(dt))]

                return(ggplot(dt1, aes(x=x, y=y)) +

                ggplot2::geom_point(shape=16, show.legend = FALSE,  alpha = 0.3, color = "paleturquoise4") +

                ggplot2::geom_abline(intercept = 0, slope = 1, col="darkred") +
                theme_minimal())}


GGscatterPlot <- function(data, mapping, ...) {


    x <- GGally::eval_data_col(data, mapping$x)
    y <- GGally::eval_data_col(data, mapping$y)

    dt <- data.table(Var1 = x, Var2 = y)

#Prepare plot
    pp <- ScatterPlot(dt = dt, x = "Var1", y = "Var2")

return(pp)
}

devtools::load_all()
fixedTerminalLymphomaGraph <- lymphomaGraph
V(fixedTerminalLymphomaGraph)$isTerminal <- FALSE
V(fixedTerminalLymphomaGraph)[nodeScore > 0]$isTerminal <- TRUE
V(fixedTerminalLymphomaGraph)$nodeScore <- -1

SFBenchDT <- rbind(fread("./results/Reproducible_steinForestSuboptimalBench_rcbc.tsv")%>%
  .[, `:=`(modulesVcount = strsplit(modulesVcount, "\\|"))],
                         fread("./results/Reproducible_steinForestSuboptimalBench_cplexAPI.tsv")%>%
  .[, `:=`(modulesVcount = strsplit(modulesVcount, "\\|"))])


SFBenchDT[, .(`TotalTime(hours)` = sum(time)/60/60), by = .(solver)]

SFBenchDT[, .(TotalTime = sum(time)/60/60), by = .(solver)][,TotalTime] %>% sum

SFBenchDT[, .(AvSise = mean(vcount)), by = .(solver)]

#
p1 <- ggplot(SFBenchDT, aes(x = solver, y = time/60)) +
  geom_boxplot() +
  geom_jitter(width = 0.1, aes(color = as.factor(trial))) +
  #scale_y_continuous( breaks = seq(0,100,5)) +
  theme_bw(base_size = 20) +
  theme(legend.position = "none")+
  labs(y = "Time (min)",
       title = "Time taken to solve Steiner forest problem")

p2 <- ggplot(SFBenchDT, aes(x = solver, y = time/60)) +
  geom_boxplot() +
  geom_jitter(width = 0.1, aes(color = as.factor(trial))) +
  #scale_y_continuous( breaks = seq(0,100,5)) +
  theme_bw(base_size = 20) +
  theme(legend.position = "none")+
  scale_y_log10()+
  labs(y = "Time (logScale)",
       title = "Time taken to solve Steiner forest problem")


#
p3 <- ggplot(SFBenchDT, aes(x = solver, y = vcount))+
  geom_boxplot() +
  geom_jitter(width = 0.1, aes(color = as.factor(trial))) +
  theme_bw(base_size = 20) +
  #scale_color_hue()+
  theme(legend.position = "none")+
  labs( y = "Aggregated module size",
        title = "Size of Steiner forest solution")


p4 <- SFBenchDT[,.(modulesVcount = unlist(modulesVcount)), by = .(solver, time, vcount, Niteration, trial)] %>%

  ggplot(., aes(x = solver, y = as.numeric(modulesVcount))) +
  geom_jitter(width = 0.1, aes(color = as.factor(trial)), alpha = 0.2) +
  geom_boxplot() +
  theme_bw(base_size = 20) +
  theme(legend.position = "none")+
  labs(y = "vCount of all modules in the pool",
       title = "Individual tree solutions to resampled terminals")


p5 <- ggplot(SFBenchDT, aes(x = solver, y = Niteration)) +
  geom_boxplot() +
  geom_jitter(width = 0.1, aes(color = as.factor(trial))) +
  theme_bw(base_size = 20) +
  theme(legend.position = "none")+
  labs(y = "Number of solutions",
       title = "Size of solution pool")


plot <- egg::ggarrange(p1,p3,p4,p5)



#ggsave(plot = plot, filename = "./results/CBCvsCPLEXvcGPLKPerformanceMStT.png", width = 16, height = 15)

allSolutions <- list.files("./results", pattern = ".RData", full.names = TRUE)

allSolutions <- allSolutions[grepl("Reproducible", allSolutions)]

allSolutionDT <- data.table(files = allSolutions)

allSolutionDT[,`:=`(graphPool = map(files, ~{ load(.x)
                                               SolutionPool[-length(SolutionPool)] } )), by = files]

#allSolutionDT[,`:=`(graphPool = map(files, ~{ load(.x)
#                                               SolutionPool } )), by = files]

allSolutionDT[,`:=`(poolSize = map(graphPool, length) %>% unlist,
                    collapsedGraph = map(graphPool, ~{ induced.subgraph(fixedTerminalLymphomaGraph, V(fixedTerminalLymphomaGraph)[map(.x, ~{ V(.x)$name }) %>% unlist() %>% unique()]) }),
                    Type =  str_match(files, "Version_(.*?)_")%>% .[1,2],
                    trial = str_match(files, "_Trial(.*?).RData")%>% .[1,2],
                    samplingMethod = str_extract(files, "MStT|steinForest")), by = files]

allSolutionDT[is.na(Type), `:=`(Type =  str_match(files, "Solver_(.*?)_")%>% .[1,2]), by = files]

allSolutionDT[,`:=`(eigenCentralityVector = map(collapsedGraph, ~{ eigen_centrality(.x)$vector  } )), by = files]


allSolutionDT[,`:=`(nodesVectorCollapsedGraph = map(collapsedGraph, ~{ V(.x)$name  } )), by = files]

allSolutionDT[,`:=`(nodesVectors = map(graphPool, ~{ map(.x, ~{ V(.x)$name  } ) } )), by = files]

eigenCentralityDT <- allSolutionDT[,.(eigenCentrality = unlist(eigenCentralityVector),
                 nodeID = names(unlist(eigenCentralityVector))), by = .(trial, Type, samplingMethod)]


eigenCentralityDT <- dcast(eigenCentralityDT,  nodeID+trial+samplingMethod~Type, value.var = "eigenCentrality", fill = 0) 

# New facet label names for dose variable
trial.labs <- str_c("trial ", eigenCentralityDT[,trial] %>% unique)
names(trial.labs) <- eigenCentralityDT[,trial] %>% unique


ScatterPlot(dt = eigenCentralityDT[samplingMethod == "steinForest"], x = "rcbc", y = "cplexAPI") +
  facet_wrap(~trial, ncol = 5, labeller = labeller(trial = trial.labs))+
  labs(y = "cplexAPI", x = "rcbc", title = "eigen_centrality") +
  theme_minimal(base_size = 20) 

#ggsave("./results/EigenCentralityCompareSolvers_CBC_CPLEX.png", width = 14, height = 10)

eigenCentralityDT <- allSolutionDT[,.(eigenCentrality = unlist(eigenCentralityVector),
                 nodeID = names(unlist(eigenCentralityVector))), by = .(trial, Type, samplingMethod)]

eigenCentralityDT <- dcast(eigenCentralityDT,  nodeID+Type+samplingMethod~trial, value.var = "eigenCentrality", fill = 0) 

SM <- "steinForest"

trialPlots <- map(eigenCentralityDT[samplingMethod == SM, Type] %>% unique, ~{

  if(SM == "MStT"){ DT <- eigenCentralityDT[Type == .x & samplingMethod == SM,.(nodeID, Type, samplingMethod, `1` = `1`, `2` = `2`, `3` = `3`)] }else{ DT <- eigenCentralityDT[Type == .x & samplingMethod == SM] }

ggpairs(DT, 
        columns = 3:ncol(DT), 
        lower = list(continuous = wrap(GGscatterPlot, method="pearson")),
        upper = list(continuous = wrap("cor", method= "pearson")),
        diag = "blank",
        columnLabels = str_c("Trial ", colnames(DT[, 3:ncol(DT)])),
        title = paste("- ", .x, "- Eigen-centrality (collapsed Solutions)"))

} )


#trialPlots[[1]]

#trialPlots[[5]]

#ggsave("./results/EigenCentralityCompareTrials_CBC.png", width = 17, height = 15)

distMethod <- "manhattan"

nodesDT <- allSolutionDT[,.(nodeID = unlist(nodesVectors, recursive = FALSE),
                           solution = seq(1:poolSize)), by = .(trial, Type, poolSize, samplingMethod)] %>%
  .[,.(nodeID = unlist(nodeID), value = 1), by = .(trial, Type, poolSize, samplingMethod, solution)]


nodesDT[,.N, by  = .(trial, Type, poolSize, samplingMethod, solution)]

SM <- "steinForest"

DistanceResults <- foreach::foreach(type = nodesDT[samplingMethod == SM,Type] %>% unique)%do%{

nodeVectorDT <- dcast(nodesDT[Type == type & samplingMethod == SM],  Type+trial +solution +poolSize+ samplingMethod~nodeID, fill = 0) 

distanceDT <- map(nodeVectorDT[,trial] %>% unique, ~{ 

nodesMatrix <- as.matrix(nodeVectorDT[trial == .x, 5:(ncol(nodeVectorDT))])

distance <- dist(nodesMatrix, method = distMethod, diag = FALSE, upper =FALSE) 

distanceDT <-  reshape::melt(as.matrix(distance), varnames = c("row", "col")) %>%
  as.data.table()

distanceDT <- distanceDT[row > col]

distanceDT[,`:=`(trial = .x, Type = type)]

return(distanceDT)

} ) %>% rbindlist()

return(distanceDT)

} %>% rbindlist()


ggplot(DistanceResults[grepl("1", trial)], aes(x=value, color = Type, fill = Type))+ 
  geom_histogram(bins =100, alpha=0.5,show.legend = T)+
  theme_minimal()+
  labs(y='', x= paste(distMethod, 'distance'))+
  facet_wrap(~trial, ncol = 5, labeller = labeller(trial = trial.labs))

ggplot(DistanceResults[!grepl("1", Type)], aes(x=value, color = trial, fill = trial))+ 
  geom_histogram(alpha=0.5, bins =100, show.legend = T)+
  theme_minimal()+
  labs(y='', x= paste(distMethod, 'distance'))+
  facet_wrap(~Type)