SFDesign (version 0.1.3)

customLHD: Generate a Latin-hypercube design (LHD) based on a custom criterion

Description

This function generates a LHD by minimizing a user-specified design criterion.

Usage

customLHD(
  compute.distance.matrix,
  compute.criterion,
  update.distance.matrix,
  n,
  p,
  design = NULL,
  max.sa.iter = 1e+06,
  temp = 0,
  decay = 0.95,
  no.update.iter.max = 400,
  num.passes = 10,
  max.det.iter = 1e+06,
  method = "full",
  scaled = TRUE
)

Value

design

optimized LHD.

total.iter

total number of swaps in the optimization.

criterion

optimized criterion.

crit.hist

criterion history during the optimization process.

Arguments

compute.distance.matrix

a function to calculate pairwise distance

compute.criterion

a function to calculate the criterion based on the pairwise distance

update.distance.matrix

a function to update the distance matrix after swapping one column of two design points

n

design size.

p

design dimension.

design

an initial LHD. If design=NULL, a random LHD is generated.

max.sa.iter

maximum number of swapping involved in the simulated annealing (SA) algorithm.

temp

initial temperature of the simulated annealing algorithm. If temp=0, it will be automatically determined.

decay

the temperature decay rate of simulated annealing.

no.update.iter.max

the maximum number of iterations where there is no update to the global optimum before SA stops.

num.passes

the maximum number of passes of the whole design matrix if deterministic swapping is used.

max.det.iter

maximum number of swapping involved in the deterministic swapping algorithm.

method

choice of "deterministic", "sa", or "full". See details for the description of each choice.

scaled

whether the design is scaled to unit hypercube. If scaled=FALSE, the design is represented by integer numbers from 1 to design size. Leave it as TRUE when no initial design is provided.

Details

customLHD generates a LHD by minimizing a user-specified design criterion.

• If method='sa', then a simulated annealing algorithm is used to optimize the LHD. To custom the optimization process, you can change the default values for max.sa.iter, temp, decay, no.update.iter.max. In this optimization step, two design points are randomly chosen and their coordinate along one dimension are swaped. If the new design improves the criterion, then it is accepted; otherwise, it is accepted with some probability.

• If method='deterministic', then a deterministic swap algorithm is used to optimize the LHD. To custom the optimization process, you can change the default values for num.passes, max.det.iter. In this optimization step, we swap the coordinates of all pairs of design points (start with design point 1 with design point 2, then 1 with 3, ... 1 with n, then 2 with 3 until n-1 with n). Only accept the change if the swap leads to an improvement.

• If method='full', then optimization goes through the above two stages.

Examples

# Below is an example showing how to create functions needed to generate
# MaxPro LHD manually by customLHD without using the maxproLHD function in
# the package.
compute.distance.matrix <- function(A){
   s = 2
   log_prod_metric = function(x, y) s * sum(log(abs(x-y)))
   return (c(proxy::dist(A, log_prod_metric)))
}
compute.criterion <- function(n, p, d) {
  s = 2
  dim <- as.integer(n * (n - 1) / 2)
  # Find the minimum distance
  Dmin <- min(d)
  # Compute the exponential summation
  avgdist <- sum(exp(Dmin - d))
  # Apply the logarithmic transformation and scaling
  avgdist <- log(avgdist) - Dmin
  avgdist <- exp((avgdist - log(dim)) * (p * s) ^ (-1))
  return(avgdist)
}

update.distance.matrix <- function(A, col, selrow1, selrow2, d) {
  s = 2
  n = nrow(A)
  # transform from c++ idx to r idx
  selrow1 = selrow1 + 1
  selrow2 = selrow2 + 1
  col = col + 1
  # A is the updated matrix
  row1 <- min(selrow1, selrow2)
  row2 <- max(selrow1, selrow2)

  compute_position <- function(row, h, n) {
    n*(h-1) - h*(h-1)/2 + row-h
  }

  # Update for rows less than row1
  if (row1 > 1) {
    for (h in 1:(row1-1)) {
      position1 <- compute_position(row1, h, n)
      position2 <- compute_position(row2, h, n)
      d[position1] <- d[position1] + s * log(abs(A[row1, col] - A[h, col])) -
          s * log(abs(A[row2, col] - A[h, col]))
      d[position2] <- d[position2] + s * log(abs(A[row2, col] - A[h, col])) -
          s * log(abs(A[row1, col] - A[h, col]))
    }
  }

  # Update for rows between row1 and row2
  if ((row2-row1) > 1){
    for (h in (row1+1):(row2-1)) {
      position1 <- compute_position(h, row1, n)
      position2 <- compute_position(row2, h, n)
      d[position1] <- d[position1] + s * log(abs(A[row1, col] - A[h, col])) -
          s * log(abs(A[row2, col] - A[h, col]))
      d[position2] <- d[position2] + s * log(abs(A[row2, col] - A[h, col])) -
          s * log(abs(A[row1, col] - A[h, col]))
    }
  }

  # Update for rows greater than row2
  if (row2 < n) {
    for (h in (row2+1):n) {
      position1 <- compute_position(h, row1, n)
      position2 <- compute_position(h, row2, n)
      d[position1] <- d[position1] + s * log(abs(A[row1, col] - A[h, col])) -
          s * log(abs(A[row2, col] - A[h, col]))
      d[position2] <- d[position2] + s * log(abs(A[row2, col] - A[h, col])) -
          s * log(abs(A[row1, col] - A[h, col]))
    }
  }
  return (d)
}

n = 6
p = 2
# Find an appropriate initial temperature
crit1 = 1 / (n-1)
crit2 = (1 / ((n-1)^(p-1) * (n-2))) ^ (1/p)
delta = crit2 - crit1
temp = - delta / log(0.99)
result_custom = customLHD(compute.distance.matrix,
function(d) compute.criterion(n, p, d),
update.distance.matrix, n, p, temp = temp)