landscape.R

#' Make a 2D landscape for an Ising network
#'
#' Calculate the potential value \eqn{U(n)} for each system state, represented by the
#' number of active nodes \eqn{n}. The potential value is determined so that the Boltzmann
#' distribution is preserved. The Boltzmann distribution is the basis and the
#' steady-state distribution of all dynamic methods for Ising models, including
#' those used in [IsingSampler::IsingSampler()] and Glauber dynamics. This means
#' that if you assume the real-life system has the same steady-state distribution
#' as the Boltzmann distribution of the Ising model, then possibility that their
#' are \eqn{n} active nodes in the system is proportional to \eqn{e^{U(n)}}.
#' Because of this property of \eqn{e^{U(n)}}, it is aligned with the potential
#' landscape definition by Wang et al. (2008) and can quantitatively represent
#' the stability of different system states.
#'
#' The potential function \eqn{U(n)} is calculated by the following equation:
#' \deqn{U(n) = -\log(\sum_{v}^{a(v)=n} e^{-\beta H(v)})/\beta,}
#' where \eqn{v} represent a specific activation state of the network,
#' \eqn{a(v)} is the number of active nodes for \eqn{v}, and \eqn{H} is the
#' Hamiltonian function for Ising networks.
#'
#' @param thresholds,weiadj The thresholds and the weighted adjacency matrix
#' of the Ising network. If you have an `IsingFit` object estimated using
#' [IsingFit::IsingFit()], you can find those two parameters in its components
#' (`<IsingFit>$thresholds` and `<IsingFit>$weiadj`).
#' @param beta The \eqn{\beta} value for calculating the Hamiltonian.
#' @param transform By default, this function considers the Ising network
#' to use `-1` and `1` for two states. Set `transform = TRUE` if the Ising
#' network uses `0` and `1` for two states, *which is often the case for the
#' Ising networks estimated using* [IsingFit::IsingFit()].
#'
#' @return A `2d_Isingland` object that contains the following components:
#' \itemize{
#'   \item `dist_raw`,`dist` Two tibbles containing the probability
#'   distribution and the potential values for different states.
#'   \item `thresholds`,`weiadj`,`beta` The parameters supplied to the function.
#'   \item `Nvar` The number of variables (nodes) in the Ising network.
#' }
#' @seealso
#' [make_3d_Isingland()] if you have two groups of nodes that you want
#' to count the number of active ones separately.
#' @references
#' Wang, J., Xu, L., & Wang, E. (2008). Potential landscape and flux framework of nonequilibrium networks: Robustness, dissipation, and coherence of biochemical oscillations. Proceedings of the National Academy of Sciences, 105(34), 12271-12276. https://doi.org/10.1073/pnas.0800579105
#' Sacha Epskamp (2020). IsingSampler: Sampling methods and distribution functions for the Ising model. R package version 0.2.1. https://CRAN.R-project.org/package=IsingSampler
#' Glauber, R. J. (1963). Time-dependent statistics of the Ising model. Journal of Mathematical Physics, 4(2), 294-307. https://doi.org/10.1063/1.1703954
#'
#' @examples
#' Nvar <- 10
#' m <- rep(0, Nvar)
#' w <- matrix(0.1, Nvar, Nvar)
#' diag(w) <- 0
#' result1 <- make_2d_Isingland(m, w)
#' plot(result1)
#' @export
make_2d_Isingland <- function(thresholds, weiadj, beta = 1, transform = FALSE) {
  Nvar <- length(thresholds)

  # transformation
  if (transform) {
    vector.one <- matrix(1, Nvar, 1)
    thresholds <- 0.5 * thresholds + 0.25 * weiadj %*% vector.one
    weiadj <- 0.25 * weiadj
  }

  ## generate all states
  d <- tibble::tibble(v = all_combination(c(-1, 1), len = Nvar))

  ## calculate their H (or U)
  d <- d %>%
    dplyr::rowwise() %>%
    dplyr::mutate(
      U = as.numeric(ham(v, thresholds, weiadj)),
      sum_state = sum(v),
      freq = exp(-beta * U),
      n_active = (sum_state + Nvar) / 2
    ) %>%
    dplyr::ungroup()

  ## summarize based on the number of symptoms
  d_sum <- d %>%
    dplyr::group_by(sum_state) %>%
    dplyr::summarize(sum_freq = sum(freq)) %>%
    dplyr::mutate(
      n_active = (sum_state + Nvar) / 2,
      U = -log(sum_freq) / beta
    )

  result <- list(
    dist_raw = d, dist = d_sum, thresholds = thresholds,
    weiadj = weiadj, beta = beta, Nvar = Nvar
  )
  class(result) <- c("2d_Isingland", "Isingland", "landscape")
  return(result)
}

