## IMPORT PACKAGES ##
library(lattice)
library(grid)
library(XML)
library(igraph)
library(gridExtra)
library(ggplot2)
library(plyr)
library(dplyr)
library(vegan)
library(ggplot2)
library(reshape2)

## Import Core QPress functions
Graph_dir <- "/Users/mpmarzlo/Documents/Collaborations/201703_projetREMAIC/QualModels_final/"
Graph_dir <- "/Users/mpmarzlo/Documents/Collaborations/201703_projetREMAIC/QualModels_Frontiers/"
source(file.path(Graph_dir, 'PlotFunctions.R'), echo = FALSE)

graph_list <- list.files(Graph_dir, pattern = '2021.dia' )
setwd(Graph_dir)
invisible(list.files(pattern = '2021.dia'))

B1 = model.dia(("B1_english2021.dia")) 
B2 = model.dia(("B2_english2021.dia")) 
B3 = model.dia(("B3_english2021.dia")) 
N1 = model.dia(("N1_english2021.dia")) 
N2 = model.dia(("N2_english2021.dia")) 
N3 = model.dia(("N3_english2021.dia")) 
Re = model.dia(("Researcher_english2021.dia")) 
Sy = model.dia(("Synthetic_english2021.dia")) 

models_label = list(labelsB1 = levels(B1$From), labelsB2 = levels(B2$From), labelsB3 =levels(B3$From), labelsN1 = levels(N1$From),labelsN2 = levels(N2$From),labelsN3 = levels(N3$From), labelsSy = levels(Sy$From), labelsRe = levels(Re$From))
mod_list <- list('B1' = B1, 'B2' =B2, 'B3' = B3, 'N1'=N1, 'N2'=N2, 'N3'=N3, 'Sy'= Sy,'Re'=Re)

## PLOT MODELS FOR VISUAL CHECK
plot_graph = invisible(lapply(1:length(mod_list), function(x){
    g_B1 <- graph_from_data_frame(mod_list[[x]], directed = F)
    plot(g_B1,
    main = paste('model',names(mod_list[x])) )
} ))

## ADD DENSITY DEPENDENCE TO EACH MODEL NODE
B1 <- enforce.limitation(B1)
B2 <- enforce.limitation(B2)
B3 <- enforce.limitation(B3)
N1 <- enforce.limitation(N1)
N2 <- enforce.limitation(N2)
N3 <- enforce.limitation(N3)
Re <- enforce.limitation(Re)
Sy <- enforce.limitation(Sy)

mod_list <- list('B1' = B1, 'B2' =B2, 'B3' = B3, 'N1'=N1, 'N2'=N2, 'N3'=N3, 'Sy'= Sy,'Re'=Re)

scenariosList = list(
     list("Increase in Lease Area" = c("Lease_area" = +1)),
     list("Reduction in Farming and land uses" = c("Farming_and_land_uses" = -1)),
     list("Increase in Waste collection" = c("Waste_collection" = +1)),
     list("Improved shellfish farming management" = c("Shellfish_aqua_management" = +1)),
     list("Improved integrated management" = c("Integrated_management" = +1)),
     list("Decrease in Stock density" = c("Rearing_density" = -1))    
     )

scenarios = list(
     "Lease_area" = c("Lease_area" = +1),
     "Farming_and_land_uses" = c("Farming_and_land_uses" = -1),
     "Waste_collection" = c("Waste_collection" = +1),
     "Shellfish_aqua_management" = c("Shellfish_aqua_management" = +1),
     "Integrated_management" = c("Integrated_management" = +1),
     "Rearing_density" = c("Rearing_density" = -1)) 

# Model subset per scenarios
scen_mod = do.call('rbind', lapply(models_label, function(lab){
  names(scenarios) %in% lab 
}) )
row.names(scen_mod) = names(mod_list)
colnames(scen_mod) = names(scenarios)

