R/optimization.R

In MiloMonnier/supplynet: Supply network design and robustness assessment tools

# FIRST IMPROVEMENT ALGORITHM -------------------------------------------------

#' Graph Rewiring First Improvement Optimization Algorithm
#'
#' Minimize a given graph metric (function) by iterative rewiring process.
#' At each time step, a candidate is generated by rewiring
#'
#' @param g an igraph object
#' @param rewiring.FUN function of rewiring to apply to the graph at each time step
#' @param fun2min function, metric to minimize
#' @param niter maximal number of iterations. Default is Infinite
#' @param score.tolerance a numeric threshold of fun2min under which graph is
#' accepted as result. Default is 0.01 (1\% of a normalized metric)
#' @param max.succ.reject maximal number of successive rejects of candidates,
#' beyond which we consider that the graph cannot be improved anymore, and is
#' returned. Default is 100
#' @param verbose Wether or not displaying fun2min evolution. Default is TRUE
#' @param ... additionnal arguments passed to rewiring.FUN
#'
#' @return a list of 5 elements:
#' \describe{
#'   \item{g}{the optimized graph}
#'   \item{value}{the fun2min value obtained at the end of the process}
#'   \item{g_list}{list of igraphs object. History of successive improvements}
#'   \item{iter}{number of iterations elapsed}
#'   \item{acceptation_rate}{Acceptation rate during the process}
#' }
#' @export
#'
#' @examples
#' library(igraph)
#' g = sample_pa(n=100, power=1, m=1,  directed=FALSE)
#' is_dag(g)
#' res = firstImprovement(g, rewiring.FUN=rewire, with=each_edge(prob=1),
#' fun2min=transitivity, niter=10)
#'
firstImprovement = function(g,
                            rewiring.FUN      = NULL,
                            fun2min           = NULL,
                            niter             = Inf,
                            score.tolerance   = 0.01,
                            max.succ.reject   = 100,
                            verbose           = TRUE,
                            ...)
{
  # g = graph; f = function value; c = current; n = neighbor
  g_c = g_n = g
  f_c = f_n = fun2min(g)
  g_list = list() # List to save graph when improved only
  i = 0           # Nb of iterations
  k = 0           # Nb of successive rejects
  n_accepted = 0  # For output acceptation rate

  if (verbose) {
    message("Iter\tCurrent\tNeigh\tAcc.rate")
    message(sprintf("%i\t%.4f\t%.3f\t%.1f", i, f_c, f_n, 0, k))
  }

  repeat {
    # Rewire the current graph to generate a neighbour
    g_n = rewiring.FUN(g_c, ...)
    f_c = fun2min(g_c)
    f_n = fun2min(g_n)
    # If neighbour score is closer to objective, keep it
    if (f_n < f_c) {
      g_c = g_n
      k = 0
      n_accepted = n_accepted + 1
    }
    i = i + 1
    k = k + 1
    acceptation_rate = n_accepted / i
    if (verbose)
      message(sprintf("%i\t%.4f\t%.4f\t%.3f", i, f_c, f_n, acceptation_rate))
    # Stops when expected score is (almost) reached OR excessive nb of successive rejects
    if (f_c<score.tolerance  ||  k==max.succ.reject || i>=niter)
      break
  }
  # Return 1 graph and the dataframe of computed indices through iterations
  res = list(g=g_c,
             value=f_c,
             g_list=g_list,
             niter=i,
             acceptation_rate=acceptation_rate)
}
# TODO Compute indices to follow during evolution
# df_ind = rbind(df_ind, sapply(indices.FUN.list, mapply, list(g)))
# colnames(df_ind) = names(indices.FUN.list)
# https://stackoverflow.com/questions/30759367/apply-list-of-functions-to-list-of-values



# SIMULATED ANNEALING -----------------------------------------------------