#' Make a matrix of landscapes for multiple Ising networks
#'
#' Make multiple landscapes together for different parameters.
#'
#' @param Ising_grid Parameter values generated by [make_Ising_grid()].
#' @inheritParams make_2d_Isingland
#'
#' @return A `2d_Isingland_matrix` object that contains the following
#' components:
#' \itemize{
#'   \item `dist_raw`,`dist` Two tibbles containing the probability
#'   distribution and the potential values for different states.
#'   \item `Nvar` The number of variables (nodes) in the Ising network.
#' }
#' @examples
#' Nvar <- 10
#' m <- rep(0, Nvar)
#' w <- matrix(0.1, Nvar, Nvar)
#' diag(w) <- 0
#' result4 <- make_Ising_grid(
#' all_thresholds(seq(-0.1, 0.1, 0.1), .f = `+`),
#' whole_weiadj(seq(0.5, 1.5, 0.5)),
#' m, w
#' ) %>% make_2d_Isingland_matrix()
#' plot(result4)
#' @export
make_2d_Isingland_matrix <- function(Ising_grid, transform = FALSE) {
  dist_raw <- Ising_grid %>%
    dplyr::rowwise() %>%
    dplyr::mutate(
      landscape = list(make_2d_Isingland(
        thresholds = thresholds_list,
        weiadj = weiadj_list,
        beta = beta_list,
        transform = transform
      )),
      dist = list(get_dist(landscape))
    ) %>%
    dplyr::ungroup()

  dist_tidy <- dist_raw %>%
    dplyr::select(dist, dplyr::all_of(attr(Ising_grid, "par_name"))) %>%
    tidyr::unnest(dist)

  return(structure(
    list(
      dist_raw = dist_raw,
      dist = dist_tidy,
      Nvar = dist_raw$landscape[[1]]$Nvar
    ),
    class = c("2d_Isingland_matrix", "Isingland", "landscape"),
    par_name = attr(Ising_grid, "par_name")
  ))
}

#' @export
print.landscape <- function(x, ...) {
	cat("A landscape object of the class", class(x)[1], "was estimated. Use `plot()` to draw the landscape plot.")
}

#' @export
plot.2d_Isingland <- function(x, ...) {
  ggplot2::ggplot(data = x$dist, ggplot2::aes(x = n_active, y = U)) +
    ggplot2::geom_point() +
    ggplot2::geom_line() +
    ggplot2::theme_bw() +
    ggplot2::xlab("Number of active nodes") +
		ggplot2::scale_x_continuous(breaks = seq(from = 0, to = x$Nvar, by = 3), minor_breaks = 1:x$Nvar)
}

