R/noisyCE2.R

In noisyCE2: Cross-Entropy Optimisation of Noisy Functions

#' Cross-Entropy Optimisation of Noisy Functions
#' 
#' Unconstraint optimisation of noisy functions through the cross-entropy
#' algorithm.
#'
#' @param f objective function which takes the vector of optimisation variables
#'   as first argument.
#' @param domain a `list` (or other coercible objects) where each component
#'   specifies the domain of each variable of the objective function `f`.
#'   The components of the list may be either objects of `typevar` class (see
#'   [type_variable]) or strings identifying one of [type_variable] functions
#'   (for example `"real"` for function [type_real()]). See § Examples.
#' @param ... other arguments to be passed to `f` or to other methods (for
#'   `print` and `plot`).
#' @param rho parameter \eqn{\rho} of the Cross-Entropy algorithm. This argument
#'   may be passed either as a numeric value in \eqn{(0,1)} or as an unevaluated
#'   expression which may include the number of current iteration `n`, or the
#'   argument `N`.
#' @param N parameter \eqn{N} of the Cross-Entropy algorithm. This argument
#'   may be passed either as a positive integer or as an unevaluated expression
#'   which may include the number of current iteration `n`.
#' @param smooth list of unevaluated expressions to be used as smoothing rules
#'   for the parameters of the sampling probability distributions of **all
#'   variables**. If not `NULL`, all default or set smoothing rules of all
#'   variables will be overwritten. See [type_variable] for details and examples.
#' @param stopwindow unevaluated expression returning the object to be passed to
#'   the stopping rule. Symbol `gam` permits the time series \eqn{\gamma_t} to
#'   be used (as a `numeric` vector).
#' @param stoprule stopping rule passed as an unevaluated expression including
#'   `x` as the object returned by evaluation of argument `stopwindow`. The
#'   algorithm is stopped when zero is returned by the evaluation of `stoprule`.
#'   If returned object has attribute `mess`, this is used as a message.
#'   Currently, built-in stopping rules are [ts_change()] and [geweke()], others
#'   may be defined by user.
#' @param maxiter maximum number of iteration. When it is reached, algorithm is
#'   stopped whether or not the stopping criterion is satisfied. If the maximum
#'   number of iteration is reached, the `code` and the `message` components of
#'   `noisyCE` object are overwritten.
#' @param maximise if `TRUE` (default) `f` is maximised, otherwise a
#'   minimisation of `f` is performed.
#' @param x,object object of class `noisyCE2`, as returned by `noisyCE2`.
#' @param verbose algorithm verbosity (values `v`, `vv` and `vvv` are admitted).
#' @param what type of plot should be drawn. If `what = "x"` (default), values
#'   of the variables are plotted as time series; if `what = "gam"`, time series
#'   of statistics \eqn{\gamma} is plotted; if `what = "param"`, time series of
#'   parameters of the sampling distributions are plotted.
#' @param start,end first and last value to be plotted. If `NULL`, all values
#'   are plotted.
#'
#' @return
#' An object of class `noisyCE2` structured as a list with the following
#' components:
#' \item{f}{argument `f`.}
#' \item{fobj}{objective function `f` where possible arguments passed through
#'   argument `...` have been substituted. Thus, the value of the objective
#'   function maximised by `noisyCE` in `x0` can be computed as `fobj(x0)`. If
#'   a minimisation has been performed, `fobj` returns `f` with sign inverted.}
#' \item{xopt}{`numeric` vector with solution.}
#' \item{hxopt}{matrix of `niter` rows and `length(xopt)` columns with values of
#'   variables generated by the optimisation algorithm.}
#' \item{param}{`list` of `length(xopt)` components where time series of
#'   parameters (vectors \eqn{v_t}) are stored for each variable as `data.frame`
#'   objects with `niter+1` rows (the first rows are the starting values set
#'   through function `noisyCEcontrol`).}
#' \item{gam}{vector of values \eqn{\gamma_t}.}
#' \item{niter}{number of iterations.}
#' \item{code}{convergence code of the algorithm. Value `0` means that algorithm
#'   has converged; other values are defined according to the stopping rule.}
#' \item{convMess}{textual message associated to the convergence code (if any).}
#' \item{compTimes}{named vector computation times of each phase.}
#'
#' @examples
#' library(magrittr)
#' # Optimisation of the 4-dimensional function:
#' # f(x1,x2,x3,x4)=-(x1-1)^2-(x2-2)^2-(x3-3)^2-(x4-4)^2
#' sol <- noisyCE2(function(x) -sum((x - (1:4))^2), domain = rep('real', 4))
#' # Representation of the convergence process:
#' plot(sol, what = 'x')
#' plot(sol, what = 'gam')
#'
#' @export
noisyCE2 <- function(f, domain, ..., rho = 0.05, N = 1000,
  smooth = NULL,
  stopwindow = tail(gam, (n > 20) * n / 2), stoprule = ts_change(x),
  maxiter = 1000, maximise = TRUE, verbose = 'v') {
  
  # Read date and time
  compTimes <- Sys.time()
  
  # Redefine the objective function including the other parameters
  fobj <- function(x) { (2 * maximise - 1) * f(x, ...) }
  
  # Save variable names
  namesX <- names(domain)
  if (is.null(namesX)) {
  	namesX <- paste0('X', seq_along(domain))
  }
  
  # Define the block of variables
  if (is.null(smooth)) {
  	block <- generate_block(domain)
  } else {
  	block <- generate_block(domain, smooth = smooth)
  }
  
  # Initialisations
  n <- 0
  vt <- lapply(block, function(x) { x$init })
  vhist <- lapply(vt, function(x) as.data.frame(matrix(x, 1)))
  gam <- NULL
  convCode <- -2
  
  # Cross-entropy algorithm
  compTimes %<>% c(Sys.time())
  
  while((convCode != 0) & (n < maxiter)) {
  	
    # Update the counter
    n %<>% add(1)

    # Simulate the N values of the independent variables
    X <- mapply(
      FUN = function(varj, theta) { varj$randomXj(eval(N), theta) },
      block, vt, SIMPLIFY = FALSE
    )

    # Simulate the N values of the objective function
    X %<>% Reduce(cbind, ., NULL)
    fx <- apply(X, 1, fobj)
    
    # Compute and store the new gamma
    fx %>%
      quantile(1 - eval(rho), type = 1, na.rm = TRUE) %>%
      c(gam, .) -> gam
    
    # Update the parameters
    X[fx >= gam[n], ] %>%
      data.frame %>%
      mapply(
        FUN = function(varj, xj) { varj$x2v(xj) },
        block, ., SIMPLIFY = FALSE
      ) -> vt

    # Smooth the parameters
    xt <- NULL # Only set for avoiding NOTE in package check
    vt %<>% smooth_block(block, ., vhist)
    vhist %<>% concat_block_vect(vt)
    
    # Progress message
    vmessage(verbose, 1,
      '\r     * iter. ',
      format(n, digits = 4, nsmall = 0, width = 4),
      ' - gamma: ',
      format(gam[n], digits = 4, nsmall = 4),
      '       ',
      sep = '', appendLF = FALSE
    )
    
    # Check the stopping criterion
    x <- eval(substitute(stopwindow))
    if (length(x) > 0) {
      convCode <- unname(eval(substitute(stoprule)))
    }
  }
  
  compTimes %<>% c(Sys.time())
  
  if((convCode != 0) & (n >= maxiter)) {
  	convCode <- -1
  	attr(convCode, 'convMess') <- 'the maximum number of iterations has been reached'
  }
  
  # Output
  mess <- attr(convCode, 'convMess')
  vmessage(verbose, 1, '\n     * The algorithm has ', mess)

  #  
  mapply(
    FUN = function(varj, vj) { unlist(varj$v2x(vj)) },
    block, vhist, SIMPLIFY = FALSE
  ) %>%
    Reduce(cbind, .) %>%
    data.frame %>%
    set_rownames(NULL) -> hxopt

  if(any(is.null(namesX))) {
    colnames(hxopt) <- paste0('X', 1:ncol(hxopt))	
  } else {
  	colnames(hxopt) <- namesX
  }
  
  tail(hxopt, 1) %>%
    as.numeric %>%
    set_names(colnames(hxopt)) -> xopt
  
  vhist %<>% set_names(namesX)
  
  #
  compTimes %<>%
    c(Sys.time()) %>%
    diff %>%
    set_names(c('Initialisation', 'Optimisation', 'Preparing output'))
  
  # Output
  list(
    f = f, fobj = fobj, xopt = xopt, hxopt = hxopt, param = vhist,
    gam = unname(gam), niter = n, code = convCode, convMess = mess,
    compTimes = compTimes
  ) %>%
    structure(class = 'noisyCE2') %>%
    return()
}



