# ========5.topological=======

#' Extract each sample network from the whole network

#'

#' @param whole_net the whole network

#' @param otutab otutab, columns are samples, these columns will be extract

#' @param threads threads, default: 1

#' @param save_net should save these sub_nets? FALSE or a filename

#' @param fast less indexes for faster calculate ?

#' @param verbose verbose

#' @param remove_negative remove negative edge or not? default: FALSE

#'

#' @return a dataframe contains all sub_net parameters

#' @export

#' @family topological

#' @examples

#' data(otutab, package = "pcutils")

#' extract_sample_net(co_net, otutab) -> sub_net_pars

extract_sample_net <- function(whole_net, otutab, threads = 1, save_net = FALSE, fast = TRUE, remove_negative = FALSE, verbose = TRUE) {

  i <- NULL

  V(whole_net)$name -> v_name

  reps <- ncol(otutab)

  if (verbose) message("extracting")

  sub_nets <- lapply(1:reps, \(i){

    rownames(otutab)[otutab[, i] > 0] -> exist_sp

    subgraph(whole_net, which(v_name %in% exist_sp)) -> spe_sub

    class(spe_sub) <- c("metanet", "igraph")

    return(spe_sub)

  })

  names(sub_nets) <- colnames(otutab)

  if (verbose) message("calculating topological indexes")

  # parallel

  # main function

  loop <- function(i) {

    spe_sub <- sub_nets[[i]]

    indexs <- net_par(spe_sub, mode = "n", fast = fast, remove_negative = remove_negative)[["n_index"]]

    wc <- igraph::cluster_fast_greedy(spe_sub, weights = abs(igraph::E(spe_sub)$weight))

    indexs$modularity <- igraph::modularity(wc)

    indexs

  }

  {

    if (threads > 1) {

      pcutils::lib_ps("foreach", "doSNOW", "snow", library = FALSE)

      if (verbose) {

        pb <- utils::txtProgressBar(max = reps, style = 3)

        opts <- list(progress = function(n) utils::setTxtProgressBar(pb, n))

      } else {

        opts <- NULL

      }

      cl <- snow::makeCluster(threads)

      doSNOW::registerDoSNOW(cl)

      res <- foreach::`%dopar%`(

        foreach::foreach(i = 1:reps, .options.snow = opts),

        loop(i)

      )

      snow::stopCluster(cl)

      gc()

    } else {

      res <- lapply(1:reps, loop)

    }

  }

  # simplify method

  sub_net_pars <- do.call(rbind, res)

  rownames(sub_net_pars) <- colnames(otutab)

  if (is.logical(save_net)) {

    if (save_net) save_net <- paste0("sub_net_", date())

  }

  if (is.character(save_net)) {

    saveRDS(sub_nets, file = paste0(save_net, ".RDS"))

  }

  sub_net_pars

}

#' Calculate natural_connectivity

#'

#' @param p an igraph or metanet object

#' @return natural_connectivity (numeric)

#' @export

#' @references \code{`nc` in `ggClusterNet`}

#' @family topological

#' @examples

#' igraph::make_ring(10) %>% nc()

nc <- function(p) {

  adj_matrix <- as.matrix(igraph::as_adj(p, sparse = FALSE))

  adj_matrix[abs(adj_matrix) != 0] <- 1

  lambda <- eigen(adj_matrix, only.values = TRUE)$values

  lambda <- sort(lambda, decreasing = TRUE)

  lambda_sum <- 0

  N <- length(lambda)

  for (i in 1:N) lambda_sum <- lambda_sum + exp(lambda[i])

  lambda_average <- log(lambda_sum / N, base = exp(1))

  lambda_average

}

#' Calculate all topological indexes of a network

#'

#' @param go an igraph or metanet object

#' @param fast less indexes for faster calculate ?

#' @param mode calculate what? c("v", "e", "n", "all")

#' @param remove_negative remove negative edge or not? default: FALSE

#'

#' @return a 3-elements list

#' \item{n_index}{indexs of the whole network}

#' \item{v_index}{indexs of each vertex}

#' \item{e_index}{indexs of each edge}

#' @export

#' @family topological

#' @examples

#' igraph::make_graph("Walther") %>% net_par()

#' c_net_index(co_net) -> co_net_with_par