#' @export
plot.2d_Isingland_matrix <- function(x, ...) {
  p <- ggplot2::ggplot(data = x$dist, ggplot2::aes(x = n_active, y = U)) +
    ggplot2::geom_point() +
    ggplot2::geom_line() +
    ggplot2::theme_bw() +
    ggplot2::xlab("Number of active nodes") +
  	ggplot2::scale_x_continuous(breaks = seq(from = 0, to = x$Nvar, by = 3), minor_breaks = 1:x$Nvar)

  if (length(attr(x, "par_name")) == 1) {
    p <- p + ggplot2::facet_wrap(attr(x, "par_name"))
  } else if (length(attr(x, "par_name")) == 2) {
    p <- p + ggplot2::facet_grid(attr(x, "par_name")) +
      ggplot2::scale_x_continuous(sec.axis = ggplot2::sec_axis(~., name = attr(x, "par_name")[2], breaks = NULL, labels = NULL)) +
      ggplot2::scale_y_continuous(sec.axis = ggplot2::sec_axis(~., name = attr(x, "par_name")[1], breaks = NULL, labels = NULL))
  }

  return(p)
}

#' Make a 3D landscape for an Ising network
#'
#' Similar to [make_2d_Isingland()], but two categories of nodes can
#' be specified so the number of active nodes can be calculated
#' separately.
#'
#' @inheritParams make_2d_Isingland
#' @param x,y Two vectors specifying the indices or the names of the
#' nodes for two categories. If they are character vectors, the names
#' should match the row names of the `thresholds` matrix.
#'
#' @seealso [make_2d_Isingland()] for the algorithm.
#'
#' @return A `3d_Isingland` object that contains the following components:
#' \itemize{
#'   \item `dist_raw`,`dist` Two tibbles containing the probability
#'   distribution and the potential values for different states.
#'   \item `thresholds`,`weiadj`,`beta` The parameters supplied to the function.
#'   \item `Nvar` The number of variables (nodes) in the Ising network.
#' }
#'
#' @export
make_3d_Isingland <- function(thresholds, weiadj, x, y, beta = 1, transform = FALSE) {
  l_2d <- make_2d_Isingland(thresholds, weiadj, beta, transform)
  d <- l_2d$dist_raw

  if (is.character(x)) {
  	if (!all(x %in% row.names(thresholds))) stop("The names in x are not found in the row names of the thresholds matrix.")
    x <- which(row.names(thresholds) %in% x)
  }
  if (is.character(y)) {
  	if (!all(y %in% row.names(thresholds))) stop("The names in y are not found in the row names of the thresholds matrix.")
    y <- which(row.names(thresholds) %in% y)
  }

  ## summarize based on the number of symptoms in the groups x & y
  d_sum <- d %>%
    dplyr::rowwise() %>%
    dplyr::mutate(
      sum_state_x = sum(v[x]),
      sum_state_y = sum(v[y])
    ) %>%
    dplyr::ungroup() %>%
    dplyr::group_by(sum_state_x, sum_state_y) %>%
    dplyr::summarize(sum_freq = sum(freq)) %>%
    dplyr::mutate(
      n_active_x = (sum_state_x + length(x)) / 2,
      n_active_y = (sum_state_y + length(y)) / 2,
      U = -log(sum_freq) / beta
    )

  result <- l_2d
  result$dist <- d_sum
  class(result) <- c("3d_Isingland", "Isingland", "landscape")
  return(result)
}

make_dist_matrix <- function(d, x_name, y_name, z_name) {
  x <- sort(unique(d[[x_name]]))
  y <- sort(unique(d[[y_name]]))
  z <- matrix(nrow = length(x), ncol = length(y))
  for (i in 1:nrow(d)) {
    z[
      which(x == as.numeric(d[i, x_name])),
      which(y == as.numeric(d[i, y_name]))
    ] <-
      as.numeric(d[i, z_name])
  }
  return(list(x = x, y = y, z = z))
}

#' @export
plot.3d_Isingland <- function(x, ...) {
  dm <- make_dist_matrix(x$dist, "sum_state_x", "sum_state_y", "U")
  p <- plotly::plot_ly(
    x = dm$x, y = dm$y,
    z = dm$z, type = "surface"
  )

  p <- plotly::layout(p, scene = list(
    xaxis = list(title = "Number of active nodes (x)"),
    yaxis = list(title = "Number of active nodes (y)"),
    zaxis = list(title = "U")
  )) %>%
    plotly::colorbar(title = "U")
  return(p)
}