##################################################
## 1) Load data and packages
##################################################
rm(list = ls())
library(bipartite)
library(igraph)
library(networkD3)
library(readxl)
library(reshape2)

# Load the matrix from Excel
datos <- read_excel("data-METANETWORKS-QS.xlsx", sheet = "aggregated", col_names = TRUE)
agg_net <- as.matrix(datos[,-1])
rownames(agg_net) <- datos[[1]]
storage.mode(agg_net) <- "numeric"

cat("Matrix dimensions:", dim(agg_net), "\n")

# Create a bipartite graph with igraph
g <- graph_from_incidence_matrix(agg_net, weighted = TRUE)

# Create node data.frame
sites <- rownames(agg_net)
species <- colnames(agg_net)
nodes <- data.frame(nName = c(sites, species))
nodes$ID <- 0:(nrow(nodes)-1)
# Label type: "Site" or "Species"
nodes$type <- c(rep("Site", length(sites)), rep("Species", length(species)))
# Calculate degree
nodes$nodeDegree <- degree(g, mode = "all")

# Create edgeList for networkD3
edgeList <- melt(agg_net)
edgeList <- subset(edgeList, value > 0)
colnames(edgeList) <- c("Site", "Species", "Weight")
edgeList$SourceID <- match(edgeList$Site, nodes$nName) - 1
edgeList$TargetID <- match(edgeList$Species, nodes$nName) - 1

##################################################
## 2) Compute modularity with DIRT_LPA_wb_plus 100 times
##################################################
max_modularity <- -Inf
best_result <- NULL

for (i in 1:100) {
  res <- DIRT_LPA_wb_plus(agg_net)  # Compute modularity
  mod_temp <- convert2moduleWeb(agg_net, res)
  if (res$modularity > max_modularity) {
    max_modularity <- res$modularity
    best_result <- mod_temp
  }
}
cat("Best modularity:", max_modularity, "\n")

# Visualize partition (optional)
par(mar = c(10,10,8,5))
plotModuleWeb(best_result)

##################################################
## 3) Extract module information with listModuleInformation
##################################################
info_mod <- listModuleInformation(best_result)
# The structure of info_mod is:
#   [[1]]  (general details)
#   [[2]]  (list of modules)
str(info_mod)

# In your case, info_mod[[2]] is a list of 'k' modules.
# Each module is list(siteVector, speciesVector).

modules_list <- info_mod[[2]]
num_modules <- length(modules_list)
cat("Number of detected modules:", num_modules, "\n")

# Create membership vector for species
mod_membership_species <- rep(NA, length(species))
names(mod_membership_species) <- species

# Iterate over each module
for (mod_index in seq_along(modules_list)) {
  mod_info <- modules_list[[mod_index]]
  # mod_info[[1]] -> sites
  # mod_info[[2]] -> species in this module
  sp_vec <- mod_info[[2]]
  for (sp in sp_vec) {
    if (sp %in% names(mod_membership_species)) {
      mod_membership_species[sp] <- mod_index  # Assign module number
    }
  }
}

# Show module assignment table for species
cat("Module table (species only):\n")
print(table(mod_membership_species, useNA="ifany"))

##################################################
## 4) Assign modules to nodes in 'nodes'
##################################################
# For sites, use "Site"
membership_all <- rep(NA, nrow(nodes))
membership_all[nodes$type == "Site"] <- "Site"
# For species, assign what we set in mod_membership_species
for (i in 1:nrow(nodes)) {
  if (nodes$type[i] == "Species") {
    sp_name <- nodes$nName[i]
    membership_all[i] <- mod_membership_species[sp_name]
  }
}

nodes$Module <- membership_all

cat("Final module table (including 'Site'):\n")
print(table(nodes$Module, useNA="ifany"))

##################################################
## 5) Create interactive plot with networkD3
##################################################
# Color palette
unique_mods <- sort(unique(nodes$Module))
color_palette <- c("red", "blue", "green", "purple", "orange", "pink", "brown", "cyan", "gray")
if (length(unique_mods) > length(color_palette)) {
  color_palette <- colorRampPalette(color_palette)(length(unique_mods))
}
colourScale <- paste0("d3.scaleOrdinal().domain([\"",
                      paste(unique_mods, collapse = '", "'),
                      "\"]).range([\"",
                      paste(color_palette, collapse = '", "'),
                      "\"]);")

D3_network <- forceNetwork(Links = edgeList,
                           Nodes = nodes,
                           Source = "SourceID",
                           Target = "TargetID",
                           Value = "Weight",
                           NodeID = "nName",
                           Nodesize = "nodeDegree",
                           Group = "Module",
                           height = 1000,
                           width = 1000,
                           fontSize = 20,
                           linkWidth = JS("function(d){ return d.value/5; }"),
                           opacity = 0.85,
                           zoom = TRUE,
                           legend = TRUE,
                           charge = -30,
                           colourScale = JS(colourScale))

D3_network
saveNetwork(D3_network, "D3_verdadero.html", selfcontained = TRUE)