#' @describeIn noisyCE2 display synthetic information about a `noisyCE2` object
#' @export
print.noisyCE2 <- function(x, ...) {
  cat('Object of class "noisyCE2"\n')
  cat('--------------------------------------\n')
  cat('Status               :', x$convMess, '\n')
  cat('Computation time     :', 
    format(unclass(sum(x$compTime)), digits = 4),
    attr(x$compTime, 'units'), '\n'
  )
  cat('Number of iterations :', x$niter, '\n')
  cat('Time per iteration   :',
    format(unclass(sum(x$compTime) / x$niter), digits = 4),
    attr(x$compTime / x$niter, 'units'), '\n'
  )
  cat('Last gamma value     :', tail(x$gam, 1), '\n')
  cat('Solution x*          :', x$xopt, '\n')
  #cat('Optimal value f(x*)  :', x$fobj(x$xopt), '\n')
  invisible(x)
}



#' @describeIn noisyCE2 display summary information about a `noisyCE2` object
#' @export
summary.noisyCE2 <- function(object, ...) {
  print(object, ...)
}



#' @describeIn noisyCE2 plot various components of a `noisyCE2` object
#' @export
plot.noisyCE2 <- function(x, what = c('x', 'gam', 'param'),
  start = NULL, end = NULL, ...) {
  	
  what %<>% match.arg(c('x', 'gam', 'param'))
  
  if (what == 'x') {
    stats::ts(x$hxopt, 1) %>%
      stats::window(start, end) %>%
      stats::plot.ts(...)
  }
  if (what == 'gam') {
    stats::ts(x$gam, 1) %>%
      stats::window(start, end) %>%
      stats::plot.ts(...)
  }
  if (what == 'param') {
    for(j in seq_along(x$param)) {
      colnames(x$param[[j]]) %<>% paste0(names(x$param)[j], '_', .)
    }
    x$param %>%
      Reduce(cbind, .) %>%
      stats::ts(0) %>%
      stats::window(start, end) %>%
      stats::plot.ts(...)
  }
  
  invisible(NULL)
}



#' @describeIn noisyCE2 get the solution of the optimisation
#' @export
coef.noisyCE2 <- function(object, ...) {
  object$xopt
}