output_var <-as.factor(c("Production", "Economic_return", "Pollutants_bacteria", "Social_acceptability", "Undersized_organic_waste", "Plankton", "Recreationnal_fishing", "Commercial_fishing"))

response_var <- as.factor(c("Plankton", "Nutrients", "Pollutants_bacteria","Undersized_organic_waste", "Production", "Rearing_density","Shellfish_quality","Economic_return","Lease_area","Social_acceptability" ))#,"Patrimonial_values"

response_3var <- as.factor(c("Undersized_organic_waste", "Economic_return","Social_acceptability" ))

include_synth <- F
include_research <- F
nb_list_var <- F

## AGGREGATE ALL GRAPH OBJECTS FOR MODEL STRUCTURE COMPARISON

for(i in 1:length(mod_list)){
  mod_list[[i]]$Model <- as.factor(names(mod_list)[i])
}
Agg_Model <- do.call('rbind', mod_list)
#summary(Agg_Model)

# Remove synthetic and/or scientists model for plot
if(! include_synth){
Agg_Model <- Agg_Model[Agg_Model$Model!='Sy',]}
if(! include_research){
Agg_Model <- Agg_Model[Agg_Model$Model!='Re',]}
Agg_Model$Model <- factor(Agg_Model$Model) 
#summary(Agg_Model)
nb_model = nlevels(Agg_Model$Model)
nb_model

# Barplot of Frequency of Variables
var_count  = as.data.frame(xtabs(~From, data = Agg_Model[Agg_Model[,'From']==Agg_Model[,'To'], ]))
var_count$Freq = as.numeric(var_count$Freq)
var_count = var_count[order(var_count$Freq, decreasing = T),]
var_count$From = factor(var_count$From, levels = as.character(var_count$From))
var_count$Var_Nb = seq(1, nrow(var_count))
var_count$Name =  as.factor(c(do.call('rbind', lapply(strsplit(as.character(var_count$From),  split='_'),function(x){paste(x, collapse = ' ')}))))
var_count$Name = factor(var_count$Name, levels = as.character(var_count$Name))

var_count$category = c('environmental')
# GREEN for following variables: Nutrients, Plankton, Pollutants_bacteria, Predation, Production, Unanticipated_losses, Undersized_organic_waste, Turbidity, Wild_shellfish
green_var =  c('Nutrients', 'Plankton', 'Pollutants_bacteria', 'Predation', 'Production', 'Unanticipated_losses', 'Undersized_organic_waste', 'Turbidity', 'Wild_shellfish')
var_count[var_count$From %in% green_var,'category'] = 'environmental'
# BLUE for following variables: Economic_return, Stock_density, Shellfish_quality, Area_concessions, Certification, Commercial_fishing, Market_price, Recreationnal_fishing, Running_costs
blue_var = c('Economic_return', 'Rearing_density', 'Shellfish_quality', 'Lease_area', 'Certification',  'Market_price', 'Running_costs','General_demand','External_production')
var_count[var_count$From %in% blue_var,'category'] = 'economic'
# RED for following variables:  Farming_and_land_uses, Shellfish_aqua_management, Social_acceptability, Integrated_management, Waste_collection, Patrimonial_values
red_var = c( 'Farming_and_land_uses', 'Shellfish_aqua_management', 'Social_acceptability', 'Integrated_management', 'Waste_collection', 'Patrimonial_values','Commercial_fishing','Recreationnal_fishing','Solid_waste')
var_count[var_count$From %in% red_var,'category'] = 'social'
var_count$category = factor(var_count$category, levels = c('social','environmental','economic'))

var_count_mod  = as.data.frame(xtabs(~From+Model, data = Agg_Model[Agg_Model[,'From']==Agg_Model[,'To'], ]))
var_count_mod$From = factor(var_count_mod$From, levels = as.character(var_count$From))
var_count_mod = var_count_mod[var_count_mod$Freq!=0,]
var_synth = merge (var_count, var_count_mod, by = 'From')