net_par <- function(go, mode = c("v", "e", "n", "all"), fast = TRUE, remove_negative = FALSE) {

  from <- to <- NULL

  stopifnot(is_igraph(go))

  if ("all" %in% mode) mode <- c("v", "e", "n")

  n_index <- NULL

  v_index <- NULL

  e_index <- NULL

  Negative_percentage <- ifelse(!is.null(E(go)$cor), sum(igraph::E(go)$cor < 0) / length(igraph::E(go)), NA)

  # remove negative weight

  if (remove_negative) {

    if (!is.null(E(go)$cor)) {

      # message("Remove negative correlation edges")

      c_net_filter(go, cor > 0, mode = "e") -> go

    }

  }

  # non-weighted network

  up <- go

  if (!is.null(igraph::edge_attr(up)[["weight"]])) up <- igraph::delete_edge_attr(up, "weight")

  if ("n" %in% mode) {

    # Calculate Network Parameters

    n_index <- data.frame(

      check.names = F,

      `Node_number` = length(igraph::V(go)), # number of nodes

      `Edge_number` = length(igraph::E(go)), # number of edges

      `Edge_density` = igraph::edge_density(go), # density of network, connectance

      `Negative_percentage` = Negative_percentage, # negative edges percentage

      `Average_path_length` = igraph::average.path.length(up), # Average path length

      `Global_efficiency` = igraph::global_efficiency(up),

      `Average_degree` = mean(igraph::degree(go)), # Average degree

      `Average_weighted_degree` = ifelse(is.null(igraph::E(go)$weight), mean(igraph::degree(go)), sum(igraph::E(go)$weight) / length(igraph::V(go))), # weighted degree

      Diameter = igraph::diameter(up), # network diameter

      `Clustering_coefficient` = igraph::transitivity(go), # Clustering coefficient

      `Centralized_betweenness` = igraph::centralization.betweenness(go)$centralization, # Betweenness centralization

      `Natural_connectivity` = nc(go) # natural

    )

    if (!fast) {

      # mean_dist=mean_distance(go)#

      # w_mean_dist=ifelse(is.null(E(go)$weight),mean_dist,mean_distance(go))

      # v_conn= vertex.connectivity(go) #

      # e_conn= edge.connectivity(go) #

      # components= count_components(go) #

      modularity <- igraph::modularity(igraph::cluster_fast_greedy(go)) #

      rand.g <- igraph::erdos.renyi.game(length(V(go)), length(E(go)), type = "gnm")

      rand_m <- igraph::modularity(igraph::cluster_fast_greedy(rand.g))

      relative_modularity <- (modularity - rand_m) / rand_m #

      n_index <- data.frame(

        check.names = F,

        n_index,

        Modularity = modularity,

        `Relative_modularity` = relative_modularity,

        `Centralized_closeness` = igraph::centralization.closeness(go)$centralization, # Closeness centralization

        `Centralized_degree` = igraph::centralization.degree(go)$centralization, # Degree centralization

        `Centralized_eigenvector` = igraph::centralization.evcent(go)$centralization # eigenvector centralization

      )

    }

    n_index <- apply(n_index, 1, FUN = \(x)replace(x, is.nan(x), 0)) %>%

      t() %>%

      as.data.frame()

    n_index <- cbind_new(get_n(go, simple = TRUE), n_index)

  }

  if ("v" %in% mode) {

    # Calculate Vertices Parameters

    v_index <- data.frame(

      check.names = F,

      Degree = igraph::degree(go),

      `Clustering_coefficient` = igraph::transitivity(go, type = "local"), # local clustering coefficient

      Betweenness = igraph::betweenness(go), # betweenness

      Eccentricity = igraph::eccentricity(go),

      Closeness = igraph::closeness(go),

      `Hub_score` = igraph::hub_score(go)[["vector"]]

      # page_rank = page.rank(go)$vector

      # igraph::evcent(go)[["vector"]]

      # igraph::local_efficiency(go)

    )

    # weighted degree

    if (!is.null(E(go)$cor)) {

      get_e(go) -> edge_list

      edge_list %>%

        dplyr::select(from, cor) %>%

        rbind(., dplyr::select(edge_list, to, cor) %>% dplyr::rename(from = to)) %>%

        dplyr::group_by(from) %>%

        dplyr::summarise(w_degree = sum(cor)) -> w_degree

      v_index$`Average_weighted_degree` <- w_degree[match(rownames(v_index), w_degree$from), "w_degree"] %>% unlist()

    }

    v_index <- apply(v_index, 1, FUN = \(x)replace(x, is.nan(x), 0)) %>%

      t() %>%

      as.data.frame()

    v_index <- cbind_new(get_v(go), v_index)

  }

  if ("e" %in% mode) {

    # Calculate Edges Parameters

    e_index <- get_e(go)

    # if(!(edge_attr(go)%>%unlist()%>%is.null()))e_index=data.frame(edge_attr(go),e_index)

  }

  return(list(n_index = n_index, v_index = v_index, e_index = e_index))

}

