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

In SFDesign: Space-Filling Designs

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
)

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.

Value

design	optimized LHD.
total.iter	total number of swaps in the optimization.
criterion	optimized criterion.
crit.hist	criterion history during the optimization process.

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)

SFDesign documentation built on June 22, 2025, 1:06 a.m.