# Coloured Barplot
var_plot = var_synth[var_synth$Freq.x >1,]
var_plot$Name = factor(var_plot$Name)
levels(var_plot$Name) <- paste(seq(1,nlevels(var_plot$Name)), levels(var_plot$Name))


p <- ggplot(data = var_plot, aes(x = Name, y = Model, fill = Freq.y, col= category ) )+
  theme(legend.position="none", axis.text.x=element_text(angle=45,hjust=1)) +
  geom_tile(width = .8, height = .4, size = 2)+ 
  scale_fill_gradient(low='darkgrey')+
  xlab("Variables")+ylab("Models")+
  ggtitle("Variables Inclusion In Stakeholders' Models")
p

p <- ggplot(data = var_plot, aes(x = Model, y = Name, fill = Freq.y, col= category ) )+
  theme(legend.position="none") +
  geom_tile(width = .8, height = .4,  size = 1)+ 
  scale_fill_gradient(low='darkgrey')+
  xlab("Models")+ylab("Variables")+
  ggtitle("Variables Inclusion In Stakeholders' Models")+
  scale_y_discrete(limits = rev(levels(var_plot$Name)))

p0 <- p + theme(
plot.title = element_text(color="black", size=14, face="italic"),
axis.title.x = element_text(color="black", size=14, face="bold"),
axis.title.y = element_text(color="black", size=14, face="bold"),
plot.margin=unit(c(1,1,1,1.5),"cm") )
p0

# Distance between models

dist_mod = t(xtabs(~From+Model, data = Agg_Model[Agg_Model[,'From']==Agg_Model[,'To'], ]))
rownames(dist_mod)

dist_mod=(dist(t(xtabs(~From+Model, data = Agg_Model[Agg_Model[,'From']==Agg_Model[,'To'], ])), method = 'binary',diag = T, upper = T))
dist_mod = data.matrix(dist_mod)
dist_mod = melt(dist_mod)

ggplot(data = dist_mod, aes(x = Var1, y = Var2, fill = value )  )+  
  theme(legend.position="none", axis.text.x=element_text(angle=45,hjust=1)) +
  geom_tile(width = .8, height = .4)+ 
  scale_fill_gradient(low='white', high= 'black')

# Other Barplots about Variables Inclusion in Models

p <- ggplot( 
  data = var_count, aes(x = From, y = Freq ) ) +
  theme(axis.text.x=element_text(angle=45,hjust=1)) +
  geom_bar (stat= 'identity')+
  xlab("Variables")+ylab("Frequency")+
  ggtitle("Frequency of Variables Inclusion Across All Models")

p1 <- p + theme(
plot.title = element_text(color="red", size=14, face="italic"),
axis.title.x = element_text(color="blue", size=14, face="bold"),
axis.title.y = element_text(color="#993333", size=14, face="bold"),
plot.margin=unit(c(1,1,1,1.5),"cm") )

p <- ggplot( 
  data = var_count[var_count$Freq>1,], aes(x = From, y = Freq ) ) +
  theme(axis.text.x=element_text(angle=45,hjust=1)) +
  geom_bar (stat= 'identity')+
  xlab("Variables")+ylab("Frequency")+
  ggtitle("Frequency of Variables Inclusion Across All Models")

p2 <- p + theme(
plot.title = element_text(color="red", size=14, face="italic"),
axis.title.x = element_text(color="blue", size=14, face="bold"),
axis.title.y = element_text(color="#993333", size=14, face="bold"),
plot.margin=unit(c(1,1,1,1.5),"cm") )

gridExtra::grid.arrange(p1,p2, nrow = 2)

include_research <- F
nb_list_var <- F

## AGGREGATE ALL GRAPH OBJECTS FOR MODEL STRUCTURE COMPARISON

for(i in 1:length(mod_list)){
  mod_list[[i]]$Model <- as.factor(names(mod_list)[i])
}
Agg_Model <- do.call('rbind', mod_list)
#summary(Agg_Model)

