R/OLHD.Y1998.R

In LHD: Latin Hypercube Designs (LHDs)

#' Orthogonal Latin Hypercube Design
#'
#' \code{OLHD.Y1998} returns a \code{2^m+1} by \code{2m-2} orthogonal Latin hypercube design generated by the construction method of Ye (1998)
#'
#' @param m A positive integer, and it must be greater than or equal to 2.
#'
#' @return If all inputs are logical, then the output will be an orthogonal LHD with the following run size: \code{n=2^m+1} and the following factor size: \code{k=2m-2}.
#'
#' @references Ye, K.Q. (1998) Orthogonal column Latin hypercubes and their application in computer experiments. \emph{Journal of the American Statistical Association}, \strong{93}(444), 1430-1439.
#'
#' @examples
#' #create an orthogonal LHD with m=3. So n=2^m+1=9 and k=2*m-2=4
#' OLHD.Y1998(m=3)
#'
#' #create an orthogonal LHD with m=4. So n=2^m+1=17 and k=2*m-2=6
#' OLHD.Y1998(m=4)
#'
#' @export

OLHD.Y1998=function(m){

  if(m<2){
    stop("m must be greater than or equal to 2")
  }

  #m >= 2

  q=2^(m-1)

  #construction of M starts

  e=sample(seq(1,q),q)        #first column
  e=matrix(e,ncol=1)

  I=matrix(c(1,0,0,1),ncol=2,nrow=2,byrow=T)
  R=matrix(c(0,1,1,0),ncol=2,nrow=2,byrow=T)

  AL=array(0,dim=c(q,q,m-1))  #there are m-1 of AL's

  if(m==2){

    AL[,,m-1]=R

    M=cbind(e,AL[,,m-1]%*%e)

  }

  if(m>=3){

    for (i in 1:(m-1-1)) {    #This is the loop for L. However, When L=m-1, m-1-L=0, and AL cannot be
      #calcultated within the loop. So L=1,...,m-2, and AL with L=m-1 will be calculated separately
      a=1

      for (j in 1:(m-1-i)) {  #This is the loop for m-1-L, the first half of AL
        a=kronecker(a,I)
      }

      b=1

      for (k in 1:i) {
        b=kronecker(b,R)
      }

      AL[,,i]=kronecker(a,b)

    }

    c=1

    for (l in 1:(m-1)) {
      c=kronecker(c,R)
    }

    AL[,,m-1]=c

    M=e

    for (i in 1:(m-1)) {
      M=cbind(M,AL[,,i]%*%e)
    }

    for (i in 1:(m-2)) {
      M=cbind(M,AL[,,i]%*%AL[,,m-1]%*%e)
    }
  }

  #construction of M ends

  #construction of S starts
  j=rep(1,q)
  j=matrix(j,ncol=1)

  ak=array(0,dim=c(q,1,m-1))  #there are m-1 of ak's

  B=array(1,dim=c(2,1,m-1))

  if(m==2){
    B[1,,m-1]=-1           #B_{m-k}
    ak[,,m-1]=B[,,1]
    S=cbind(j,ak[,,m-1])
  }

  if(m>=3){

    for (i in 1:(m-1)) {   #This is the loop for k.
      temp=B
      temp[1,,m-i]=-1      #B_{m-k}

      d=1

      for (k in 1:(m-1)) {
        d=kronecker(d,temp[,,k])
      }

      ak[,,i]=d
    }

    S=j

    for (i in 1:(m-1)) {
      S=cbind(S,ak[,,i])
    }

    for (i in 2:(m-1)) {
      S=cbind(S,ak[,,1]*ak[,,i])
    }


  }

  #construction of S ends


  #construction of T starts

  T0=matrix(0,nrow=q,ncol=(2*m-2))

  for (i in 1:q) {
    for (k in 1:(2*m-2)) {
      T0[i,k]=M[i,k]*S[i,k]
    }
  }

  #construction of T ends

  CP=matrix(0,ncol=(2*m-2),nrow=1)   #center point

  X=rbind(T0,CP,-T0)
  X
}