#' Add topological indexes for a network

#' @param go igraph or metanet

#' @param force replace existed net_par

#'

#' @export

#' @rdname net_par

c_net_index <- function(go, force = FALSE) {

  if (!force) {

    if (!is.null(graph_attr(go)[["net_par"]])) stop("Already calculated net_pars, set `force = TRUE to replace existed net_par")

  }

  net_par(go, fast = FALSE) -> res

  graph_attr(go) <- as.list(res$n_index)

  graph_attr(go)[["net_par"]] <- TRUE

  vertex_attr(go) <- as.list(res$v_index)

  edge_attr(go) <- as.list(res$e_index)

  go

}

#' Fit power-law distribution for an igraph

#'

#' @param go igraph

#' @param p.value calculate p.value

#'

#' @return ggplot

#' @export

#' @family topological

#' @examples

#' fit_power(co_net)

fit_power <- function(go, p.value = FALSE) {

  x <- y <- formula <- NULL

  # igraph::degree distribution

  degree_dist <- table(igraph::degree(go))

  dat <- data.frame(degree = as.numeric(names(degree_dist)), count = as.numeric(degree_dist))

  # fit, set the original a & b

  mod <- stats::nls(count ~ a * degree^b, data = dat, start = list(a = 2, b = 1.5))

  summary(mod)

  # extract the coefficient

  a <- round(coef(mod)[1], 3)

  b <- round(coef(mod)[2], 3)

  fit <- fitted(mod)

  SSre <- sum((dat$count - fit)^2)

  SStot <- sum((dat$count - mean(dat$count))^2)

  R2 <- round(1 - SSre / SStot, 3)

  # bootstrap t get p.value

  if (p.value) {

    dat_rand <- dat

    p_num <- lapply(seq_len(999), \(i){

      dat_rand$count <- sample(dat_rand$count)

      SSre_rand <- sum((dat_rand$count - fit)^2)

      SStot_rand <- sum((dat_rand$count - mean(dat_rand$count))^2)

      R2_rand <- 1 - SSre_rand / SStot_rand

      R2_rand > R2

    })

    p_value <- (sum(unlist(p_num)) + 1) / (999 + 1)

  }

  p <- ggplot(dat, aes(x = degree, y = count)) +

    geom_point() +

    theme_bw() +

    stat_smooth(method = "nls", formula = y ~ a * x^b, method.args = list(start = list(a = 2, b = 1.5)), se = FALSE) +

    labs(x = "Degree", y = "Count")

  if (p.value) {

    label <- data.frame(

      x = 0.8 * max(dat$degree),

      y = c(0.9, 0.8, 0.7) * max(dat$count),

      formula = c(

        sprintf("italic(Y) == %.3f*italic(X)^%.3f", a, b),

        sprintf("italic(R^2) == %.3f", R2),

        sprintf("italic(P) < %.3f", p_value)

      )

    )

  } else {

    label <- data.frame(

      x = 0.8 * max(dat$degree),

      y = c(0.9, 0.8) * max(dat$count),

      formula = c(

        sprintf("italic(Y) == %.3f*italic(X)^%.3f", a, b),

        sprintf("italic(R^2) == %.3f", R2)

      )

    )

  }

  p + geom_text(aes(x = x, y = y, label = formula), data = label, parse = TRUE)

}

#' Degree distribution comparison with random network

#'

#' @param go igraph object

#' @param plot plot or not

#'

#' @return ggplot

#' @export

#' @family topological

#' @examples

#' rand_net(co_net)

