R/OptimizeModels.R

In fcampelo/MetaTuner: MetaTuner: Metaheuristics Tuner

#' Optimize regression models to estimate new configurations
#'
#' @param parameters data frame containing information about the tunable
#'                   parameters
#' @param models list containing models generated by [FitModels()].
#' @param optimization.method optimization method to be used
#' @param ndigits number of significant digits to use for each parameter.
#'
#' @return list containing new candidate configurations found by optimization
#'
#' @author Felipe Campelo (\email{fcampelo@@ufmg.br}),
#'         Athila Trindade (\email{rochaathila@@gmail.com})
#'

OptimizeModels <- function(parameters,
                           models,
                           optimization.method = "Nelder-Mead",
                           ndigits){

  ## ==============
  ## Error checking done in the calling routine
  ## ==============

  # ========== Prepare matrix of new configurations
  mysample <- matrix(as.numeric(NA),
                     nrow = length(models),
                     ncol = nrow(parameters))


  # ========== Prepare optimization parameters
  # Objective function
  myobjfun <- function(x, mymodel, mypars){
    newX        <- as.data.frame(matrix(rep(x, 2), nrow  = 2, byrow = TRUE))
    names(newX) <- mypars$name
    modclass    <- class(mymodel$model)

    if (modclass == "lm"){
      y <- stats::predict.lm(object  = mymodel$model,
                             newdata = newX)[1]
    } else if (modclass == "rq"){
      y <- quantreg::predict.rq(object  = mymodel$model,
                                newdata = newX)[1]
    } else if (modclass == "hqreg"){
      #y <- hqreg::predict.hqreg(object = mymodel$model,
      #                          X      = x,
      #                          lambda = mymodel$model$lambda[2])[1]

      coefs         <- mymodel$model$beta[,2]
      expcoefs      <- names(coefs)
      listexponents <- mapply(function(exps, pattern){
        stringi::stri_split_fixed(exps, pattern)},
        exps = expcoefs[2:length(expcoefs)], pattern = ".")

      Xaux          <- rep(1, length(listexponents))

      for(i in seq(listexponents))
        for (j in 1:length(x)) {
          Xaux[i] <- Xaux[i] * (x[j] ^ as.numeric(listexponents[[i]][[j]]))
        }

      ###This last two lines I have to do since the function predict didn't work
      # with an hqreg object when running it in parallel (I have no clue why
      # it happened)
      y <- coefs[1]          ##Intercept
      for(i in seq(listexponents)) y <- y + Xaux[i] * coefs[i + 1]
    }
    else {
      stop("Model class", modclass, "not recognized by function OptimizeModels")
    }

    return(sign(y) * min(1e32, abs(y))) # <-- prevents infinity
  }

  # Box constraints
  ci <- rep(c(0, -1), each = nrow(parameters))
  Ui <- diag(nrow = nrow(parameters),
             ncol = nrow(parameters))
  Ui <- rbind(Ui, -Ui)

  #========== Optimize (minimize) models in a parallel environment
  i <- NULL
  optparams <- foreach::foreach(i = models, .combine = 'rbind') %dopar% {
    theta <- stats::runif(nrow(parameters))
    Y <- stats::constrOptim(theta    = theta,
                            f        = myobjfun,
                            grad     = NULL,
                            ui       = Ui,
                            ci       = ci,
                            method   = optimization.method,
                            mymodel  = i,
                            mypars   = parameters)
    return(Y$par)
  }

  mysample <- optparams
  rownames(mysample) <- NULL

  for (j in 1:ncol(mysample)){
    mysample[, j] <- round(mysample[, j], digits = ndigits[j])
  }

  # # PARALLEL-IZE HERE
  # # VVVVVVVVVVVVVVVV
  #   for (i in seq(models)){
  #   # Initial point for optimization (random, feasible)
  #   theta <- runif(nrow(parameters))
  #   # cat("\n", theta)
  #
  #   Y <- stats::constrOptim(theta      = theta,
  #                           f          = myobjfun,
  #                           grad       = NULL,
  #                           ui         = Ui,
  #                           ci         = ci,
  #                           method     = optimization.method,
  #                           mymodel    = models[[i]],
  #                           mypars     = parameters)
  #   mysample[i, ] <- Y$par
  #   cat(".")
  # }
  #
  #
  # for (j in 1:ncol(mysample)){
  #   mysample[, j] <- round(mysample[, j], digits = ndigits[j])
  # }

  # ========== Return new configurations
  newconfs <- apply(X      = mysample,
                    MARGIN = 1,
                    FUN    = function(x){
                      list(config = x,
                           Yij    = data.frame(instance.ID = character(),
                                               y           = numeric(),
                                               stringsAsFactors = FALSE),
                           perf   = NA)})
  return(newconfs)

}