R/Search.R

In MOODE: Multi-Objective Optimal Design of Experiments

#' Searching for a multi-objective optimal completely randomised design.
#' @description Performing search for a (nearly) optimum factorial design, optimising with respect to a specified compound criterion.
#' 
#' @param mood.object Object of class `mood`, generated by `mood` function, containing the parameters of the experiment, the compound criterion and search parameters 
#' @param algorithm Parameter specifying the search algorithm. If `ptex` (default for \eqn{K<=4}), the point-exchange algorithm is used, and if `coordex` (default for \eqn{K>4}), the coordinate-exchange.
#' @param parallel If `TRUE` use the `doFuture` package to run independent iterations of the algorithm in parallel using `foreach`. Requires `doFuture` library to be installed and a `Future` `plan` to be specified. See examples. 
#' @param verbose If `TRUE`, progress messages through the search iterations are shown.
#' @export
#' @details `Search` takes the mood object as an input with all the parameters of the experiment. Runs a point-exchange or a coordinate-exchange algorithm, returns design and model matrices, computation time and criteria values.
#' See \insertCite{KoutraMOODE}{MOODE} for examples of using `parallel = TRUE`.
#' 
#' @return Updated object of class `mood` containing the outputs generated by the search:
#' \itemize{
#' \item `X.design` Design matrix.
#' \item `df` The number of pure error degrees of freedom.
#' \item `X1` Primary model matrix for the found (nearly-) optimum design.
#' \item `X2` Model matrix of potential terms for the found (nearly-) optimum design.  
#' \item `compound.value` The compound criterion value of the (nearly-) optimum design.
#' \item `criteria.values` Component criteria values of the (nearly-) optimum design.
#' \item `path` The "path" of compound criterion values of the optimum designs obtained after for each random start.
#' \item `time` Computation time.
#' \item `algorithm` Point exchange or coordinate exchange used to find the design?
#' \item `parallel` Were different runs of the algorithm performed across different CPU cores (`TRUE`/`FALSE`) 
#' }
#' 
#' @seealso [mood]
#' @references 
#'    \insertAllCited
#' @examples
#' 
#'example1 <- mood(K = 2, Levels = 3, Nruns = 10, criterion.choice = "GDP", 
#'                 kappa = list(kappa.Ds = 1./3, kappa.DP = 1./3, kappa.LoF = 1./3), 
#'                 control = list(tau2 = 0.1), 
#'                 model_terms = list(primary.model = "first_order", 
#'                 potential.terms = c("x12", "x22", "x1x2")))
#' # Using point exchange
#'Search_point <- Search(example1, algorithm = 'ptex')
#'Search_point
#' # Using coordinate exchange (the default for K>4)
#'Search_coord <- Search(example1)
#'Search_coord


