R/maximin.R

In maximin: Space-Filling Design under Maximin Distance

#*******************************************************************************
#
# Space-filling Design under Maximin Distance
# Copyright (C) 2018, Virginia Tech
#
# This library is a free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License as published by the Free Software Foundation; either
# version 2.1 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with this library; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
#
# Questions? Contact Furong Sun (furongs@vt.edu) and Robert B. Gramacy (rbg@vt.edu)
#
#*******************************************************************************

## maximin:
##
## generates a space-filling design in a CONTINUOUS but BOUNDED DESIGN SPACE under the criterion of maximum-minimum distance
# Xorig is for sequential maximin design;
# Two questions to be explored:
##    1. monitor convergence to avoid too big T
##    2. perform simulated annealing to avoid local maxima
maximin <- function(n, p, T=10*n, Xorig=NULL, Xinit=NULL, verb=FALSE, plot=FALSE, boundary=FALSE){ 
  
  ####################################### sanity checks #######################################
  if(!is.null(Xorig)){
     # if(class(Xorig) != "matrix"){
     #    Xorig <- as.matrix(Xorig)
     # }else 
    if(ncol(Xorig) != p){
        stop("column dimension mismatch between X and Xorig") 
     }
  }
  
  if(n <= 1) stop("n must be bigger than 1.")
  
  if(T <= n) warning("T should be bigger than n.")
  
  ## the initial design: dimensionality of n * p
  if(!is.null(Xinit)){
     X <- Xinit ## from a previous experiment, OR
  }else{
     X <- matrix(runif(n*p), ncol=p) ## randomly generated within the global input space
  }
  
  ############# Below are the generic functions required for this algorithm #############
  ## a generic function to calculate the distance between the "to-be-swapped-in" location and the boundary
  d.bound <- function(xnew){
    
    xnew <- as.numeric(xnew)
    ## the distance to one boundary: should be multiplied by 4, but I want to do it later
    d.k1 <- xnew^2
    ## the distance to the other boundary: ...
    d.k2 <- (1-xnew)^2
    return(4*c(d.k1, d.k2)) ## return four times of both: we will 
                            ## calculate the 'md' inside the function, maximin.obj
  } # Why FOUR times?
  
  ## the objective function to be optimized (maximized)
  maximin.obj <- function(X, row.out, xrow.in, Xorig=NULL){
    
    ## the distance between the "to-be-swapped-in" location, and other X locations, the existing design, and the boundary
    d <- as.numeric(distance(matrix(xrow.in, nrow=1), rbind(X[-row.out,], Xorig)))
    
    if(boundary==TRUE){
       d <- c(d, d.bound(xrow.in))
    }
    return(sqrt(min(d))) # use the square root of the md as the objective function since it provides much less extreme values
  }
  
  ## the search space when optimizing: xd is the nearest location of xs
  sp <- function(xs, xd=NULL, md){
    
    if(!is.null(xd)){
       ds <- as.numeric(distance(matrix(xs, nrow=1), matrix(xd, nrow=1)))
    }else{
       ds <- md
    }
    
    ## treat 'xs' as the center and 'sqrt(ds)' as the radius to justify a hypercube
    lo <- xs - sqrt(ds) # lower bound
    up <- xs + sqrt(ds) # upper bound
    
    ## guarantee the search space to be within the global input space: [0,1]^p
    lo[lo<0] <- 0; up[up>1] <- 1
    
    return(list(lo=lo, up=up, ds=ds)) ## save the distance between 'xs' and 'xd', ds, for future use
  }
  
  ## The following function is to be used when "jumping out of the well":
  ##    to justify the start location, which can not be the same as 'xs', 
  ##                  however, it should be close to 'xs' as much as possible;
  ##    'xd' is the closest X or Xorig location to 'xs';
  ##    The start location will be in-between 'xs' and 'xd' and its distance 
  ##    to 'xs' is 5% of 'sqrt(ds)': 'sqrt(ds)' is the Euclidean distance between 'xs' and 'xd':
  st <- function(xs, xd, ds){
     
     st.prime <- rep(NA, length(xs))
     for(i in 1:length(xs)){
         
         if(xs[i] == xd[i]){
            st.prime[i] <- xs[i]
         }else{
            st.prime[i] <- ifelse(xs[i] > xd[i], 
                                  xs[i] - 0.05*sqrt(ds), ## not sure it is right
                                  xs[i] + 0.05*sqrt(ds))
         }
     }
     return(st.prime)
  }
  
  ## to check whether each element of a numeric vector is the same
  cvs <- function(x){
         diff(range(x)) < sqrt(.Machine$double.eps)
  }
  
  ########### done with generic functions, the algorithm formally starts ###########
  D <- plgp::distance(X)
  md <- min(as.numeric(D[upper.tri(D)]))
  md.ind <- which(D==md, arr.ind=TRUE)[,1] ## the list of runs among X with `md'
  
  if(!is.null(Xorig)){
     D2 <- plgp::distance(X, Xorig)  ## an nrow(X) * nrow(Xorig) matrix
     md2 <- min(D2)                  ## the 'md' between X and Xorig
     if(md2 < md){
        md <- md2
        md.ind <- which(D2==md2, arr.ind=TRUE)[,1]
     }else if(md2 == md){
        md.ind <- unique(c(md.ind, which(D2==md2, arr.ind=TRUE)[,1]))
     }
  }
  
  ## preallocate the 'md' and 'design matrix' for each iteration
  mind <- rep(NA, T+1); mind[1] <- md
  Xind <- vector("list", T+1); Xind[[1]] <- X
  
  ## build up the design iteratively
  for(t in 1:T){
    
    ## track the locations with 'non-md'
    cand.ind <- setdiff(1:n, md.ind) 
    
    ## swap out the location with 'md'
    row.out.ind <- ceiling(runif(1)*length(md.ind))
    row.out <- md.ind[row.out.ind]
    xold <- X[row.out,]
    
    ## define the search space
    sw.in <- sp(xs=xold, md=md)
    low <- sw.in$lo  ## lower bound
    upp <- sw.in$up  ## upper bound
    
    ## the start location for 'optim' must be in the class of matrix
    start <- matrix(xold, nrow=1)
    ## try to maximize the 'md' between the "to-be-swapped-in" location and other X runs & Xorig, constrained by the search space
    out <- optim(par=start, fn=maximin.obj, method="L-BFGS-B", 
                 lower=low, upper=upp, control=list(fnscale=-1), # to maximize, instead of minimizing
                 X=X, row.out=row.out, Xorig=Xorig)
    md.ind <- md.ind[!md.ind==row.out] ## update the 'md' list
    
    if(plot == TRUE){
       d <- sample(1:p, 2, replace=FALSE) ## randomly choose TWO coordinates from the input space to plot
       Xv <- X[,d]
       plot(Xv, xlab=paste0("x", d[1]), ylab=paste0("x", d[2]), 
            xlim=c(0,1), ylim=c(0,1), main=paste0("md.loc at it=", t))
       points(Xorig[,d], pch=20, col="forestgreen")
       points(xold[d[1]], xold[d[2]], pch=20, col="red")
       rect(low[d[1]], low[d[2]], upp[d[1]], upp[d[2]], col=rgb(0,1,0,alpha=0.2)) ## the search space
       if(!is.null(Xorig)){
          legend("topright", c("xold", "Xorig"), xpd=TRUE, 
                 horiz=TRUE, inset=c(-0.015, -0.045), pch=20, col=c("red","forestgreen"), bty='n')
       }else{
          legend("topright", "xold", xpd=TRUE, 
                 horiz=TRUE, inset=c(-0.015, -0.045), pch=20, col="red", bty='n')
       }
    }
    
    while(all.equal(as.numeric(start), as.numeric(out$par))==TRUE){ ## use `all.equal` instead of `identical` due to numerical precision: `identical' is more numerically precise than `all.equal`
      if(length(md.ind) >= 1){ ## just in case there are more than one location with 'md' ...    ## the 1st strategy
         row.out.ind <- ceiling(runif(1)*length(md.ind))
         row.out <- md.ind[row.out.ind]
         xold <- X[row.out,]
         
         sw.in <- sp(xs=xold, md=md)
         low <- sw.in$lo
         upp <- sw.in$up
         
         start <- matrix(xold, nrow=1)
         out <- optim(par=start, fn=maximin.obj, method="L-BFGS-B",
                      lower=low, upper=upp, control=list(fnscale=-1),
                      X=X, row.out=row.out, Xorig=Xorig)
         md.ind <- md.ind[!md.ind==row.out] ## update the 'md' list
         
         if(plot==TRUE){
            Xv <- X[,d]
            plot(Xv, xlab=paste0("x", d[1]), ylab=paste0("x", d[2]), 
                 xlim=c(0,1), ylim=c(0,1), main=paste0("md.loc at it=", t))
            points(Xorig[,d], pch=20, col="forestgreen")
            points(xold[d[1]], xold[d[2]], pch=20, col="red")
            rect(low[d[1]], low[d[2]], upp[d[1]], upp[d[2]], col=rgb(0,1,0,alpha=0.2)) # the search space
            if(!is.null(Xorig)){
               legend("topright", c("xold", "Xorig"), 
                      xpd=TRUE, horiz=TRUE, inset=c(-0.015, -0.045), pch=20, col=c("red","forestgreen"), bty='n')
            }else{
               legend("topright", "xold", 
                      xpd=TRUE, horiz=TRUE, inset=c(-0.015, -0.045), pch=20, col="red", bty='n')
            }
        }
      }else if(length(cand.ind) >= 1){ ## if there is no run with 'md', then randomly start with a run with 'non-md' ... ## the 2nd strategy
        row.out.ind <- ceiling(runif(1)*length(cand.ind))
        row.out <- cand.ind[row.out.ind]
        xold <- X[row.out,]
        
        ## combine the distance matrices
        if(is.null(Xorig)==TRUE) D2 <- NULL
        DD2 <- cbind(D, D2)
        
        ## justify 'xd'
        Dor <- order(DD2[row.out,], decreasing=FALSE)
        xd <- ifelse(rep(Dor[2]<=n, p), X[Dor[2],], Xorig[Dor[2]-n,])
        
        ## define the search space
        sw.in <- sp(xs=xold, xd=xd, md=NULL)
        low <- sw.in$lo
        upp <- sw.in$up
        
        start <- matrix(xold, nrow=1)
        out <- optim(par=start, fn=maximin.obj, method="L-BFGS-B",
                     lower=low, upper=upp, control=list(fnscale=-1),
                     X=X, row.out=row.out, Xorig=Xorig)
        cand.ind <- cand.ind[!cand.ind==row.out] ## update the 'non-md' list
        
        if(plot==TRUE){
           Xv <- X[,d]
           plot(Xv, xlab=paste0("x", d[1]), ylab=paste0("x", d[2]), 
                xlim=c(0,1), ylim=c(0,1), main=paste0("md.loc at it=", t))
           points(Xorig[,d], pch=20, col="forestgreen")
           points(xold[d[1]], xold[d[2]], pch=20, col="red")
           points(xd[d[1]], xd[d[2]], pch=20, col="blue")
           rect(low[d[1]], low[d[2]], upp[d[1]], upp[d[2]], col=rgb(0,1,0,alpha=0.2))
           if(!is.null(Xorig)){
              legend("topright", c("xold", expression(x[d]), "Xorig"), 
                     xpd=TRUE, horiz=TRUE, inset=c(-0.015, -0.045), 
                     pch=20, col=c("red","blue", "forestgreen"), bty='n')
           }else{
              legend("topright", c("xold", expression(x[d])), 
                     xpd=TRUE, horiz=TRUE, inset=c(-0.015, -0.045), 
                     pch=20, col=c("red", "blue"), bty='n')
           }
        }
      }else{## after going through each run with ``md'' and with ``non-md'': start "jumping out of the well", i.e.,
            ## let the stuck location jump to a place close to the X run having the maximum ``md'' to other X runs and Xorig. 
            ## Such a way guarantees the start location no longer be stuck, and thus, there is no necessity to be iterative;
            ## the 3rd strategy
        
        ## justify ``xs''
        ## for each location in X, calculate its ``md'' to other X runs and Xorig
        me <- apply(DD2, 1, function(x) min(x[x>0]))
        ## justify the index with the maximum `md' 
        max.ind <- which(rowSums(DD2 < max(me))==1)
        inds <- ceiling(runif(1)*length(max.ind))
        max.inds <- max.ind[inds]
        xs <- X[max.inds,]
        
        ## justify ``xd''
        Dor <- order(DD2[max.inds,], decreasing=FALSE)
        xd <- ifelse(rep(Dor[2]<=n, p), X[Dor[2],], Xorig[Dor[2]-n,])
        
        ## define the search space
        sw.in <- sp(xs=xs, xd=xd, md=NULL)
        low <- sw.in$lo
        upp <- sw.in$up
        
        ## the new start location
        st.new <- st(xs=xs, xd=xd, ds=sw.in$ds)
        
        ## update X
        X[row.out,] <- st.new ## use the "new" location to replace the stuck location
        
        start <- matrix(X[row.out,], nrow=1)
        out <- optim(par=start, fn=maximin.obj, method="L-BFGS-B", 
                     lower=low, upper=upp, control=list(fnscale=-1), 
                     X=X, row.out=row.out, Xorig=Xorig)
        
        ## plot ``xs'' and ``xd''
        if(plot==TRUE){
           Xv <- X[,d]
           plot(Xv, xlab=paste0("x", d[1]), ylab=paste0("x", d[2]), 
                xlim=c(0,1), ylim=c(0,1), main=paste0("md.loc at it=", t))
           points(Xorig[,d], pch=20, col="forestgreen")
           points(xold[[1]], xold[d[2]], pch=20, col="cyan")
           points(xs[d[1]], xs[d[2]], pch=20, col="blue")
           points(st.new[d[1]], st.new[d[2]], pch=20, col="red")
           points(xd[d[1]], xd[d[2]], pch=20, col="orange")
           rect(low[d[1]], low[d[2]], upp[d[1]], upp[d[2]], col=rgb(0,1,0,alpha=0.2))
           arrows(xold[d[[1]]], xold[d[[2]]], st.new[d[[1]]], st.new[d[[2]]], lty=2, col="gray")
           if(!is.null(Xorig)){
              legend("topright", c("st.new", "Xorig"), xpd=TRUE, 
                     horiz=TRUE, inset=c(-0.015, -0.045), pch=20, col=c("red","forestgreen"), bty='n')
           }else{
              legend("topright", "st.new", xpd=TRUE, 
                     horiz=TRUE, inset=c(-0.015, -0.045), pch=20, col="red", bty='n')
           }
        }
      }
    }
    
    if(plot==TRUE){
       xnew <- as.numeric(out$par)
       Xv <- X[,d]
       plot(Xv, xlab=paste0("x", d[1]), ylab=paste0("x", d[2]), 
            xlim=c(0,1), ylim=c(0,1), main=paste0("md.loc at it=", t))
       points(Xorig[,d], pch=20, col="forestgreen")
       points(X[row.out,][d[1]], X[row.out,][d[2]], pch=20, col="red")
       points(xnew[d[1]], xnew[d[2]], pch=20, col="blue")
       rect(low[d[1]], low[d[2]], upp[d[1]], upp[d[2]], col=rgb(0,1,0,alpha=0.2))
       arrows(X[row.out,][d[1]], X[row.out,][d[2]], 
              xnew[d[1]], xnew[d[2]], length=0.1, lty=2, col="black") ## 'start' --> 'xnew'
       if(!is.null(Xorig)){
          legend("topright", c("xnew", "Xorig"), xpd=TRUE, horiz=TRUE, 
                 inset=c(-0.015, -0.045), pch=20, col=c("blue","forestgreen"), bty='n')
       }else{
          legend("topright", "xnew", xpd=TRUE, horiz=TRUE, 
                 inset=c(-0.015, -0.045), pch=20, col="blue", bty='n')
       }
    }
    
    X[row.out,] <- out$par
    D[row.out, -row.out] <- D[-row.out, row.out] <- as.numeric(distance(matrix(X[row.out,], nrow=1), 
                                                                        X[-row.out,]))    # update D
    md <- min(as.numeric(D[upper.tri(D)]))
    md.ind <- which(D==md, arr.ind=TRUE)[,1]
    
    ## distances to the existing design
    if(!is.null(Xorig)){
       D2[row.out,] <- as.numeric(distance(matrix(X[row.out,], nrow=1), Xorig)) ## update D2
       md2 <- min(D2)
       if(md2 < md){
          md <- md2
          md.ind <- which(D2==md2, arr.ind=TRUE)[,1]
       }else if(md2==md){
          md.ind <- unique(c(md.ind, which(D2==md2, arr.ind=TRUE)[,1]))
       }
    }
    
    if(plot==TRUE){
       mtext(paste0("md = ", format(round(as.numeric(md), 5), nsmall=5)))
    }
    
    # # the 4th strategy: not optimal!
    
    # ## When t >= 1.3 * n, if there is no improvement on ``md'' during the past
    # ## (0.3 * n) iterations, then the algorithm should be pushed to make progress ...
    # ## In addition, the push should not be performed during the last (0.2 * T - 1) iterations 
    # ## since such ``push'' only improves design criterion in the long term ...
    
    # if((t >= 1.3*n) && (t <= 0.8*T) && (cvs(mind[(t+1-0.3*n):(t+1)]) == TRUE)){
    #   
    #   if(is.null(Xorig)==TRUE)
    #      D2 <- NULL
    #
    #   DD2 <- cbind(D, D2)
    #   
    #   ## justify ``xs''
    #   me <- apply(DD2, 1, function(x) min(x[x>0]))
    #   max.ind <- which(rowSums(DD2 < max(me))==1)
    #   inds <- ceiling(runif(1)*length(max.ind))
    #   max.inds <- max.ind[inds]
    #   xs <- X[max.inds,]
    #   
    #   ## justify ``xd''
    #   Dor <- order(DD2[max.inds,], decreasing=FALSE)
    #   xd <- ifelse(rep(Dor[2] <= n, p), X[Dor[2],], Xorig[Dor[2]-n,])
    #   
    #   ## justify the search space
    #   sw.in <- sp(xs=xs, xd=xd, md=NULL)
    #   low <- sw.in$lo
    #   upp <- sw.in$up
    #   
    #   ## the new start location
    #   st.new <- st(xs=xs, xd=xd, ds=sw.in$ds)
    #   
    #   ## update X
    #   X[row.out,] <- st.new
    #   
    #   ## start optimization
    #   start <- matrix(X[row.out,], nrow=1)
    #   out <- optim(par=start, fn=maximin.obj, method="L-BFGS-B", 
    #                lower=low, upper=upp, control=list(fnscale=-1), 
    #                X=X, row.out=row.out, Xorig=Xorig)
    #   
    #   X[row.out,] <- out$par
    #   D[row.out, -row.out] <- D[-row.out, row.out] <- as.numeric(distance(matrix(X[row.out,], nrow=1), 
    #                                                                       X[-row.out,])) ## update D
    #   md <- min(as.numeric(D[upper.tri(D)]))
    #   md.ind <- which(D==md, arr.ind=TRUE)[,1]
    #   
    #   ## distances to the existing design
    #   if(!is.null(Xorig)){
    #      D2[row.out,] <- as.numeric(distance(matrix(X[row.out,], nrow=1), 
    #                                          Xorig)) ## update D2
    #      md2 <- min(D2)
    #      if(md2 < md){
    #         md <- md2
    #         md.ind <- which(D2==md2, arr.ind=TRUE)[,1]
    #      }else if(md2==md){
    #         md.ind <- unique(c(md.ind, which(D2==md2, arr.ind=TRUE)[,1]))
    #      }
    #   }
    # }
    
    # save important metrics, such as `md' and `X', at each iteration: not sure whether X should be saved with each iteration
    mind[t+1] <- md
    Xind[[t+1]] <- X
    
    if((verb == TRUE) && (t%%10 == 0)) 
        cat("t=", t, "/T=", T, " is done.\n", sep="")
  }
  
  # # At the last iteration, retrieve the maximum ``md'' so far and the corresponding `X'---?
  mind[t+1] <- max(mind) # not sure it is needed, 
  X <- Xind[[t+1]]
  
  ## combine `Xorig' and `X'
  XXorig <- rbind(Xorig, X)
  
  return(list(Xf=XXorig, mi=mind))
}