R/GA.R

In LHD: Latin Hypercube Designs (LHDs)

#' Genetic Algorithm for LHD
#'
#' \code{GA} returns a \code{n} by \code{k} LHD matrix generated by genetic algorithm (GA)
#'
#' @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 m A positive even integer, which stands for the population size and it must be an even number. The default is set to be 10. A large value of \code{m} will result a high CPU time, and it is recommended to be no greater than 100.
#' @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 pmut A probability, which stands for the probability of mutation. The default is set to be 1/(\code{k} - 1).
#' @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 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.
#'
#' @references Liefvendahl, M., and Stocki, R. (2006) A study on algorithms for optimization of Latin hypercubes. \emph{Journal of Statistical Planning and Inference}, \strong{136}, 3231-3247.
#'
#' @examples
#' #generate a 5 by 3 maximin distance LHD with the default setting
#' try=GA(n=5,k=3)
#' try
#' phi_p(try)   #calculate the phi_p of "try".
#'
#' #Another example
#' #generate a 8 by 4 nearly orthogonal LHD
#' try2=GA(n=8,k=4,OC="AvgAbsCor")
#' try2
#' AvgAbsCor(try2)  #calculate the average absolute correlation.
#' @export

GA=function(n,k,m=10,N=10,pmut=1/(k-1),OC="phi_p",p=15,q=1,maxtime=5){
  #n and k are the rs and fa.
  #m: the number of population and it must be an even number.
  #N: maximum number of iterations.
  #pmut: the probability of mutation
  #OC: optimality criterion, the default is "phi_p", along with default p and q

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

  C=1  #Initialize counter index

  #step 2 starts, each X[,,i] is the L_i, i=1, ..., m
  X=rep(0,n*k*m)

  dim(X)=c(n,k,m)

  for (i in 1:m) {
    X[,,i]=rLHD(n=n,k=k)
  }
  #step 2 ends

  if(OC=="phi_p"){
    #step 3 starts
    result=rep(0,m)

    for (i in 1:m) {
      result[i]=phi_p(X[,,i],p=p,q=q)
    }
    #step 3 ends

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

    while (C<=N) {

      time0=Sys.time()

      #step 4 starts: Select survivors (or parients)

      temp=cbind(result,1:m)

      temp=temp[order(temp[,1]),]

      SI=temp[1:(m/2),2]                 #SI:Survivors Index

      Xnew=rep(0,n*k*(m/2))              #Creat survivors matrix L^{s}

      dim(Xnew)=c(n,k,m/2)

      for (i in 1:(m/2)) {

        Xnew[,,i]=X[,,SI[i]]

      }

      #step 4 ends. Xnew are the survivors matrix

      Xbest=Xnew[,,1]     #step 5 Find the best survivor

      for (i in 2:(m/2)) {    #step 6
        #step 7 and 9 starts
        rcol=sample(1:k,1)

        X[,,i]=Xbest
        X[,rcol,i]=Xnew[,rcol,i]

      } #step 8: end for

      X[,,1]=Xbest   #step 9

      for (i in 2:(m/2)) {    #step 10
        #step 11 and 13 starts
        rcol=sample(1:k,1)

        X[,,(i+(m/2))]=Xnew[,,i]
        X[,rcol,(i+(m/2))]=Xbest[,rcol]

      } #step 12: end for

      X[,,(1+(m/2))]=Xbest   #step 13

      #Until here, the X has been fully updated.

      for (i in 2:m) {      #step 14
        for (j in 1:k) {    #step 15

          z=stats::runif(1,0,1)             #step 16

          if (z<pmut){

            X[,,i]=exchange(X=X[,,i],j=j)   #step 17

          } #step 18: end if

        }   #step 19: end for
      }     #step 20: end for

      #step 21: update phi_p for all the L_i
      for (i in 1:m) {
        result[i]=phi_p(X[,,i],p=p,q=q)
      }

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

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

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

    temp=cbind(result,1:m)

    temp=temp[order(temp[,1]),]

    Xbest=X[,,temp[1,2]]

  }

  if(OC=="AvgAbsCor"){
    #step 3 starts
    result=rep(0,m)

    for (i in 1:m) {
      result[i]=AvgAbsCor(X[,,i])
    }
    #step 3 ends

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

    while (C<=N) {

      time0=Sys.time()

      #step 4 starts: Select survivors (or parients)

      temp=cbind(result,1:m)

      temp=temp[order(temp[,1]),]

      SI=temp[1:(m/2),2]                 #SI:Survivors Index

      Xnew=rep(0,n*k*(m/2))              #Creat survivors matrix L^{s}

      dim(Xnew)=c(n,k,m/2)

      for (i in 1:(m/2)) {

        Xnew[,,i]=X[,,SI[i]]

      }

      #step 4 ends. Xnew are the survivors matrix

      Xbest=Xnew[,,1]     #step 5 Find the best survivor

      for (i in 2:(m/2)) {    #step 6
        #step 7 and 9 starts
        rcol=sample(1:k,1)

        X[,,i]=Xbest
        X[,rcol,i]=Xnew[,rcol,i]

      } #step 8: end for

      X[,,1]=Xbest   #step 9

      for (i in 2:(m/2)) {    #step 10
        #step 11 and 13 starts
        rcol=sample(1:k,1)

        X[,,(i+(m/2))]=Xnew[,,i]
        X[,rcol,(i+(m/2))]=Xbest[,rcol]

      } #step 12: end for

      X[,,(1+(m/2))]=Xbest   #step 13

      #Until here, the X has been fully updated.

      for (i in 2:m) {      #step 14
        for (j in 1:k) {    #step 15

          z=stats::runif(1,0,1)             #step 16

          if (z<pmut){

            X[,,i]=exchange(X=X[,,i],j=j)   #step 17

          } #step 18: end if

        }   #step 19: end for
      }     #step 20: end for

      #step 21: update AvgAbsCor for all the L_i
      for (i in 1:m) {
        result[i]=AvgAbsCor(X[,,i])
      }

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

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

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

    temp=cbind(result,1:m)

    temp=temp[order(temp[,1]),]

    Xbest=X[,,temp[1,2]]

  }

  if(OC=="MaxAbsCor"){
    #step 3 starts
    result=rep(0,m)

    for (i in 1:m) {
      result[i]=MaxAbsCor(X[,,i])
    }
    #step 3 ends

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

    while (C<=N) {

      time0=Sys.time()

      #step 4 starts: Select survivors (or parients)

      temp=cbind(result,1:m)

      temp=temp[order(temp[,1]),]

      SI=temp[1:(m/2),2]                 #SI:Survivors Index

      Xnew=rep(0,n*k*(m/2))              #Creat survivors matrix L^{s}

      dim(Xnew)=c(n,k,m/2)

      for (i in 1:(m/2)) {

        Xnew[,,i]=X[,,SI[i]]

      }

      #step 4 ends. Xnew are the survivors matrix

      Xbest=Xnew[,,1]     #step 5 Find the best survivor

      for (i in 2:(m/2)) {    #step 6
        #step 7 and 9 starts
        rcol=sample(1:k,1)

        X[,,i]=Xbest
        X[,rcol,i]=Xnew[,rcol,i]

      } #step 8: end for

      X[,,1]=Xbest   #step 9

      for (i in 2:(m/2)) {    #step 10
        #step 11 and 13 starts
        rcol=sample(1:k,1)

        X[,,(i+(m/2))]=Xnew[,,i]
        X[,rcol,(i+(m/2))]=Xbest[,rcol]

      } #step 12: end for

      X[,,(1+(m/2))]=Xbest   #step 13

      #Until here, the X has been fully updated.

      for (i in 2:m) {      #step 14
        for (j in 1:k) {    #step 15

          z=stats::runif(1,0,1)             #step 16

          if (z<pmut){

            X[,,i]=exchange(X=X[,,i],j=j)   #step 17

          } #step 18: end if

        }   #step 19: end for
      }     #step 20: end for

      #step 21: update MaxAbsCor for all the L_i
      for (i in 1:m) {
        result[i]=MaxAbsCor(X[,,i])
      }

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

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

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

    temp=cbind(result,1:m)

    temp=temp[order(temp[,1]),]

    Xbest=X[,,temp[1,2]]

  }

  if(OC=="MaxProCriterion"){
    #step 3 starts
    result=rep(0,m)

    for (i in 1:m) {
      result[i]=MaxProCriterion(X[,,i])
    }
    #step 3 ends

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

    while (C<=N) {

      time0=Sys.time()

      #step 4 starts: Select survivors (or parients)

      temp=cbind(result,1:m)

      temp=temp[order(temp[,1]),]

      SI=temp[1:(m/2),2]                 #SI:Survivors Index

      Xnew=rep(0,n*k*(m/2))              #Creat survivors matrix L^{s}

      dim(Xnew)=c(n,k,m/2)

      for (i in 1:(m/2)) {

        Xnew[,,i]=X[,,SI[i]]

      }

      #step 4 ends. Xnew are the survivors matrix

      Xbest=Xnew[,,1]     #step 5 Find the best survivor

      for (i in 2:(m/2)) {    #step 6
        #step 7 and 9 starts
        rcol=sample(1:k,1)

        X[,,i]=Xbest
        X[,rcol,i]=Xnew[,rcol,i]

      } #step 8: end for

      X[,,1]=Xbest   #step 9

      for (i in 2:(m/2)) {    #step 10
        #step 11 and 13 starts
        rcol=sample(1:k,1)

        X[,,(i+(m/2))]=Xnew[,,i]
        X[,rcol,(i+(m/2))]=Xbest[,rcol]

      } #step 12: end for

      X[,,(1+(m/2))]=Xbest   #step 13

      #Until here, the X has been fully updated.

      for (i in 2:m) {      #step 14
        for (j in 1:k) {    #step 15

          z=stats::runif(1,0,1)             #step 16

          if (z<pmut){

            X[,,i]=exchange(X=X[,,i],j=j)   #step 17

          } #step 18: end if

        }   #step 19: end for
      }     #step 20: end for

      #step 21: update MaxProCriterion for all the L_i
      for (i in 1:m) {
        result[i]=MaxProCriterion(X[,,i])
      }

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

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

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

    temp=cbind(result,1:m)

    temp=temp[order(temp[,1]),]

    Xbest=X[,,temp[1,2]]

  }

  avgtime=round(mean(timeALL),2)
  iterations=length(timeALL)

  close(progressbar)
  print(paste0("average CPU time per iteration is: ", avgtime, " seconds"))
  print(paste0("the number of iterations completed is: ", iterations))
  print(paste0("the elements in design matrix is scaled to be 1 to ", n))

  Xbest
}