R/SLHD.R

In LHD: Latin Hypercube Designs (LHDs)

#' Sliced Latin Hypercube Design (SLHD)
#'
#' \code{SLHD} returns a \code{n} by \code{k} LHD matrix generated by improved two-stage algorithm
#'
#' @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 t A positive integer, which stands for the number of slices. \code{n}/\code{t} must be a positive integer, that is, n is divisible by t. \code{t} must be smaller or equal to \code{k} when \code{n} is 9 or larger. \code{t} must be smaller than \code{k} when \code{n} is smaller than 9. Otherwise, the funtion will never stop. The default is set to be 1.
#' @param N A positive integer, which stands for the number of iterations. The default is set to be 10. A large value of \code{N} will result a high CPU time, and it is recommended to be no greater than 500.
#' @param T0 A positive number, which stands for the user-defined initial temperature. The default is set to be 10.
#' @param rate A positive percentage, which stands for temperature decrease rate, and it should be in (0,1). For example, rate=0.25 means the temperature decreases by 25\% each time. The default is set to be 10\%.
#' @param Tmin A positive number, which stands for the minimium temperature allowed. When current temperature becomes smaller or equal to \code{Tmin}, the stopping criterion for current loop is met. The default is set to be 1.
#' @param Imax A positive integer, which stands for the maximum perturbations the algorithm will try without improvements before temperature is reduced. The default is set to be 5. For the computation complexity consideration, \code{Imax} is recommended to be smaller or equal to 5.
#' @param OC An optimality criterion. The default setting is "phi_p", and it could be one of the following: "phi_p", "AvgAbsCor", "MaxAbsCor", "MaxProCriterion".
#' @param p A positive integer, which is the parameter in the phi_p formula, and \code{p} is prefered to be large. The default is set to be 15.
#' @param q The default is set to be 1, and it could be either 1 or 2. If \code{q} is 1, \code{dij} is the Manhattan (rectangular) distance. If \code{q} is 2, \code{dij} is the Euclidean distance.
#' @param stage2 A logic input argument, and it could be either FALSE or TRUE. If \code{stage2} is FALSE (the default setting), \code{SLHD} will only implement the first stage of the algorithm. If \code{stage2} is TRUE, \code{SLHD} will implement the whole algorithm.
#' @param maxtime A positive number, which indicates the expected maximum CPU time given by user, and it is measured by minutes. For example, maxtime=3.5 indicates the CPU time will be no greater than three and half minutes. The default is set to be 5.
#'
#' @return If all inputs are logical, then the output will be a \code{n} by \code{k} LHD. As mentioned from the original paper, the first stage plays a much more important role since it optimizes the slice level. More resources should be given to the first stage if computational budgets are limited. Let m=n/t, where m is the number of rows for each slice, if (m)^k >> n, the second stage becomes optional. That is the reason why we add a \code{stage2} parameter to let users decide if the second stage is needed.
#'
#' @references Ba, S., Myers, W.R., and Brenneman, W.A. (2015) Optimal Sliced Latin Hypercube Designs. \emph{Technometrics}, \strong{57}, 479-487.
#'
#' @examples
#' #generate a 5 by 3 maximin distance LHD with the default setting
#' try=SLHD(n=5,k=3)
#' try
#' phi_p(try)   #calculate the phi_p of "try".
#'
#' #generate a 5 by 3 maximin distance LHD with stage II
#' #let stage2=TRUE and other input are the same as above
#' try2=SLHD(n=5,k=3,stage2=TRUE)
#' try2
#' phi_p(try2)   #calculate the phi_p of "try2".
#'
#' #Another example
#' #generate a 8 by 4 nearly orthogonal LHD
#' try3=SLHD(n=8,k=4,OC="AvgAbsCor",stage2=TRUE)
#' try3
#' AvgAbsCor(try3)  #calculate the average absolute correlation.
#' @export

