rewire.R

##
## wdnet: Weighted directed network
## Copyright (C) 2023  Yelie Yuan, Tiandong Wang, Jun Yan and Panpan Zhang
## Jun Yan <jun.yan@uconn.edu>
##
## This file is part of the R package wdnet.
##
## The R package wdnet is free software: You can redistribute it and/or
## modify it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or any later
## version (at your option). See the GNU General Public License at
## <https://www.gnu.org/licenses/> for details.
##
## The R package wdnet is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
##

#' @importFrom stats cor
#' @importFrom CVXR norm2
#' @importFrom utils modifyList
NULL

#' Degree preserving rewiring for directed networks
#'
#' Degree preserving rewiring towards the target structure \code{eta}.
#'
#' @param edgelist A two-column matrix, each row represents a directed edge from
#'   the first column to the second column.
#' @param eta A matrix generated by \code{wdnet::get_eta_directed()}.
#' @param iteration An integer, the number of rewiring iterations, with each
#'   iteration consisting of \code{nattempts} rewiring attempts.
#' @param nattempts An integer, the number of rewiring attempts for each
#'   iteration. Default value equals the number of rows in \code{edgelist}.
#' @param rewire.history Logical, whether the rewiring history should be
#'   returned.
#'
#' @return Rewired edgelist, degree based assortativity coefficients after each
#'   iteration, rewiring history (including the index of sampled edges and
#'   rewiring result). For each rewiring attempt, two rows are sampled form the
#'   edgelist, for example Edge1:(v_1, v_2) and Edge2:(v_3, v_4). If the
#'   rewiring attempt is accepted, the sampled edges are replaced as (v_1, v_4),
#'   (v_3, v_2).
#'
#' @keywords internal
#' 
dprewire_directed <- function(
    edgelist, eta,
    iteration = 200, nattempts,
    rewire.history = FALSE) {
  if (missing(nattempts) || is.null(nattempts)) nattempts <- nrow(edgelist)
  edgelist <- as.matrix(edgelist)
  snode <- edgelist[, 1]
  tnode <- edgelist[, 2]
  temp <- node_strength_cpp(
    snode = snode,
    tnode = tnode,
    nnode = max(edgelist),
    weight = 1,
    weighted = FALSE
  )
  outd <- temp$outs
  ind <- temp$ins

  sout <- outd[snode]
  sin <- ind[snode]
  tout <- outd[tnode]
  tin <- ind[tnode]

  df_s <- data.frame(
    type = rownames(eta),
    index = seq_len(nrow(eta)) - 1
  )
  df_t <- data.frame(
    type = colnames(eta),
    index = seq_len(ncol(eta)) - 1
  )
  type_s <- paste0(sout, "-", sin, split = "")
  type_t <- paste0(tout, "-", tin, split = "")

  index_s <- df_s[match(type_s, df_s$type), "index"]
  index_t <- df_t[match(type_t, df_t$type), "index"]
  rm(df_s, df_t, type_s, type_t, temp, outd, ind)

  ret <- dprewire_directed_cpp(
    iteration, nattempts,
    tnode,
    sout, sin,
    tout, tin,
    index_s, index_t,
    eta, rewire.history
  )
  rho <- data.frame(
    "Iteration" = c(0:iteration),
    "outout" = NA_real_,
    "outin" = NA_real_,
    "inout" = NA_real_,
    "inin" = NA_real_
  )
  rho[1, 2:5] <- c(
    "outout" = stats::cor(sout, tout),
    "outin" = stats::cor(sout, tin),
    "inout" = stats::cor(sin, tout),
    "inin" = stats::cor(sin, tin)
  )
  rho[2:(iteration + 1), 2] <- ret$outout
  rho[2:(iteration + 1), 3] <- ret$outin
  rho[2:(iteration + 1), 4] <- ret$inout
  rho[2:(iteration + 1), 5] <- ret$inin

  colnames(rho) <- c("Iteration", "outout", "outin", "inout", "inin")
  edgelist[, 2] <- ret$tnode
  result <- list(
    "assortcoef" = rho,
    "netwk" = create_wdnet(edgelist = edgelist, directed = TRUE),
    "iteration" = iteration,
    "nattempts" = nattempts
  )
  if (rewire.history) {
    colnames(ret$history) <- c("Attempt", "Edge1", "Edge2", "Accepted")
    ret$history[, 1:3] <- ret$history[, 1:3] + 1
    result$history <- ret$history
  }
  return(result)
}