# Remove scientists model for plot
if(! include_research){
Agg_Model <- Agg_Model[Agg_Model$Model!='Re',]}
Agg_Model$Model <- factor(Agg_Model$Model) 
Agg_Model$From <- factor(Agg_Model$From)
Agg_Model$To <- factor(Agg_Model$To)
#summary(Agg_Model)
nb_model = nlevels(Agg_Model$Model)
nb_model

# Change Variable Names to Numbers For Graphs
Agg_Model$NameFrom <- Agg_Model$From 
Agg_Model$NameTo <- Agg_Model$To
levels(Agg_Model$From) <- as.character(seq(1,nlevels(Agg_Model$From)))
levels(Agg_Model$To  ) <- as.character(seq(1,nlevels(Agg_Model$To)))

# Remove self limitation 
Agg_ModPlot <- Agg_Model[Agg_Model[,'From']!=Agg_Model[,'To'], ]

# Put Model Synthetic First in the list
Agg_ModPlot$Model <- factor(Agg_ModPlot$Model, levels = levels(Agg_ModPlot$Model)[c(7,1:6)] ) 
nb_model = nlevels(Agg_ModPlot$Model)


# Summarise information per model
nodes = as.data.frame( xtabs( ~ From, data= Agg_ModPlot ))
names(nodes)[2]  = "Nb_From"
nodes = cbind(nodes, list('Nb_To' = c(xtabs( ~ To, data= Agg_ModPlot ))) ) 
nodes = cbind(nodes, list('Model' = rowSums(xtabs( ~ From+Model, data= Agg_ModPlot )))) 
nodes = cbind(nodes, matrix(xtabs( ~ From+Model, data= Agg_ModPlot ), ncol= nb_model) ) 
nodes = cbind(nodes, matrix(xtabs( ~  To + Model, data= Agg_ModPlot), ncol= nb_model))   
names(nodes)[5:ncol(nodes)] <- c(paste0("From", 1:nb_model), paste0("To", 1:nb_model))

nodes$Model1 = nodes$From1 > 0 | nodes$To1 > 0
nodes$Model2 = nodes$From2 > 0 | nodes$To2 > 0
nodes$Model3 = nodes$From3 > 0 | nodes$To3 > 0
nodes$Model4 = nodes$From4 > 0 | nodes$To4 > 0
nodes$Model5 = nodes$From5 > 0 | nodes$To5 > 0
nodes$Model6 = nodes$From6 > 0 | nodes$To6 > 0
nodes$Model7 = nodes$From7 > 0 | nodes$To7 > 0
ind_start = 19
if(include_synth){
nodes$Model8 = nodes$From8 > 0 | nodes$To8 > 0
ind_start = 21}

nodes$ModelInclude = apply(nodes[,seq(ind_start,(ind_start+nb_model-1))], 1, FUN = function(x) paste(as.character(1:nb_model)[x], collapse= '_') )
nodes$Nb_Include = rowSums(nodes[,ind_start:(ind_start+nb_model-1)])
nodes$Model_graph = as.numeric(substr(nodes$ModelInclude, start=1, stop=1))
nodes$Model_graph[nodes$Nb_Include == 2] = 8

# RE-order variables in ascending order
row.names(nodes) = levels(Agg_ModPlot$NameFrom)
nodes = nodes[as.character(var_count$From[var_count$Freq >0]),]
nodes$From = 1:nrow(nodes)
#summary(nodes)
 
# Remove Duplicated Edges and Count Number of Occurrence
interactionCount = as.data.frame(xtabs( ~ From+To+Type, data= Agg_ModPlot ) ) 
Agg_ModPlot = Agg_ModPlot[! duplicated(Agg_ModPlot[,c('From','To','Type')]), ]
Agg_ModPlot = merge(Agg_ModPlot, interactionCount)
Agg_ModPlot$ModelInclus = factor(Agg_ModPlot$Model, levels= c(levels(Agg_ModPlot$Model),"Pair"))
Agg_ModPlot$ModelInclus[Agg_ModPlot$Freq==2] = "Pair"
Agg_ModPlot$ModelNb = Agg_ModPlot$ModelInclus
levels(Agg_ModPlot$ModelNb) = 1:nlevels(Agg_ModPlot$ModelNb)
#summary(Agg_ModPlot)

