R/FastMmLHD.R

In LHD: Latin Hypercube Designs (LHDs)

#' Fast Maximin Distance LHD
#'
#' \code{FastMmLHD} returns a \code{n} by \code{k} maximin distance LHD matrix generated by the construction method of Wang, L., Xiao, Q., and Xu, H. (2018)
#'
#' @param n A positive integer, which stands for the number of rows (or run size).
#' @param k A positive integer, which stands for the number of columns (or factor size).
#' @param method A distance measure method. The default setting is "manhattan", and it could be one of the following: "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski". Any unambiguous substring can be given.
#' @param t1 A tunning parameter, which determines how many repeats will be implemented to search for the optimal design. The default is set to be 10.
#'
#' @return If all inputs are logical, then the output will be a \code{n} by \code{k} maximin distance LHD under under the maximin L_1 distance criterion..
#'
#' @references Wang, L., Xiao, Q., and Xu, H. (2018)  Optimal maximin $L_{1}$-distance Latin hypercube designs based on good lattice point designs. \emph{The Annals of Statistics}, \strong{46}(6B), 3741-3766.
#'
#' @examples
#' #n by n design when 2n+1 is prime
#' try=FastMmLHD(8,8)
#' try
#' phi_p(try)   #calculate the phi_p of "try".
#'
#' #n by n design when n+1 is prime
#' try2=FastMmLHD(12,12)
#' try2
#' phi_p(try2)   #calculate the phi_p of "try2".
#'
#' #n by n-1 design when n is prime
#' try3=FastMmLHD(7,6)
#' try3
#' phi_p(try3)   #calculate the phi_p of "try3".
#'
#' #General cases
#' try4=FastMmLHD(24,8)
#' try4
#' phi_p(try4)   #calculate the phi_p of "try4".
#'
#' @export

FastMmLHD=function(n,k,method="manhattan",t1=10){
  #n: number of runs
  #k: number of columns
  #method: the distance method, with default setting is "manhattan".
  #t1: tunning parameter

  if (!(class(n) == "numeric") || !(n%%1 == 0) || (n < 3))
  {
    stop("'n' must be a postive integer no less than 3")
  }
  if (!(class(k) == "numeric") || !(k%%1 == 0) || (k < 2))
  {
    stop("'k' must be a postive integer no less than 2")
  }
  if (n < k) {stop("n must be no less than k")}


  if (!is.na(pmatch(method, "euclidian")) || (method == "l2") || (method == "L2"))
  {method = "euclidean"}
  if (!is.na(pmatch(method, "manhattan"))|| (method == "l1") || (method == "L1"))
  {method = "manhattan"}

  METHODS = c("euclidean", "manhattan")
  count = pmatch(method, METHODS)

  if (is.na(count)==TRUE)
    stop("invalid distance method")
  if (count == -1)
    stop("ambiguous distance method")
  ###

  #the first special case of n \times n design when m = 2n+1 is prime. (2.3)
  if ((n == k) && numbers::isPrime(2*n + 1)==TRUE)
  {
    m = 2*n + 1

    GD=GLP(n=m,k=m-1,h=1:(m-1))  #GLP design

    RD=GD   #result design

    for (i in 1:m){
      for (j in 1:(m-1)){

        if(GD[i,j]<(m/2)){RD[i,j]=2*GD[i,j]}
        if(GD[i,j]>=(m/2)){RD[i,j]=2*(m-GD[i,j])}

      }
    }

    resultdesign = RD/2
    return(resultdesign[1:n, 1:n]-1)
  }

  #the second special case of n \times (k<=n) design when n+1 is prime. (2.2)
  else if ((n >= k) && numbers::isPrime(n + 1)==TRUE)
  {
    finaldesign = NULL
    finalvalue = -Inf

    for(i in 1:t1)
    {

      GD=GLP(n=n+1,k=k,h=sample(seq(from=1,to=n),k))    #starting design: D_b with b=0. GLP design

      TD=WT(GD,baseline=0)[-(n+1),]   #E_b with b=0. Transformed design

      minvalue=min(stats::dist(TD, method = method))  #current minimum distance

      for (c in 1:n)      #for b=1, ... , (n+1)-1
      {
        newdesign = (GD + c) %% (n+1)    #new D_b

        newresult = WT(newdesign,baseline=0)     #new E_b

        cutv = newresult[(n+1),1]   #cut value: this is value of the constant row

        newresult=newresult[-(n+1),] #get rid of the last row since it is constant

        for(i in 1:n){    #bring one number down for remaining elements > cutv
          for(j in 1:k){

            if(newresult[i,j] > cutv){
              newresult[i,j] = newresult[i,j] - 1
            }

          }
        }

        newminvalue = min(stats::dist(newresult, method = method))

        if (newminvalue > minvalue)      #updating minimum distance
        {
          resultdesign = newresult
          minvalue = newminvalue
        }

      }

      if (minvalue > finalvalue)
      {
        finaldesign = resultdesign
        finalvalue = minvalue
      }
    }

    return(finaldesign)
  }

  #the third specical case of n \times (k<n) design when n is prime. (2.1)
  else if ((n > k) && numbers::isPrime(n)==TRUE)
  {
    finaldesign = NULL
    finalvalue = -Inf

    for(i in 1:t1)
    {

      GD=GLP(n=n,k=k,h=sample(seq(from=1,to=(n-1)),k))    #starting design: D_b with b=0. GLP design

      TD=WT(GD,baseline=0)   #E_b with b=0. Transformed design

      minvalue=min(stats::dist(TD, method = method))  #current minimum distance

      for(c in 1:(n-1))   #for b=1, ... , n-1
      {
        newdesign = (GD + c) %% n      #new D_b

        newresult = WT(newdesign,baseline=0)   #new E_b

        newminvalue = min(stats::dist(newresult, method = method))

        if (newminvalue > minvalue)          #updating minimum distance
        {
          resultdesign = newresult
          minvalue = newminvalue
        }
      }

      if (minvalue > finalvalue)
      {
        finaldesign = resultdesign
        finalvalue = minvalue
      }
    }

    return(finaldesign)
  }

  else {
    #General case, when above special cases are not true.
    #coprime vector to n
    copn=NULL

    for(i in 1:(n-1)){

      if (numbers::coprime(i,n)==TRUE)
      {
        copn=c(copn,i)
      }
    }

    #Euler function value of n
    eul=length(copn)

    if (k>eul){
      stop(paste0("the number of column for this run size must be no greater than ", eul))
    }

    else {

      GD = GLP(n=n,k=eul,h=sample(copn,eul))

      TD = WT(GD,baseline=0)

      minvalue = min(stats::dist(TD, method = method))

      for (c in 1:(n-1))
      {
        newdesign = (GD + c) %% n

        newresult = WT(newdesign,baseline=0)

        newminvalue = min(stats::dist(newresult, method = method))

        if (newminvalue > minvalue)
        {
          resultdesign = newresult
          minvalue = newminvalue
        }
      }

      #This is the searching part: search for the best subset if k!=eul.
      if(k==eul){

        return(resultdesign)

      }

      else {
        Temp=NULL
        Temp.value=Inf

        for (i in 1:t1) {

          rcol=sample(1:eul,k)   #random columns

          X=resultdesign[,rcol]

          X.value=phi_p(X)

          if(X.value<=Temp.value){
            Temp=X
            Temp.value=X.value
          }
        }

        print("Design has been generated, but it may not be the best one for this size.")
        return(Temp)
      }

    }

  }

}