#' Degree preserving rewiring for undirected networks
#'
#' Degree preserving rewiring towards the target structure \code{eta}.
#'
#' @param edgelist A two column matrix, each row represents an undirected edge.
#' @param iteration An integer, number of rewiring iterations, each iteration
#'   consists of \code{nattempts} rewiring attempts.
#' @param nattempts An integer, number of rewiring attempts for each iteration.
#'   The default value equals the number of rows in \code{edgelist}.
#' @param eta A matrix generated by \code{wdnet::get_eta_undirected()}.
#' @param rewire.history Logical, whether the rewiring history should be
#'   returned.
#' @return Rewired edgelist, assortativity coefficient after each iteration, and
#'   rewiring history (including the index of sampled edges and rewiring
#'   result). For each rewiring attempt, two rows are sampled from the
#'   \code{edgelist}, for example Edge1:\{v_1, v_2\} and Edge2:\{v_3, v_4\}, the
#'   function try to rewire the sampled edges as \{v_1, v_4\}, \{v_3, v_2\}
#'   (rewire type 1) or \{v_1, v_3\}, \{v_2, v_4\} (rewire type 2) with
#'   probability 1/2.
#'
#' @keywords internal
#' 
dprewire_undirected <- function(
    edgelist, eta,
    iteration = 200, nattempts,
    rewire.history = FALSE) {
  if (missing(nattempts) || is.null(nattempts)) nattempts <- nrow(edgelist)
  edgelist <- as.matrix(edgelist)
  degree <- data.frame(table(c(edgelist)))$Freq
  d_df <- data.frame(type = rownames(eta), index = seq_len(nrow(eta)) - 1)
  node1 <- edgelist[, 1]
  node2 <- edgelist[, 2]
  index1 <- d_df[match(degree[node1], d_df$type), "index"]
  index2 <- d_df[match(degree[node2], d_df$type), "index"]
  rm(d_df)
  degree1 <- degree[c(node1, node2)]
  degree2 <- degree[c(node2, node1)]
  ret <- dprewire_undirected_cpp(
    iteration, nattempts,
    node1, node2,
    degree1, degree2,
    index1, index2,
    eta, rewire.history
  )
  rm(node1, node2, degree1, degree2, index1, index2)
  rho <- data.frame("Iteration" = c(0:iteration), "Value" = NA_real_)
  rho[1, 2] <- assortcoef(edgelist = edgelist, directed = FALSE)
  rho[2:(iteration + 1), 2] <- ret$rho
  colnames(rho) <- c("Iteration", "Value")
  edgelist <- cbind(ret$node1, ret$node2)
  result <- list(
    "assortcoef" = rho,
    "netwk" = create_wdnet(edgelist = edgelist, directed = FALSE),
    "iteration" = iteration,
    "nattempts" = nattempts
  )
  if (rewire.history) {
    colnames(ret$history) <- c(
      "Attempt", "Edge1", "Edge2",
      "RewireType", "Accepted"
    )
    ret$history[, 1:4] <- ret$history[, 1:4] + 1
    result$history <- ret$history
  }
  return(result)
}