# Generate colors based on model
dir_graph= F
g_Agg = graph_from_data_frame(d= Agg_ModPlot, vertices= nodes, directed = dir_graph)

if(include_synth){
colrs = c("gray50", "tomato", "gold","lightblue","darkgreen","blanchedalmond","lawngreen","dodgerblue3","lightgrey")
}else{
  colrs = c("gray50", "tomato", "gold","lightblue","blanchedalmond","lawngreen","dodgerblue3","lightgrey")
}
# Set Node Colour based on Model
V(g_Agg)$color <- colrs[V(g_Agg)$Model_graph]
# Set node size based on Nb Of Inclusion
V(g_Agg)$size <- (V(g_Agg)$Nb_Include+10)*1.2
# Set Node Label colour
V(g_Agg)$label.color <- 'black'

# Set edge width based on weight:
if(dir_graph){
E(g_Agg)$width <- E(g_Agg)$Freq/2
E(g_Agg)$arrow.size <- .2
E(g_Agg)$edge.colour <- colrs[E(g_Agg)$ModelNb]
}else{
E(g_Agg)$width <- 1  
}


# Plot graphs

# Define legends positions
x_leg =  c(-1.7,1,1.4)
y_leg =  c(1.5, 1.3)  
cex_leg = c(1.2, .5)
nb_list_var = F

set.seed(29)
plot(g_Agg,
    vertex.label.cex = 0.8, 
  edge.width = E(g_Agg)$weight, 
    edge.arrow.size= .2,
    edge.curved=.1,
  layout=layout_with_lgl(graph=g_Agg, cellsize = 10e+7),
    main = paste('Comparing Model Structures (LGL projections)' ))
levels(Agg_Model$Model)

legend(x=x_leg[1], y= y_leg[1], levels(Agg_ModPlot$ModelInclus), pch=21,
       col="#777777", pt.bg=colrs, pt.cex=2, cex= cex_leg[1], bty="n", ncol=1, title = 'Model')

if(nb_list_var){
legend(x=x_leg[2], y= y_leg[2], paste(1:35,rownames(nodes)[1:35],sep =' - '), pch=NA,
       col="#777777", pt.bg=colrs, pt.cex=.5, cex= cex_leg[2], bty="n", ncol=1)
legend(x=x_leg[3], y= y_leg[2], paste(36:nrow(nodes),rownames(nodes)[36:nlevels(Agg_ModPlot$NameFrom)],sep =' - '), pch=NA,
       col="#777777", pt.bg=colrs, pt.cex=.5, cex= cex_leg[2], bty="n", ncol=1)
}

plot(g_Agg,
vertex.label.cex = 0.8, 
  edge.width = E(g_Agg)$weight, 
    edge.arrow.size= .2,
    edge.curved=.1,
  layout=layout_with_graphopt,
    main = paste('Comparing Model Structures (Graph-opt projection)' ))
legend(x=x_leg[1], y= y_leg[1], levels(Agg_ModPlot$ModelInclus), pch=21,
       col="#777777", pt.bg=colrs, pt.cex=4, cex= cex_leg[1], bty="n", ncol=1, title = 'Model')

if(nb_list_var){
legend(x=x_leg[2], y= y_leg[2], paste(1:35,rownames(nodes)[1:35],sep =' - '), pch=NA,
       col="#777777", pt.bg=colrs, pt.cex=.5, cex= cex_leg[2], bty="n", ncol=1)
legend(x=x_leg[3], y= y_leg[2], paste(36:nrow(nodes),rownames(nodes)[36:nlevels(Agg_ModPlot$NameFrom)],sep =' - '), pch=NA,
       col="#777777", pt.bg=colrs, pt.cex=.5, cex= cex_leg[2], bty="n", ncol=1)
}