rand_net <- function(go = go, plot = TRUE) {

  freq <- net <- NULL

  # generate a random network

  rand.g <- igraph::erdos.renyi.game(length(V(go)), length(E(go)), type = "gnm")

  if (!plot) {

    return(rand.g)

  }

  data1 <- data.frame(

    freq = igraph::degree_distribution(go), net = "Network",

    degree = 0:(length(degree_distribution(go)) - 1)

  )

  data2 <- data.frame(

    freq = igraph::degree_distribution(rand.g), net = "Random E-R",

    degree = 0:(length(degree_distribution(rand.g)) - 1)

  )

  # if data1[1,1]=0, it'is delete single vertex

  if (data1[1, 1] == 0) data1 <- data1[-1, ]

  data <- rbind(data1, data2)

  p1 <- ggplot(data) +

    geom_point(aes(x = degree, y = freq, group = net, fill = net), pch = 21, size = 2) +

    geom_smooth(aes(x = degree, y = freq, group = net, color = net), se = FALSE, method = "loess", formula = "y ~ x") +

    labs(x = "Degree", y = "Proportion") +

    scale_color_manual(values = c("#F58B8B", "#7AADF0")) +

    scale_fill_manual(values = c("#F58B8B", "#7AADF0")) +

    MetaNet_theme +

    theme(legend.position = c(0.8, 0.9), legend.title = element_blank())

  print(p1)

  return(rand.g)

}

#' Net_pars of many random network

#'

#' @param go igraph

#' @param reps simulation time

#' @param threads threads

#' @param verbose verbose

#'

#' @export

#' @rdname compare_rand

rand_net_par <- function(go, reps = 99, threads = 1, verbose = TRUE) {

  i <- NULL

  # parallel

  # main function

  loop <- function(i) {

    # generate a random network

    rand.g <- igraph::erdos.renyi.game(length(igraph::V(go)),

      length(igraph::E(go)),

      type = "gnm"

    )

    indexs <- net_par(rand.g, mode = "n")[["n_index"]]

    wc <- igraph::cluster_fast_greedy(rand.g)

    indexs$modularity <- igraph::modularity(wc)

    indexs

  }

  {

    if (threads > 1) {

      pcutils::lib_ps("foreach", "doSNOW", "snow", library = FALSE)

      if (verbose) {

        pb <- utils::txtProgressBar(max = reps, style = 3)

        opts <- list(progress = function(n) utils::setTxtProgressBar(pb, n))

      } else {

        opts <- NULL

      }

      cl <- snow::makeCluster(threads)

      doSNOW::registerDoSNOW(cl)

      res <- foreach::`%dopar%`(

        foreach::foreach(i = 1:reps, .options.snow = opts),

        loop(i)

      )

      snow::stopCluster(cl)

      gc()

    } else {

      res <- lapply(1:reps, loop)

    }

  }

  # simplify method

  rand_net_pars <- do.call(rbind, res)

  rand_net_pars

}

#' Compare some indexes between your net with random networks

#'

#' @param pars your net pars resulted by net_pars()

#' @param randp random networks pars resulted by rand_net_par()

#' @param index compared indexes: "Average_path_length","Clustering_coefficient" or else

#'

#' @return ggplot

#' @export

#' @family topological

#' @examples

#' data("c_net")

#' rand_net_par(co_net_rmt, reps = 30) -> randp

#' net_par(co_net_rmt, fast = FALSE) -> pars

#' compare_rand(pars, randp)

compare_rand <- function(pars, randp, index = c("Average_path_length", "Clustering_coefficient")) {

  V1 <- NULL

  labss <- t(pars$n_index[, index, drop = FALSE]) %>% as.data.frame()

  rownames(labss) -> labss$indexes

  p <- pcutils::group_box(randp[, index, drop = FALSE])

  p <- p +

    geom_hline(data = labss, aes(yintercept = V1), linetype = 2, color = "blue3") +

    geom_text(

      data = labss, aes(x = 1, y = V1 * 1.05, label = paste0("Network: ", round(V1, 3))),

      color = "blue3"

    ) +

    MetaNet_theme +

    theme(legend.position = "none", axis.text.x = element_blank())

  p

}

#' Calculate small-world coefficient

#'

#' @param go igraph or metanet

#' @param reps simulation time

#' @param threads threads

#' @param verbose verbose

#'

#' @return number

#' @export

#' @family topological

#' @examples

#' \donttest{

#' # set reps at least 99 when you run.

#' smallworldness(co_net, reps = 9)

#' }

smallworldness <- function(go, reps = 99, threads = 1, verbose = TRUE) {

  rand_net_par(go, reps = reps, threads = threads, verbose = verbose) -> rands

  small_world_coefficient <- (igraph::transitivity(go) / mean(rands$Clustering_coefficient)) /

    (igraph::average.path.length(go) / mean(rands$`Average_path_length`))

  small_world_coefficient

}