R/maximin_mixexp.R

In AHM: Additive Heredity Model: Method for the Mixture-of-Mixtures Experiments

#' The mymaximin function generates the matrix of maximin design points.
#' It uses the simplex centroid design as the base design, then in a stochastics way sample the candidate design points generated by the function partition.
#'
#' @param pool, partition the base design points provided to the function
#' @param n numeric, sample size
#' @param m numeric, 3 stands for 3 components, i.e. c1, c2, and c3
#' @param iter numeric, iterations used in the stochastic sampling
#' @param Xorig matrix, initial design points
#' @description {This method is modified based on Prof. Bobby Gramacy's Computer Experiment lecture at Virginia Tech. \href{http://bobby.gramacy.com/teaching/}{Prof. Gramacy's lecture website}}
#' @return Return a matrixt of the design points for the major components
#' @export
#' @examples
#' \donttest{
#' # The case of unconstrainted experiments
#' library(mixexp)
#' num_size = 8 # num points in the design for the major component
#' Xorig = as.matrix(SCD(3))
#' # all possible combinations sum to 1
#' pool_3d =partitions::compositions(1000, 3,include.zero = TRUE)/1000
#' res_C = mymaximin(pool=pool_3d, n=num_size, m=3, iter=1e5, Xorig=Xorig)
#' DesignPoints(res_C,cornerlabs = c("c3","c2","c1"),axislabs=c("c1","c2","c3"))
#'
#' # The case of constrainted experiments
#' library(mixexp)
#' num_size = 8 # num points in the design for the major component
#' # all possible combinations sum to 1
#' pool_3d =partitions::compositions(1000, 3,include.zero = TRUE)/1000
#' c1_min=0.2
#' c1_max=0.45
#' c2_min=0.4
#' c2_max=0.6
#' c3_min=0.1
#' c3_max=0.25
#' tmp = Xvert(nfac=3,lc=c(c1_min,c2_min,c3_min),uc =c(c1_max,c2_max,c3_max),ndm=1,pseudo=FALSE)
#' Xorig=tmp[c(1:6,13),c(1:3)]
#' colnames(Xorig)=c("V1","V2","V3")
#' pool_3d = t(dplyr::filter(as.data.frame(t(as.matrix(pool_3d))),t(pool_3d)[,1] > c1_min &
#'                      t(pool_3d)[,1] <= c1_max &
#'                      t(pool_3d)[,2] > c2_min &
#'                      t(pool_3d)[,2] <= c2_max &
#'                      t(pool_3d)[,3] > c3_min &
#'                      t(pool_3d)[,3] <= c3_max
#' )
#' )
#' res_C = mymaximin(pool=pool_3d, n=num_size, m=3, iter=1e5, Xorig=Xorig)
#' DesignPoints(res_C,cornerlabs = c("c3","c2","c1"),axislabs=c("c1","c2","c3")
#'                      ,x1lower=c1_min,x2lower=c2_min,x3lower=c3_min
#'                      ,x1upper=c1_max,x2upper=c2_max,x3upper=c3_max, pseudo=FALSE)
#' }


mymaximin <- function(pool, n=50, m=3, iter=1e5, Xorig=NULL) {
  # n = 3 for c1, c2, c3
  # m is the sample size
  # pool is from the partition results
  # library(plgp, quietly=TRUE)
  # pool = pool_3d

  # X <- matrix(runif(n*m), ncol=m)  ## initial design
  if (is.null(Xorig)) {
    # use Xorig to include certain points
    n_size = n
  } else {
      n_size = n - nrow(Xorig)
    }

  init_id = sample(c(1:dim(pool)[2]),size = n_size,replace = FALSE)
  X = t(pool[,init_id])
  X <- rbind(X, Xorig)             # combine initial base and randomly picked row
  d <- plgp::distance(X) #distance(X) #
  d <- as.numeric(d[upper.tri(d)])
  md <- min(d)

  for(t in 1:iter) {
    #print(paste("iter is ", t, sep=''))
    row <- sample(1:n, 1)
    xold <- X[row,]       ## random row selection
    X[row,] <- pool[,sample(c(1:dim(pool)[2]),size = 1,replace = FALSE)] ## runif(m)   ## random new row
    dprime <- plgp::distance(X) #distance(X)  #
    dprime <- as.numeric(dprime[upper.tri(dprime)])
    mdprime <- min(dprime)
    if(mdprime > md) { md <- mdprime  ## accept
    } else { X[row,] <- xold }        ## reject
  }

  rownames(X) = c(1:nrow(X))
  return(X) # design matrix
}