plot(g_Agg,
    vertex.label.cex = 0.8, 
  edge.width = E(g_Agg)$weight, 
    edge.arrow.size= .2,
    edge.curved=.1,
  layout=layout_nicely,
    main = paste('Comparing Model Structures (aesthetic layout)' ))
legend(x=x_leg[1], y= y_leg[1], levels(Agg_ModPlot$ModelInclus), pch=21,
       col="#777777", pt.bg=colrs, pt.cex=4, cex= cex_leg[1], bty="n", ncol=1, title = 'Model')

if(nb_list_var){
legend(x=x_leg[2], y= y_leg[2], paste(1:35,rownames(nodes)[1:35],sep =' - '), pch=NA,
       col="#777777", pt.bg=colrs, pt.cex=.5, cex= cex_leg[2], bty="n", ncol=1)
legend(x=x_leg[3], y= y_leg[2], paste(36:nrow(nodes),rownames(nodes)[36:nlevels(Agg_ModPlot$NameFrom)],sep =' - '), pch=NA,
       col="#777777", pt.bg=colrs, pt.cex=.5, cex= cex_leg[2], bty="n", ncol=1)
}

plot(g_Agg,
    vertex.label.cex = 0.8, 
  edge.width = E(g_Agg)$weight, 
    edge.arrow.size= .2,
    edge.curved=.1,
  layout=layout_on_sphere,
    main = paste('Comparing Model Structures (sphere layout)' ))
legend(x=x_leg[1], y= y_leg[1], levels(Agg_ModPlot$ModelInclus), pch=21,
       col="#777777", pt.bg=colrs, pt.cex=4, cex= cex_leg[1], bty="n", ncol=1, title = 'Model')

if(nb_list_var){
legend(x=x_leg[2], y= y_leg[2], paste(1:35,rownames(nodes)[1:35],sep =' - '), pch=NA,
       col="#777777", pt.bg=colrs, pt.cex=.5, cex= cex_leg[2], bty="n", ncol=1)
legend(x=x_leg[3], y= y_leg[2], paste(36:nrow(nodes),rownames(nodes)[36:nlevels(Agg_ModPlot$NameFrom)],sep =' - '), pch=NA,
       col="#777777", pt.bg=colrs, pt.cex=.5, cex= cex_leg[2], bty="n", ncol=1)
}


print( paste(1:nlevels(Agg_ModPlot$NameFrom), levels(Agg_ModPlot$NameFrom)) )

plot.new()
if(!nb_list_var){
legend(x=-0.1, y= .9, paste(1:35,rownames(nodes)[1:35],sep =' - '), pch=NA,
       col="#777777", pt.bg=colrs, pt.cex=.9, cex= cex_leg[2], bty="n", ncol=1)
legend(x=.5, y= .9,  paste(36:nrow(nodes),rownames(nodes)[36:nlevels(Agg_ModPlot$NameFrom)],sep =' - '), pch=NA,
       col="#777777", pt.bg=colrs, pt.cex=.9, cex= cex_leg[2], bty="n", ncol=1)
}

# Generate plots
plot_matrix = lapply(1:length(scenarios), function(i){
selecMod = scen_mod[,names(scenarios[[i]])]

plotResults(mod_list = mod_list[selecMod], 
            scenario_list = scenariosList[i],
            label_list=  models_label[selecMod],
            output_var = response_var,
            label_x = names(mod_list[selecMod]),
            Scenario_name = names(scenariosList[[i]])) 
})


# Display results

plot_matrix[[1]]

plot_matrix[[2]]

plot_matrix[[3]]

plot_matrix[[4]]

plot_matrix[[5]]

plot_matrix[[6]]

sessionInfo()