#' Degree preserving rewiring.
#'
#' Rewire a given network to have predetermined assortativity coefficient(s)
#' while preserving node degree.
#'
#' The algorithm first solves for an appropriate \code{eta} using
#' \code{target.assortcoef}, \code{eta.obj}, and \code{cvxr_control}, then
#' proceeds to the rewiring process and rewire the network towards the solved
#' \code{eta}. If \code{eta} is given, the algorithm will skip the first step.
#' This function only works for unweighted networks.
#'
#' Each rewiring attempt samples two rows from \code{edgelist}, for instance
#' Edge 1:(v_1, v_2) and Edge 2:(v_3, v_4). For directed networks, if the
#' rewiring attempt is accepted, the sampled edges are rewired as (v_1, v_4),
#' (v_3, v_2); for undirected networks, the algorithm try to rewire the sampled
#' edges as \{v_1, v_4\}, \{v_3, v_2\} (type 1) or \{v_1, v_3\}, \{v_2, v_4\}
#' (type 2), each with probability 1/2.
#'
#' @param netwk A \code{wdnet} object representing an unweighted network. If
#'   \code{NULL}, the function will construct a network using either
#'   \code{edgelist}, or \code{adj}.
#' @param edgelist A two column matrix, each row represents an edge of the
#'   network.
#' @param adj An adjacency matrix of an unweighted network.
#' @param directed Logical, whether the network is directed or not. It will be
#'   ignored if \code{netwk} is provided.
#' @param target.assortcoef For directed networks, it is a list represents the
#'   predetermined value or range of assortativity coefficients. For undirected
#'   networks, it is a constant between -1 to 1. It will be ignored if
#'   \code{eta} is provided.
#' @param control A list of parameters for controlling the rewiring process and
#'   the process for solving \code{eta}. \itemize{ \item{\code{iteration}} {An
#'   integer, represents the number of rewiring iterations. Each iteration
#'   consists of \code{nattempts} rewiring attempts. The assortativity
#'   coefficient(s) of the network will be recorded after each iteration.}
#'   \item{\code{nattempts}} {An integer representing the number of rewiring
#'   attempts for each
#'   iteration. Default value equals the number of rows of \code{edgelist}}.
#'   \item{\code{history}} {Logical, whether the rewiring attempts should be
#'   recorded and returned.} \item{\code{eta.obj}} {A convex function of
#'   \code{eta} to be minimized when solving for \code{eta} with given
#'   \code{target.assortcoef}. Defaults to 0. It will be ignored if \code{eta}
#'   is provided.} \item{\code{cvxr_control} {A list of parameters passed to
#'   \code{CVXR::solve()} for solving \code{eta} with given
#'   \code{target.assortcoef}. It will be ignored if \code{eta} is provided.}}}
#' @param eta A matrix represents the target network structure. If specified,
#'   \code{target.assortcoef} will be ignored. For directed networks, the
#'   element at row "i-j" and column "k-l" represents the proportion of directed
#'   edges linking a source node with out-degree i and in-degree j to a target
#'   node with out-degree k and in-degree l. For undirected networks, \code{eta}
#'   is symmetric, the summation of the elements at row "i", column "j" and row
#'   "j", column "i" represents the proportion of edges linking to a node with
#'   degree i and a node with degree j.
#'
#' @return Rewired network; assortativity coefficient(s) after each iteration;
#'   rewiring history (including the index of sampled edges and rewiring result)
#'   and solver results.
#'
#' @export
#'
#' @examples
#' \donttest{
#' set.seed(123)
#' netwk1 <- rpanet(1e4, control = rpa_control_scenario(
#'   alpha = 0.4, beta = 0.3, gamma = 0.3
#' ))
#' ## rewire a directed network
#' target.assortcoef <- list("outout" = -0.2, "outin" = 0.2)
#' ret1 <- dprewire(
#'   netwk = netwk1,
#'   target.assortcoef = target.assortcoef,
#'   control = list(iteration = 200)
#' )
#' plot(ret1$assortcoef$Iteration, ret1$assortcoef$"outout")
#' plot(ret1$assortcoef$Iteration, ret1$assortcoef$"outin")
#'
#' ## rewire an undirected network
#' netwk2 <- rpanet(1e4,
#'   control = rpa_control_scenario(
#'     alpha = 0.3, beta = 0.1, gamma = 0.3, xi = 0.3
#'   ),
#'   initial.network = list(
#'     directed = FALSE)
#' )
#' ret2 <- dprewire(
#'   netwk = netwk2,
#'   target.assortcoef = 0.3,
#'   control = list(
#'     iteration = 300, eta.obj = CVXR::norm2,
#'     history = TRUE
#'   )
#' )
#' plot(ret2$assortcoef$Iteration, ret2$assortcoef$Value)
#' }
#' 
dprewire <- function(
    netwk,
    edgelist, directed, adj,
    target.assortcoef = list(
      "outout" = NULL,
      "outin" = NULL,
      "inout" = NULL,
      "inin" = NULL
    ),
    control = list(
      "iteration" = 200,
      "nattempts" = NULL,
      "history" = FALSE,
      "cvxr_control" = cvxr_control(),
      "eta.obj" = function(x) 0
    ),
    eta) {
  netwk <- create_wdnet(
    netwk = netwk,
    edgelist = edgelist,
    edgeweight = NULL,
    directed = directed,
    adj = adj,
    weighted = FALSE
  )
  
  if (netwk$weighted) {
    warning("Edge weights are omitted")
  }

  control <- utils::modifyList(
    list(
      "iteration" = 200,
      "nattempts" = NULL,
      "history" = FALSE,
      "cvxr_control" = cvxr_control(),
      "eta.obj" = function(x) 0
    ), control,
    keep.null = TRUE
  )

  solver.result <- NULL
  if (missing(eta)) {
    if (netwk$directed) {
      solver.result <- get_eta_directed(
        edgelist = netwk$edgelist,
        target.assortcoef = target.assortcoef,
        eta.obj = control$eta.obj,
        control = control$cvxr_control
      )
    } else {
      stopifnot(
        is.numeric(target.assortcoef) && target.assortcoef >= -1 &&
          target.assortcoef <= 1
      )
      solver.result <- get_eta_undirected(
        edgelist = netwk$edgelist,
        target.assortcoef = target.assortcoef,
        eta.obj = control$eta.obj,
        control = control$cvxr_control
      )
    }
    if (solver.result$status == "solver_error" ||
      solver.result$status == "infeasible") {
      return(list("solver.result" = solver.result))
    }
    eta <- solver.result$eta
  }
  if (netwk$directed) {
    ret <- dprewire_directed(
      edgelist = netwk$edgelist,
      eta = eta,
      iteration = control$iteration,
      nattempts = control$nattempts,
      rewire.history = control$history
    )
  } else {
    ret <- dprewire_undirected(
      edgelist = netwk$edgelist,
      eta = eta,
      iteration = control$iteration,
      nattempts = control$nattempts,
      rewire.history = control$history
    )
  }
  ret$"solver.result" <- solver.result
  ret
}

