R/OLHD.C2007.R

In LHD: Latin Hypercube Designs (LHDs)

#' Orthogonal Latin Hypercube Design
#'
#' \code{OLHD.C2007} returns a \code{2^m+1} by \code{m+{m-1 \\choose 2}} orthogonal Latin hypercube design generated by the construction method of Cioppa and Lucas (2007)
#'
#' @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=m+{m-1 \\choose 2}}.
#'
#' @references Cioppa, T.M., and Lucas, T.W. (2007) Efficient nearly orthogonal and space-filling Latin hypercubes. \emph{Technometrics}, \strong{49}(1), 45-55.
#'
#' @examples
#' #create an orthogonal LHD with m=4. So n=2^m+1=17 and k=4+3=7
#' OLHD.C2007(m=4)
#'
#' #create an orthogonal LHD with m=5. So n=2^m+1=33 and k=5+6=11
#' OLHD.C2007(m=5)
#'
#' @export

OLHD.C2007=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=seq(1,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)) {
      for (j in (i+1):(m-1)) {

        M=cbind(M,AL[,,i]%*%AL[,,j]%*%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 1:(m-2)) {
      for (j in (i+1):(m-1)) {

        S=cbind(S,ak[,,i]*ak[,,j])

      }
    }
  }

  #construction of S ends


  #construction of T starts

  if(m==2){
    T0=matrix(0,nrow=q,ncol=2)

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

    #construction of T ends

    CP=matrix(0,ncol=2,nrow=1)   #center point
  }

  if(m>=3){
    T0=matrix(0,nrow=q,ncol=(m+choose((m-1),2)))

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

    #construction of T ends

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

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