Search <- function(mood.object, algorithm = c("ptex", "coordex"), parallel = FALSE, verbose = TRUE)
{
  if (!inherits(mood.object, "mood")) cli::cli_abort(c("First argument must be an object of class 'mood'"))
  
  K <- mood.object$K
  if(K > 4 && length(algorithm) > 1 && !mood.object$orth) {
    algorithm <- "coordex"
    cli::cli_warn(c("Changed algorithm to coordinate exchange", 
                   "!" = "For K > 4, algorithm is set to coordinate exchange (\"coordex\") for performance."), 
                 immediate. = TRUE)
  } else {
    algorithm <- match.arg(algorithm)
  }
  
  
  if(identical(algorithm, "ptex")) {
    pointex <- TRUE
  } else {
    pointex <- FALSE
  }
  Levels <- mood.object$Levels
  Klev <- mood.object$Klev
  
  P <- mood.object$P
  Q <- mood.object$Q
  primary.terms <- mood.object$primary.terms
  potential.terms <- mood.object$potential.terms
  
  Nruns <- mood.object$Nruns
  Nstarts <- mood.object$Nstarts
  Biter <- mood.object$Biter
  
  W <- mood.object$W; Z0 <- mood.object$Z0
  tau2 <- mood.object$tau2; tau <- mood.object$tau
  
  criterion.choice <- mood.object$criterion.choice
  Cubic <- mood.object$Cubic
  orth<- mood.object$orth
  
  kappa.L <- mood.object$kappa.L
  kappa.LP <- mood.object$kappa.LP
  kappa.Ds <- mood.object$kappa.Ds
  kappa.DP <- mood.object$kappa.DP
  kappa.LoF <- mood.object$kappa.LoF
  kappa.bias <- mood.object$kappa.bias
  kappa.mse <- mood.object$kappa.mse
  
  alpha.DP <- mood.object$alpha.DP
  alpha.LP <- mood.object$alpha.LP
  alpha.LoF <- mood.object$alpha.LoF
  alpha.LoFL <- mood.object$alpha.LoFL
  
  search.object <- list("Nruns" = Nruns, "P" = P, "Q" = Q,
                        "criterion.choice" = criterion.choice,
                        "primary.terms" = primary.terms,
                        "potential.terms" = potential.terms,
                        "Biter" = Biter, "tau2" = tau2,"tau" = tau, 
                        "W" = W, "Z0" = Z0,
                        "alpha.DP" = alpha.DP,"alpha.LP" = alpha.LP,
                        "alpha.LoF" = alpha.LoF, "alpha.LoFL" = alpha.LoFL,
                        "kappa.Ds" = kappa.Ds, "kappa.DP" = kappa.DP,
                        "kappa.L" = kappa.L, "kappa.LP" = kappa.LP,
                        "kappa.LoF" = kappa.LoF, "kappa.bias" = kappa.bias, 
                        "kappa.mse" = kappa.mse)
  
  start_time <- Sys.time()
  
  cand.trt <- candidate_trt_set(Levels, K, Cubic)   # form the candidate set of treatments
  cand.full <- candidate_set_full(cand.trt, K)      # build candidate set, with potential terms
  
  if (orth) {
    cand.full <- candidate_set_orth(cand.full,
                                    primary.terms, potential.terms)
    cli::cli_alert_info("Point exchange algorithm will be used for othonormalised candidate sets")
    pointex <- TRUE
  }

  ###### if coordinate exchange, need to change labels as well
  if (!pointex){
    # list of factors' levels scaled to [-1,1]
    levels_scaled <- lapply(Levels, Transform) 
    levels_steps <- rep(1, K)  
    
    if(K > 1) {
      for (i in (K-1):1)
      { # how much label changes per one level change of each of the factors
        # based on how the labelling is organised in the label() function
        levels_steps[i] <- levels_steps[i+1]*length(Levels[[i+1]])
      }
    }
    steps = 2./(sapply(Levels, length) - 1) # how long is one step for each factor
    
    search.object$levels_scaled <- levels_scaled # list of levels, scaled to [-1,1]
    search.object$levels_steps <- levels_steps   # difference in labels when levels moved to the next
    search.object$steps <- steps
    search.object$levels_lengths <- sapply(Levels, length)
  }
  
  kx <- 1:Nstarts
  progressr::with_progress({
    p <- progressr::progressor(along = kx)
  
    if(isTRUE(parallel) & !requireNamespace("doFuture", quietly = TRUE)) {
      warning(
        "Package \"doFuture\" required for parallel computing, running in squential mode."
        )
      parallel <- FALSE
    }
    
    if(isTRUE(parallel)) {
      `%dofuture%` <- doFuture::`%dofuture%`
      designs <- foreach::foreach(k = kx, .options.future = list(seed = TRUE)) %dofuture% {
        if(isTRUE(verbose)) p(message = sprintf("Current iteration: %i out of %i", k, Nstarts))
      
        initial <- initial.full(cand.full, Nruns, primary.terms, potential.terms)
        X1 <- initial$X1
        X2 <- initial$X2
    
        if (k==1){
          crit.opt <- objfun(X1, X2, search.object)$compound
          X1.opt <- X1; X2.opt <- X2
        }
        s <- 1
        while (s==1)
        {
          if(pointex) {
            Xs <- point.swap(X1, X2, cand.full, search.object)
          } else {
            Xs <- coord.swap(X1, X2, K, Levels, search.object)
          }
          X1 <- Xs$X1; X2 <- Xs$X2;
          s <- Xs$search
        }
        list(X1 = X1, X2 = X2, value = Xs$compound)
      }
      crit_values <- sapply(designs, \(x) x$value)
      final <- designs[[which.min(crit_values)]]
      X1.opt <- final$X1
      X2.opt <- final$X2
      crit.opt <- final$value
    } else {
      crit_values <- rep(0, Nstarts)
      for (k in 1:Nstarts){
        if(isTRUE(verbose)) p(message = sprintf("Current iteration: %i out of %i", k, Nstarts))
        
        initial <- initial.full(cand.full, Nruns, primary.terms, potential.terms)
        X1 <- initial$X1
        X2 <- initial$X2
        
        if (k==1){
          crit.opt <- objfun(X1, X2, search.object)$compound
          X1.opt <- X1; X2.opt <- X2
        }
        s <- 1
        while (s==1)
        {
          if(pointex) {
            Xs <- point.swap(X1, X2, cand.full, search.object)
          } else {
            Xs <- coord.swap(X1, X2, K, Levels, search.object)
          }
          X1 <- Xs$X1; X2 <- Xs$X2;
          s <- Xs$search
        }
        crit_values[k] <- Xs$compound # track the change of criterion values
        if (crit.opt > Xs$compound)
        {
          X1.opt <- Xs$X1; X2.opt<-Xs$X2
          crit.opt <- Xs$compound
        }
      }
    }
  })
  
  finish_time <- Sys.time()
  time <- finish_time - start_time

  criteria.opt <- objfun(X1.opt, X2.opt, search.object)
  X.design <- X1.opt[,3:(K+2)]
  
  cli::cli_alert_success("Design search complete. Final compound objective function value = {round(criteria.opt$compound, digits = 5)}")
  
  out <- list (X.design = X.design, df = criteria.opt$df, 
        X1 = X1.opt, X2 = X2.opt, 
        compound.value = criteria.opt$compound, 
        criteria.values = criteria.opt, path = crit_values,
        time = time, algorithm = algorithm, parallel = parallel)
  mood.object$searched <- TRUE
  structure(c(mood.object, out), class = "mood")
}