#' Range of assortativity coefficients.
#'
#' The assortativity coefficient of a given network may not reach all the values
#' between -1 and 1 via degree preserving rewiring. This function calculates the
#' range of assortativity coefficients achievable through degree preserving
#' rewiring. The algorithm is designed for unweighted networks.
#'
#' The ranges are computed using convex optimization. The optimization problems
#' are defined and solved via the \code{R} package \code{CVXR}. For undirected
#' networks, the function returns the range of the assortativity coefficient.
#' For directed networks, the function computes the range of \code{which.range}
#' while other assortativity coefficients are restricted through
#' \code{target.assortcoef}.
#'
#' @param netwk A \code{wdnet} object representing an unweighted network. If
#'   \code{NULL}, the function will construct a network using either
#'   \code{edgelist} or \code{adj}.
#' @param edgelist A two-column matrix, where each row represents an edge of the
#'   network.
#' @param adj An adjacency matrix of an unweighted network.
#' @param directed Logical, whether the network is directed or not. It will be
#'   ignored if \code{netwk} is provided.
#' @param which.range The type of interested assortativity coefficient. For
#'   directed networks, it takes one of the values: "outout", "outin", "inout"
#'   and "inin". It will be ignored if the network is undirected.
#' @param target.assortcoef A list of constraints, it contains the predetermined
#'   value or range imposed on assortativity coefficients other than
#'   \code{which.range}. It will be ignored if the network is undirected.
#' @param control A list of parameters passed to \code{CVXR::solve()} for
#'   solving an appropriate \code{eta}, given the constraints
#'   \code{target.assortcoef}.
#'
#' @return Returns the range of the selected assortativity coefficient and the
#'   results from the solver.
#'
#' @export
#'
#' @examples
#' \donttest{
#' set.seed(123)
#' netwk <- rpanet(5e3,
#'   control =
#'     rpa_control_scenario(alpha = 0.5, beta = 0.5)
#' )
#' ret1 <- dprewire.range(
#'   netwk = netwk, which.range = "outin",
#'   target.assortcoef = list("outout" = c(-0.3, 0.3), "inout" = 0.1)
#' )
#' ret1$range
#' }
#' 
dprewire.range <- function(
    netwk,
    edgelist,
    adj,
    directed,
    which.range = c("outout", "outin", "inout", "inin"),
    control = cvxr_control(),
    target.assortcoef = list(
      "outout" = NULL,
      "outin" = NULL,
      "inout" = NULL,
      "inin" = NULL
    )) {
  netwk <- create_wdnet(
    netwk = netwk,
    edgelist = edgelist,
    edgeweight = NULL,
    directed = directed,
    adj = adj,
    weighted = FALSE
  )
  if (netwk$weighted) {
    warning("Edge weights are ignored.")
  }

  if (netwk$directed) {
    which.range <- match.arg(which.range)
    result <- get_eta_directed(
      edgelist = netwk$edgelist,
      target.assortcoef = target.assortcoef,
      which.range = which.range,
      control = control
    )
  } else {
    result <- get_eta_undirected(
      edgelist = netwk$edgelist,
      control = control
    )
  }
  result
}