R/OLHD.S2010.R

In LHD: Latin Hypercube Designs (LHDs)

#' Orthogonal Latin Hypercube Design
#'
#' \code{OLHD.S2010} returns a \code{r2^{C+1}+1} or \code{r2^{C+1}} by \code{2^C} orthogonal Latin hypercube design generated by the construction method of Sun et al. (2010)
#'
#' @param C A positive integer.
#' @param r A positive integer.
#' @param type Design run size type, and it could be either odd or even. If \code{type} is odd (the default setting), \code{OLHD.S2010} returns an OLHD with run size \code{r2^{C+1}+1}. If \code{type} is even, \code{OLHD.S2010} returns an OLHD with run size \code{r2^{C+1}}.
#'
#' @return If all inputs are logical, then the output will be an orthogonal LHD with the following run size: \code{n=r2^{C+1}+1} or \code{n=r2^{C+1}} and the following factor size: \code{k=2^C}.
#'
#' @references Sun, F., Liu, M.Q., and Lin, D.K. (2010) Construction of orthogonal Latin hypercube designs with flexible run sizes. \emph{Journal of Statistical Planning and Inference}, \strong{140}(11), 3236-3242.
#'
#' @examples
#' #create an orthogonal LHD with C=3, r=3, type="odd".
#' #So n=3*2^4+1=49 and k=2^3=8
#' OLHD.S2010(C=3,r=3,type="odd")
#'
#' #create an orthogonal LHD with C=3, r=3, type="even".
#' #So n=3*2^4=48 and k=2^3=8
#' OLHD.S2010(C=3,r=3,type="even")
#'
#' @export

OLHD.S2010=function(C,r,type="odd"){

  Sc=matrix(c(1,1,1,-1),ncol=2,nrow=2,byrow=T) #S1
  Tc=matrix(c(1,2,2,-1),ncol=2,nrow=2,byrow=T) #T1

  if(C>=2){
    counter=2

    while(counter<=C) {

      Sc.star=Sc
      Tc.star=Tc

      index=nrow(Sc.star)/2

      for (i in 1:index) {

        Sc.star[i,]=Sc.star[i,]*(-1)
        Tc.star[i,]=Tc.star[i,]*(-1)

      }

      a=rbind(Sc,Sc)   #first half of Sc

      b=rbind(-Sc.star,Sc.star)   #second half of Sc

      c=rbind(Tc,Tc+2^(counter-1)*Sc)   #first half of Tc

      d=rbind(-1*(Tc.star+2^(counter-1)*Sc.star),Tc.star)   #second half of Tc

      Sc=cbind(a,b)
      Tc=cbind(c,d)


      counter=counter+1
    }
  }
  #until here, Sc and Tc are fully constructed

  if(type=="odd"){

    A=NULL

    for (i in 1:r) {
      Tc.i=Tc+(i-1)*2^C*Sc

      A=cbind(A,t(Tc.i))
    }

    A=t(A)

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

    X=rbind(A,CP,-A)

  }

  if(type=="even"){

    B=NULL

    Hc=Tc-Sc/2

    for (i in 1:r) {
      Hc.i=Hc+(i-1)*2^C*Sc

      B=cbind(B,t(Hc.i))
    }

    B=t(B)

    X=rbind(B,-B)

  }

  X
}