#' Graph Rewiring Simulated Annealing Improvement Optimization Algorithm
#'
#' Try to minimize a given graph metric (function) through  rewiring process.
#' Metaheuristic trying to find a global optimum rather than local optimum.
#'
#' @param g an igraph object.
#' @param rewiring.FUN function of rewiring to apply to the graph at each time step.
#' @param fun2min function, metric to minimize.
#' @param niter maximal number of iterations. Default is Inf.
#' @param temp0 initial temperature. Default is 1.
#' @param cooling.method method used to progressively decrease temperature. Defaut is "exp".
#' @param cooling.factor decimal between 0 and 1, used if cooling.method="exp".
#' Lower cooling factor make temperature (and candidate acceptation probability)
#' decreasing faster. Default is 0.95.
#' @param verbose Wether or not displaying fun2min evolution. Default is FALSE.
#' @param ... additionnal arguments passed to rewiring.FUN
#'
#' @return a list of 7 elements:
#' \describe{
#'   \item{g}{the optimized graph}
#'   \item{value}{the fun2min value obtained at the end of the process}
#'   \item{g_list}{list of igraphs object. History of successive improvements}
#'   \item{iter}{number of iterations elapsed}
#'   \item{acceptation_rate}{Acceptation rate during the process}
#'   \item{temp0}{initial temperature passed as argument}
#'   \item{cooling.method}{cooling method passed as argument}
#' }
#' @export
#'
#' @importFrom stats runif
#' @examples
#' library(igraph)
#' g = sample_pa(n=100, power=1, m=1,  directed=FALSE)
#' is_dag(g)
#' res = simulatedAnnealing(g, rewiring.FUN=rewire, with=each_edge(prob=1),
#'  fun2min=transitivity, niter=10)
#'
simulatedAnnealing = function(g,
                              rewiring.FUN   = NULL,
                              fun2min        = NULL,
                              niter          = Inf,
                              temp0          = 1,
                              cooling.method = c("exp","fast","boltz"),
                              cooling.factor = 0.99,
                              verbose        = TRUE,
                              ...)
{
  cooling.method = match.arg(cooling.method)

  # g = graph; f = function value;  b = best; c = current; n = neighbor
  g_b = g_c = g_n = g
  f_b = f_c = f_n = fun2min(g_c)
  g_list = list() # List to save only improved graphs
  i = 0
  k = 0           # Annealing param = nb of iter until reannealing
  temp = temp0    # Initial temperature
  n_accepted = 0  # For saving acceptation rate

  if (verbose) {
    message("Iter\tBest\tCurrent\tNeigh\tTemp")
    message(sprintf("%i\t%.4f\t%.4f\t%.4f\t%.6f", i, f_b, f_c, f_n, temp))
  }

  repeat {
    # Rewire the current graph to generate a neighbour
    g_n = rewiring.FUN(g_c, ...)
    f_c = fun2min(g_c)
    f_n = fun2min(g_n)
    # If candidate is better than current OR probability to accept is high enough, keep it
    if (f_n < f_c  ||  exp(-(f_n-f_c)/temp) > runif(1)) {
      g_c = g_n
      k = 1   # Reset annealing parameter
      n_accepted = n_accepted + 1
      # Update best graph state
      if (f_c < f_b) {
        f_b = f_c
        g_b = g_c
        g_list = c(g_list, list(g_b))   # Store also into a list
        names(g_list)[length(g_list)] = i
      }
    }
    i = i + 1
    k = k + 1
    if (verbose)
      message(sprintf("%i\t%.4f\t%.4f\t%.4f\t%.6f", i, f_b, f_c, f_n, temp))
    # Lower the temperature = decreases probability to accept a worst neighbour
    if (cooling.method=="exp") {
      temp = temp * cooling.factor ^ k
    } else if (cooling.method=="fast") {
      temp = temp / k
    } else if (cooling.method=="boltz") {
      temp = temp / log(k) # Warning: 1/log(1) = Inf
    }
    # Stop when temperature reached 0 or max niter reached
    if (temp==1*10^-10 || i>=niter)
      break
  }
  # Return best graph and complementary elements
  res = list(g=g_b,
             value=f_b,
             g_list=g_list,
             niter=i,
             acceptation_rate = n_accepted/niter,
             temp0=temp0,
             cooling.method=cooling.method)
}