SLHD=function(n,k,t=1,N=10,T0=10,rate=0.1,Tmin=1,Imax=3,OC="phi_p",p=15,q=1,stage2=FALSE,maxtime=5){
  #n and k are the rs and fa.
  #t: number of slices. n/t must be an integer, that is, n is divisible by t.
  #N: maximum number of iterations.
  #T0: initial temperature
  #rate: temperature decrease rate. 0<rate<1
  #Tmin: minumum temperature for each itertaion,TPmin > 0
  #Imax:# of perturbations the algorithm will try without improvements before Temperature is reduced
  #OC: optimality criterion, the default is "phi_p", along with default p and q
  #stage2: if stage II of the algorithm will be performed. The default is FALSE (not-run).
  #Note that stage I plays a much more important role in the whole algorithm. More resources should
  #be allocated for stage I if computational budgets are tight.

  maxtime=maxtime*60  #convert minutes to seconds
  timeALL1=NULL        #record all cpu time
  timeALL2=NULL        #record all cpu time

  C=1  #Initialize counter index

  m=n/t     #the rs for each slice.

  #step 1 on page 481 starts
  Y=rep(0,n*k)

  dim(Y)=c(m,k,t)

  #Independtly generate t small LHD for t slices. Each slice is an m by k LHD.
  for (i in 1:t) {
    Y[,,i]=rLHD(m,k)
  }

  #stack the t slices to form an n by k LHD matrix.
  X=NULL

  for (i in 1:t) {
    X=rbind(X,Y[,,i])
  }

  #step 1 on page 481 ends


  #The improved two-stage algorithms starts

  #Stage I starts

  SX=0      #the S(X) in stage I
  DR=NULL   #this is the duplicated rows index, which is used to record row numbers.

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

    for (j in (i+1):n) {
      I=dij(X=X,i=i,j=j,q=q)
      if (I==0){SX=SX+1;DR=c(DR,i,j)}
    }

  }

  DR=unique(DR)

  #step ii starts: if S(X)>0
  while (SX>0){
    rrow1=sample(DR,1)
    slice=ceiling(rrow1/m)     #locate the rrow1's slice

    #select another row within the same slice
    rrow2=sample(seq(from=slice*m,by=-1,length.out=m)[seq(from=slice*m,by=-1,length.out=m)!=rrow1],1)

    rcol=sample(1:k,1)

    Xnew=X

    e1=Xnew[rrow1,rcol]            #exchange 2 elements to form Xnew
    e2=Xnew[rrow2,rcol]

    Xnew[rrow1,rcol]=e2
    Xnew[rrow2,rcol]=e1

    SXnew=0      #S(Xnew)
    DRnew=NULL   #DR for Xnew

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

      for (j in (i+1):n) {
        I=dij(X=Xnew,i=i,j=j,q=q)
        if (I==0){SXnew=SXnew+1;DRnew=c(DRnew,i,j)}
      }

    }

    DRnew=unique(DRnew)

    if (SXnew<SX){X=Xnew;SX=SXnew;DR=DRnew}    #this is step iii

  }
  #step ii and iii ends

  Xbest=X;TP=T0;Flag=1

  if(OC=="phi_p"){

    progressbar = utils::txtProgressBar(min = 0, max = N, style = 3)

    while (C<=N) {

      time0=Sys.time()

      while(Flag==1 & TP>Tmin){
        Flag=0;I=1

        while (I<=Imax) {

          #at this point, S(X)==0 already.

          #step iv starts

          rt=sample(1:t,1)        #randomly select a slice from the current X

          rcol=sample(1:k,1)      #randomly choose a column

          rrow=sample(seq(from=rt*m,by=-1,length.out=m),2)    #randomly choose 2 rows from slice rt

          Xnew=X

          e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
          e2=Xnew[rrow[2],rcol]

          Xnew[rrow[1],rcol]=e2
          Xnew[rrow[2],rcol]=e1
          #step iv ends

          SXnew=0         #S(Xnew)

          #step v starts
          for (i in 1:(n-1)) {

            for (j in (i+1):n) {
              I=dij(X=Xnew,i=i,j=j,q=q)
              if (I==0){SXnew=SXnew+1}
            }

          }

          while (SXnew>0) {

            rt=sample(1:t,1)        #randomly select a slice from the current X

            rcol=sample(1:k,1)      #randomly choose a column

            rrow=sample(seq(from=rt*m,by=-1,length.out=m),2)    #randomly choose 2 rows from slice rt

            Xnew=X

            e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
            e2=Xnew[rrow[2],rcol]

            Xnew[rrow[1],rcol]=e2
            Xnew[rrow[2],rcol]=e1

            SXnew=0

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

              for (j in (i+1):n) {
                I=dij(X=Xnew,i=i,j=j,q=q)
                if (I==0){SXnew=SXnew+1}
              }
            }

          }
          #step v ends


          #step vi starts

          a=phi_p(X=Xnew,p=p,q=q)
          b=phi_p(X=X,p=p,q=q)
          if (a<b){X=Xnew;Flag=1}
          if (a>=b){
            prob=exp((b-a)/TP)
            draw=sample(c(0,1),1,prob=c(1-prob,prob))    #draw==1 means replace
            if(draw==1){X=Xnew;Flag=1}
          }                         #step 5 ends here

          c=phi_p(X=Xbest,p=p,q=q)
          if (a<c){Xbest=Xnew;I=1}
          if (a>=c){I=I+1}

          #step vi ends

        }

        TP=TP*(1-rate)
      }

      time1=Sys.time()
      timediff=time1-time0
      timeALL1=c(timeALL1,timediff)

      ##########progress bar codes
      utils::setTxtProgressBar(progressbar, C)
      ##########

      timeALL=sum(timeALL1)+sum(timeALL2)

      if(as.numeric(timeALL+timediff)<=maxtime){C=C+1}
      if(as.numeric(timeALL+timediff)>maxtime){C=N+1}

      TP=T0;Flag=1
    }

    #Stage I ends

    #step 2 on page 481 starts
    if(t>1){

      pi_l=rep(0,n)   #the \Pi_l for l=1, ... , m

      dim(pi_l)=c(t,1,m)

      for (i in 1:m) {
        pi_l[,,i]=seq(from=(i-1)*t+1,to=i*t,1)
      }

      for (j in 1:k) {
        for (i in 1:m) {
          Xbest[,j][Xbest[,j]==i]=sample(pi_l[,,i])*100
        }
      }
      Xbest=Xbest/100
    }

    #step 2 on page 481 ends


    #Stage II starts

    if (stage2==TRUE){
      C=1

      X=Xbest

      TP=T0;Flag=1

      while (C<=N) {

        time0=Sys.time()

        while(Flag==1 & TP>Tmin){
          Flag=0;I=1

          while (I<=Imax) {
            z=stats::runif(1,0,1)     #step i

            #step ii
            if (z<=0){

              rt=sample(1:t,1)        #randomly select a slice from X

              rcol=sample(1:k,1)      #randomly choose a column

              rrow=sample(seq(from=rt*m,by=-1,length.out=m),2)    #randomly choose 2 rows from slice rt

              Xnew=X

              e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
              e2=Xnew[rrow[2],rcol]

              Xnew[rrow[1],rcol]=e2
              Xnew[rrow[2],rcol]=e1
            }

            #step iii
            if (z>0){

              rcol=sample(1:k,1)      #randomly choose a column

              if(t>1){

                rl=sample(1:m,1)        #randomly choose a l, where l=1, ..., m

                re=sample(pi_l[,,rl],2)    #randomly choose 2 elements from pi_rl that will be exchanged

                Xnew=X

                rrow=c(which(Xnew[,rcol]==re[1]),which(Xnew[,rcol]==re[2]))

                e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
                e2=Xnew[rrow[2],rcol]

                Xnew[rrow[1],rcol]=e2
                Xnew[rrow[2],rcol]=e1

              }

              if(t==1){
                rrow=sample(1:n,2)       #if there is only one slice

                Xnew=X

                e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
                e2=Xnew[rrow[2],rcol]

                Xnew[rrow[1],rcol]=e2
                Xnew[rrow[2],rcol]=e1
              }


            }

            #step iv

            a=phi_p(X=Xnew,p=p,q=q)
            b=phi_p(X=X,p=p,q=q)
            if (a<b){X=Xnew;Flag=1}
            if (a>=b){
              prob=exp((b-a)/TP)
              draw=sample(c(0,1),1,prob=c(1-prob,prob))    #draw==1 means replace
              if(draw==1){X=Xnew;Flag=1}
            }                         #step 5 ends here

            c=phi_p(X=Xbest,p=p,q=q)
            if (a<c){Xbest=Xnew;I=1}
            if (a>=c){I=I+1}

          }

          TP=TP*(1-rate)
        }

        time1=Sys.time()
        timediff=time1-time0
        timeALL2=c(timeALL2,timediff)

        ##########progress bar codes
        utils::setTxtProgressBar(progressbar, C)
        ##########


        timeALL=sum(timeALL1)+sum(timeALL2)

        if(as.numeric(timeALL+timediff)<=maxtime){C=C+1}
        if(as.numeric(timeALL+timediff)>maxtime){C=N+1}

        TP=T0;Flag=1

      }

    }

    #Stage II ends
  }

  if(OC=="AvgAbsCor"){

    progressbar = utils::txtProgressBar(min = 0, max = N, style = 3)

    while (C<=N) {

      time0=Sys.time()

      while(Flag==1 & TP>Tmin){
        Flag=0;I=1

        while (I<=Imax) {

          #at this point, S(X)==0 already.

          #step iv starts

          rt=sample(1:t,1)        #randomly select a slice from the current X

          rcol=sample(1:k,1)      #randomly choose a column

          rrow=sample(seq(from=rt*m,by=-1,length.out=m),2)    #randomly choose 2 rows from slice rt

          Xnew=X

          e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
          e2=Xnew[rrow[2],rcol]

          Xnew[rrow[1],rcol]=e2
          Xnew[rrow[2],rcol]=e1
          #step iv ends

          SXnew=0         #S(Xnew)

          #step v starts
          for (i in 1:(n-1)) {

            for (j in (i+1):n) {
              I=dij(X=Xnew,i=i,j=j,q=q)
              if (I==0){SXnew=SXnew+1}
            }

          }

          while (SXnew>0) {

            rt=sample(1:t,1)        #randomly select a slice from the current X

            rcol=sample(1:k,1)      #randomly choose a column

            rrow=sample(seq(from=rt*m,by=-1,length.out=m),2)    #randomly choose 2 rows from slice rt

            Xnew=X

            e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
            e2=Xnew[rrow[2],rcol]

            Xnew[rrow[1],rcol]=e2
            Xnew[rrow[2],rcol]=e1

            SXnew=0

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

              for (j in (i+1):n) {
                I=dij(X=Xnew,i=i,j=j,q=q)
                if (I==0){SXnew=SXnew+1}
              }
            }

          }
          #step v ends


          #step vi starts

          a=AvgAbsCor(X=Xnew)
          b=AvgAbsCor(X=X)
          if (a<b){X=Xnew;Flag=1}
          if (a>=b){
            prob=exp((b-a)/TP)
            draw=sample(c(0,1),1,prob=c(1-prob,prob))    #draw==1 means replace
            if(draw==1){X=Xnew;Flag=1}
          }                         #step 5 ends here

          c=AvgAbsCor(X=Xbest)
          if (a<c){Xbest=Xnew;I=1}
          if (a>=c){I=I+1}

          #step vi ends

        }

        TP=TP*(1-rate)
      }

      time1=Sys.time()
      timediff=time1-time0
      timeALL1=c(timeALL1,timediff)

      ##########progress bar codes
      utils::setTxtProgressBar(progressbar, C)
      ##########

      timeALL=sum(timeALL1)+sum(timeALL2)

      if(as.numeric(timeALL+timediff)<=maxtime){C=C+1}
      if(as.numeric(timeALL+timediff)>maxtime){C=N+1}

      TP=T0;Flag=1
    }

    #Stage I ends

    #step 2 on page 481 starts
    if(t>1){

      pi_l=rep(0,n)   #the \Pi_l for l=1, ... , m

      dim(pi_l)=c(t,1,m)

      for (i in 1:m) {
        pi_l[,,i]=seq(from=(i-1)*t+1,to=i*t,1)
      }

      for (j in 1:k) {
        for (i in 1:m) {
          Xbest[,j][Xbest[,j]==i]=sample(pi_l[,,i])*100
        }
      }
      Xbest=Xbest/100
    }

    #step 2 on page 481 ends


    #Stage II starts

    if (stage2==TRUE){
      C=1

      X=Xbest

      TP=T0;Flag=1

      while (C<=N) {

        time0=Sys.time()

        while(Flag==1 & TP>Tmin){
          Flag=0;I=1

          while (I<=Imax) {
            z=stats::runif(1,0,1)     #step i

            #step ii
            if (z<=0){

              rt=sample(1:t,1)        #randomly select a slice from X

              rcol=sample(1:k,1)      #randomly choose a column

              rrow=sample(seq(from=rt*m,by=-1,length.out=m),2)    #randomly choose 2 rows from slice rt

              Xnew=X

              e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
              e2=Xnew[rrow[2],rcol]

              Xnew[rrow[1],rcol]=e2
              Xnew[rrow[2],rcol]=e1
            }

            #step iii
            if (z>0){

              rcol=sample(1:k,1)      #randomly choose a column

              if(t>1){

                rl=sample(1:m,1)        #randomly choose a l, where l=1, ..., m

                re=sample(pi_l[,,rl],2)    #randomly choose 2 elements from pi_rl that will be exchanged

                Xnew=X

                rrow=c(which(Xnew[,rcol]==re[1]),which(Xnew[,rcol]==re[2]))

                e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
                e2=Xnew[rrow[2],rcol]

                Xnew[rrow[1],rcol]=e2
                Xnew[rrow[2],rcol]=e1

              }

              if(t==1){
                rrow=sample(1:n,2)       #if there is only one slice

                Xnew=X

                e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
                e2=Xnew[rrow[2],rcol]

                Xnew[rrow[1],rcol]=e2
                Xnew[rrow[2],rcol]=e1
              }


            }

            #step iv

            a=AvgAbsCor(X=Xnew)
            b=AvgAbsCor(X=X)
            if (a<b){X=Xnew;Flag=1}
            if (a>=b){
              prob=exp((b-a)/TP)
              draw=sample(c(0,1),1,prob=c(1-prob,prob))    #draw==1 means replace
              if(draw==1){X=Xnew;Flag=1}
            }                         #step 5 ends here

            c=AvgAbsCor(X=Xbest)
            if (a<c){Xbest=Xnew;I=1}
            if (a>=c){I=I+1}

          }

          TP=TP*(1-rate)
        }

        time1=Sys.time()
        timediff=time1-time0
        timeALL2=c(timeALL2,timediff)

        ##########progress bar codes
        utils::setTxtProgressBar(progressbar, C)
        ##########


        timeALL=sum(timeALL1)+sum(timeALL2)

        if(as.numeric(timeALL+timediff)<=maxtime){C=C+1}
        if(as.numeric(timeALL+timediff)>maxtime){C=N+1}

        TP=T0;Flag=1

      }

    }

    #Stage II ends
  }

  if(OC=="MaxAbsCor"){

    progressbar = utils::txtProgressBar(min = 0, max = N, style = 3)

    while (C<=N) {

      time0=Sys.time()

      while(Flag==1 & TP>Tmin){
        Flag=0;I=1

        while (I<=Imax) {

          #at this point, S(X)==0 already.

          #step iv starts

          rt=sample(1:t,1)        #randomly select a slice from the current X

          rcol=sample(1:k,1)      #randomly choose a column

          rrow=sample(seq(from=rt*m,by=-1,length.out=m),2)    #randomly choose 2 rows from slice rt

          Xnew=X

          e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
          e2=Xnew[rrow[2],rcol]

          Xnew[rrow[1],rcol]=e2
          Xnew[rrow[2],rcol]=e1
          #step iv ends

          SXnew=0         #S(Xnew)

          #step v starts
          for (i in 1:(n-1)) {

            for (j in (i+1):n) {
              I=dij(X=Xnew,i=i,j=j,q=q)
              if (I==0){SXnew=SXnew+1}
            }

          }

          while (SXnew>0) {

            rt=sample(1:t,1)        #randomly select a slice from the current X

            rcol=sample(1:k,1)      #randomly choose a column

            rrow=sample(seq(from=rt*m,by=-1,length.out=m),2)    #randomly choose 2 rows from slice rt

            Xnew=X

            e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
            e2=Xnew[rrow[2],rcol]

            Xnew[rrow[1],rcol]=e2
            Xnew[rrow[2],rcol]=e1

            SXnew=0

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

              for (j in (i+1):n) {
                I=dij(X=Xnew,i=i,j=j,q=q)
                if (I==0){SXnew=SXnew+1}
              }
            }

          }
          #step v ends


          #step vi starts

          a=MaxAbsCor(X=Xnew)
          b=MaxAbsCor(X=X)
          if (a<b){X=Xnew;Flag=1}
          if (a>=b){
            prob=exp((b-a)/TP)
            draw=sample(c(0,1),1,prob=c(1-prob,prob))    #draw==1 means replace
            if(draw==1){X=Xnew;Flag=1}
          }                         #step 5 ends here

          c=MaxAbsCor(X=Xbest)
          if (a<c){Xbest=Xnew;I=1}
          if (a>=c){I=I+1}

          #step vi ends

        }

        TP=TP*(1-rate)
      }

      time1=Sys.time()
      timediff=time1-time0
      timeALL1=c(timeALL1,timediff)

      ##########progress bar codes
      utils::setTxtProgressBar(progressbar, C)
      ##########

      timeALL=sum(timeALL1)+sum(timeALL2)

      if(as.numeric(timeALL+timediff)<=maxtime){C=C+1}
      if(as.numeric(timeALL+timediff)>maxtime){C=N+1}

      TP=T0;Flag=1
    }

    #Stage I ends

    #step 2 on page 481 starts
    if(t>1){

      pi_l=rep(0,n)   #the \Pi_l for l=1, ... , m

      dim(pi_l)=c(t,1,m)

      for (i in 1:m) {
        pi_l[,,i]=seq(from=(i-1)*t+1,to=i*t,1)
      }

      for (j in 1:k) {
        for (i in 1:m) {
          Xbest[,j][Xbest[,j]==i]=sample(pi_l[,,i])*100
        }
      }
      Xbest=Xbest/100
    }

    #step 2 on page 481 ends


    #Stage II starts

    if (stage2==TRUE){
      C=1

      X=Xbest

      TP=T0;Flag=1

      while (C<=N) {

        time0=Sys.time()

        while(Flag==1 & TP>Tmin){
          Flag=0;I=1

          while (I<=Imax) {
            z=stats::runif(1,0,1)     #step i

            #step ii
            if (z<=0){

              rt=sample(1:t,1)        #randomly select a slice from X

              rcol=sample(1:k,1)      #randomly choose a column

              rrow=sample(seq(from=rt*m,by=-1,length.out=m),2)    #randomly choose 2 rows from slice rt

              Xnew=X

              e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
              e2=Xnew[rrow[2],rcol]

              Xnew[rrow[1],rcol]=e2
              Xnew[rrow[2],rcol]=e1
            }

            #step iii
            if (z>0){

              rcol=sample(1:k,1)      #randomly choose a column

              if(t>1){

                rl=sample(1:m,1)        #randomly choose a l, where l=1, ..., m

                re=sample(pi_l[,,rl],2)    #randomly choose 2 elements from pi_rl that will be exchanged

                Xnew=X

                rrow=c(which(Xnew[,rcol]==re[1]),which(Xnew[,rcol]==re[2]))

                e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
                e2=Xnew[rrow[2],rcol]

                Xnew[rrow[1],rcol]=e2
                Xnew[rrow[2],rcol]=e1

              }

              if(t==1){
                rrow=sample(1:n,2)       #if there is only one slice

                Xnew=X

                e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
                e2=Xnew[rrow[2],rcol]

                Xnew[rrow[1],rcol]=e2
                Xnew[rrow[2],rcol]=e1
              }


            }

            #step iv

            a=MaxAbsCor(X=Xnew)
            b=MaxAbsCor(X=X)
            if (a<b){X=Xnew;Flag=1}
            if (a>=b){
              prob=exp((b-a)/TP)
              draw=sample(c(0,1),1,prob=c(1-prob,prob))    #draw==1 means replace
              if(draw==1){X=Xnew;Flag=1}
            }                         #step 5 ends here

            c=MaxAbsCor(X=Xbest)
            if (a<c){Xbest=Xnew;I=1}
            if (a>=c){I=I+1}

          }

          TP=TP*(1-rate)
        }

        time1=Sys.time()
        timediff=time1-time0
        timeALL2=c(timeALL2,timediff)

        ##########progress bar codes
        utils::setTxtProgressBar(progressbar, C)
        ##########


        timeALL=sum(timeALL1)+sum(timeALL2)

        if(as.numeric(timeALL+timediff)<=maxtime){C=C+1}
        if(as.numeric(timeALL+timediff)>maxtime){C=N+1}

        TP=T0;Flag=1

      }

    }

    #Stage II ends
  }

  if(OC=="MaxProCriterion"){

    progressbar = utils::txtProgressBar(min = 0, max = N, style = 3)

    while (C<=N) {

      time0=Sys.time()

      while(Flag==1 & TP>Tmin){
        Flag=0;I=1

        while (I<=Imax) {

          #at this point, S(X)==0 already.

          #step iv starts

          rt=sample(1:t,1)        #randomly select a slice from the current X

          rcol=sample(1:k,1)      #randomly choose a column

          rrow=sample(seq(from=rt*m,by=-1,length.out=m),2)    #randomly choose 2 rows from slice rt

          Xnew=X

          e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
          e2=Xnew[rrow[2],rcol]

          Xnew[rrow[1],rcol]=e2
          Xnew[rrow[2],rcol]=e1
          #step iv ends

          SXnew=0         #S(Xnew)

          #step v starts
          for (i in 1:(n-1)) {

            for (j in (i+1):n) {
              I=dij(X=Xnew,i=i,j=j,q=q)
              if (I==0){SXnew=SXnew+1}
            }

          }

          while (SXnew>0) {

            rt=sample(1:t,1)        #randomly select a slice from the current X

            rcol=sample(1:k,1)      #randomly choose a column

            rrow=sample(seq(from=rt*m,by=-1,length.out=m),2)    #randomly choose 2 rows from slice rt

            Xnew=X

            e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
            e2=Xnew[rrow[2],rcol]

            Xnew[rrow[1],rcol]=e2
            Xnew[rrow[2],rcol]=e1

            SXnew=0

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

              for (j in (i+1):n) {
                I=dij(X=Xnew,i=i,j=j,q=q)
                if (I==0){SXnew=SXnew+1}
              }
            }

          }
          #step v ends


          #step vi starts

          a=MaxProCriterion(X=Xnew)
          b=MaxProCriterion(X=X)
          if (a<b){X=Xnew;Flag=1}
          if (a>=b){
            prob=exp((b-a)/TP)
            draw=sample(c(0,1),1,prob=c(1-prob,prob))    #draw==1 means replace
            if(draw==1){X=Xnew;Flag=1}
          }                         #step 5 ends here

          c=MaxProCriterion(X=Xbest)
          if (a<c){Xbest=Xnew;I=1}
          if (a>=c){I=I+1}

          #step vi ends

        }

        TP=TP*(1-rate)
      }

      time1=Sys.time()
      timediff=time1-time0
      timeALL1=c(timeALL1,timediff)

      ##########progress bar codes
      utils::setTxtProgressBar(progressbar, C)
      ##########

      timeALL=sum(timeALL1)+sum(timeALL2)

      if(as.numeric(timeALL+timediff)<=maxtime){C=C+1}
      if(as.numeric(timeALL+timediff)>maxtime){C=N+1}

      TP=T0;Flag=1
    }

    #Stage I ends

    #step 2 on page 481 starts
    if(t>1){

      pi_l=rep(0,n)   #the \Pi_l for l=1, ... , m

      dim(pi_l)=c(t,1,m)

      for (i in 1:m) {
        pi_l[,,i]=seq(from=(i-1)*t+1,to=i*t,1)
      }

      for (j in 1:k) {
        for (i in 1:m) {
          Xbest[,j][Xbest[,j]==i]=sample(pi_l[,,i])*100
        }
      }
      Xbest=Xbest/100
    }

    #step 2 on page 481 ends


    #Stage II starts

    if (stage2==TRUE){
      C=1

      X=Xbest

      TP=T0;Flag=1

      while (C<=N) {

        time0=Sys.time()

        while(Flag==1 & TP>Tmin){
          Flag=0;I=1

          while (I<=Imax) {
            z=stats::runif(1,0,1)     #step i

            #step ii
            if (z<=0){

              rt=sample(1:t,1)        #randomly select a slice from X

              rcol=sample(1:k,1)      #randomly choose a column

              rrow=sample(seq(from=rt*m,by=-1,length.out=m),2)    #randomly choose 2 rows from slice rt

              Xnew=X

              e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
              e2=Xnew[rrow[2],rcol]

              Xnew[rrow[1],rcol]=e2
              Xnew[rrow[2],rcol]=e1
            }

            #step iii
            if (z>0){

              rcol=sample(1:k,1)      #randomly choose a column

              if(t>1){

                rl=sample(1:m,1)        #randomly choose a l, where l=1, ..., m

                re=sample(pi_l[,,rl],2)    #randomly choose 2 elements from pi_rl that will be exchanged

                Xnew=X

                rrow=c(which(Xnew[,rcol]==re[1]),which(Xnew[,rcol]==re[2]))

                e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
                e2=Xnew[rrow[2],rcol]

                Xnew[rrow[1],rcol]=e2
                Xnew[rrow[2],rcol]=e1

              }

              if(t==1){
                rrow=sample(1:n,2)       #if there is only one slice

                Xnew=X

                e1=Xnew[rrow[1],rcol]            #exchange 2 elements to form Xnew
                e2=Xnew[rrow[2],rcol]

                Xnew[rrow[1],rcol]=e2
                Xnew[rrow[2],rcol]=e1
              }


            }

            #step iv

            a=MaxProCriterion(X=Xnew)
            b=MaxProCriterion(X=X)
            if (a<b){X=Xnew;Flag=1}
            if (a>=b){
              prob=exp((b-a)/TP)
              draw=sample(c(0,1),1,prob=c(1-prob,prob))    #draw==1 means replace
              if(draw==1){X=Xnew;Flag=1}
            }                         #step 5 ends here

            c=MaxProCriterion(X=Xbest)
            if (a<c){Xbest=Xnew;I=1}
            if (a>=c){I=I+1}

          }

          TP=TP*(1-rate)
        }

        time1=Sys.time()
        timediff=time1-time0
        timeALL2=c(timeALL2,timediff)

        ##########progress bar codes
        utils::setTxtProgressBar(progressbar, C)
        ##########


        timeALL=sum(timeALL1)+sum(timeALL2)

        if(as.numeric(timeALL+timediff)<=maxtime){C=C+1}
        if(as.numeric(timeALL+timediff)>maxtime){C=N+1}

        TP=T0;Flag=1

      }

    }

    #Stage II ends
  }

  avgtime1=round(mean(timeALL1),2)
  iterations1=length(timeALL1)

  if(stage2==TRUE){
  avgtime2=round(mean(timeALL2),2)
  iterations2=length(timeALL2)
  }

  close(progressbar)
  print(paste0("average CPU time per iteration for Stage I is: ", avgtime1, " seconds"))
  if(stage2==TRUE){print(paste0("average CPU time per iteration for Stage II is: ", avgtime2, " seconds"))}
  print(paste0("the number of iterations completed for Stage I is: ", iterations1))
  if(stage2==TRUE){print(paste0("the number of iterations completed for Stage II is: ", iterations2))}
  print(paste0("the elements in design matrix is scaled to be 1 to ", n))

  Xbest
}