R/bdf3.mep.R

In bdf3: Efficient Block Designs for 3-Level Factorial Experiments in Block Size 3

#' @importFrom dplyr %>%
#' @importFrom stats na.omit
#'
bdf3.mep <- function(n_factors = 3, show_efficiency = TRUE) {

  label_effect <- function(vec, factors) {
    label <- ""
    superscripts <- c("0" = "", "1" = "", "2" = "\u00B2")
    for (i in seq_along(vec)) {
      if (vec[i] != 0) {
        label <- paste0(label, factors[i], superscripts[as.character(vec[i])])
      }
    }
    if (label == "") "I" else label
  }

  generate_effects <- function(n) {
    total <- 3^n
    effects <- matrix(NA, nrow = total - 1, ncol = n)
    for (i in 1:(total - 1)) {
      vec <- integer(n); num <- i
      for (j in 1:n) {
        vec[j] <- num %% 3
        num <- num %/% 3
      }
      effects[i, ] <- vec
    }
    effects
  }

  is_confounded <- function(block, effect) {
    res <- (block %*% effect) %% 3
    length(unique(res)) == 1
  }

  normalize_effect <- function(eff) {
    first_nonzero <- which(eff != 0)[1]
    if (length(first_nonzero) == 0) return(eff)
    if (eff[first_nonzero] == 2) eff <- (3 - eff) %% 3
    eff
  }

  classify_effect <- function(eff) {
    nonzeros <- sum(eff != 0)
    if (nonzeros == 1) "Main"
    else if (nonzeros == 2) "2FI"
    else "Other"
  }

  get_principal_blocks <- function(n) {
    all_runs <- expand.grid(replicate(n, 0:2, simplify = FALSE))
    all_runs <- as.matrix(all_runs)
    zero_run <- rep(0, n)
    blocks <- list()
    for (i in 1:nrow(all_runs)) {
      for (j in 1:nrow(all_runs)) {
        if (i == j) next
        r1 <- all_runs[i, ]; r2 <- all_runs[j, ]
        block <- rbind(zero_run, r1, r2)

        valid <- all(apply(block, 2, function(col) {
          all(col == 0) || setequal(sort(unique(col)), 0:2)
        }))
        if (valid) blocks[[length(blocks) + 1]] <- block
      }
    }
    unique_blocks <- unique(lapply(blocks, function(x)
      paste(apply(x, 1, paste, collapse = ""), collapse = "|")))
    lapply(unique_blocks, function(key) {
      rows_str <- unlist(strsplit(key, "\\|"))
      mat <- matrix(as.integer(unlist(strsplit(rows_str, ""))),
                    ncol = n, byrow = TRUE)
      colnames(mat) <- LETTERS[1:n]
      mat
    })
  }

  compute_efficiency_factors <- function(blocks, normalized_effects, effect_labels) {
    num_blocks <- length(blocks)
    estimable_count <- numeric(length(effect_labels))
    for (blk in blocks) {
      confounded <- apply(normalized_effects, 1,
                          function(eff) is_confounded(blk, eff))
      estimable_count <- estimable_count + as.numeric(!confounded)
    }
    eff_factor <- round(estimable_count / num_blocks, 2)
    names(eff_factor) <- effect_labels
    eff_factor
  }

  generate_unique_blocks <- function(principal_block) {
    n_factors <- ncol(principal_block)
    used_treatments_matrix <- matrix(NA, nrow = 0, ncol = n_factors)
    blocks <- list(principal_block)
    used_treatments_matrix <- rbind(used_treatments_matrix, principal_block)
    shifts <- expand.grid(replicate(n_factors, 0:2, simplify = FALSE))
    shifts <- shifts[!apply(shifts, 1, function(x) all(x == 0)), ]
    for (i in 1:nrow(shifts)) {
      shift <- as.numeric(shifts[i, ])
      derived_block <- (principal_block +
                          matrix(rep(shift, each = nrow(principal_block)),
                                 nrow = nrow(principal_block), byrow = FALSE)) %% 3
      already_used <- apply(derived_block, 1, function(row)
        any(apply(used_treatments_matrix, 1,
                  function(used_row) all(used_row == row))))
      if (!any(already_used)) {
        blocks[[length(blocks) + 1]] <- derived_block
        used_treatments_matrix <- rbind(used_treatments_matrix, derived_block)
      }
      if (nrow(used_treatments_matrix) == 3^n_factors) break
    }
    blocks
  }

  factors <- LETTERS[1:n_factors]
  all_blocks <- get_principal_blocks(n_factors)
  all_effects <- generate_effects(n_factors)
  classified <- apply(all_effects, 1, classify_effect)
  target_effects <- all_effects[classified %in% c("Main", "2FI"), , drop = FALSE]
  normalized_effects <- unique(t(apply(target_effects, 1, normalize_effect)))
  valid_rows <- apply(normalized_effects, 1, function(eff) {
    first <- which(eff != 0)[1]
    !is.na(first) && eff[first] == 1
  })
  normalized_effects <- normalized_effects[valid_rows, , drop = FALSE]
  effect_labels <- apply(normalized_effects, 1,
                         function(eff) label_effect(eff, factors))

  block_estimable <- lapply(all_blocks, function(blk) {
    apply(normalized_effects, 1, function(eff) !is_confounded(blk, eff))
  })
  covered <- c(); chosen_blocks <- list(); used_indices <- c()
  while (length(covered) < nrow(normalized_effects)) {
    best_block <- NULL; max_new <- 0
    for (i in seq_along(block_estimable)) {
      if (i %in% used_indices) next
      new_cov <- which(block_estimable[[i]])[!(which(block_estimable[[i]]) %in% covered)]
      if (length(new_cov) > max_new) {
        max_new <- length(new_cov); best_block <- i
      }
    }
    if (is.null(best_block)) break
    covered <- union(covered, which(block_estimable[[best_block]]))
    chosen_blocks[[length(chosen_blocks) + 1]] <- all_blocks[[best_block]]
    used_indices <- c(used_indices, best_block)
  }

  eff_factors <- if (show_efficiency)
    compute_efficiency_factors(chosen_blocks, normalized_effects, effect_labels)
  else NULL
  replications <- lapply(chosen_blocks, generate_unique_blocks)

  structure(list(
    n_factors = n_factors,
    factors = factors,
    normalized_effects = normalized_effects,
    effect_labels = effect_labels,
    chosen_principal_blocks = chosen_blocks,
    efficiency_factors = eff_factors,
    replications = replications
  ), class = "bdf3mep")
}
print.bdf3mep <- function(x, ...) {
  is_confounded <- function(block, effect) {
    res <- (block %*% effect) %% 3
    length(unique(res)) == 1
  }
  label_effect <- function(vec, factors) {
    label <- ""
    superscripts <- c("0" = "", "1" = "", "2" = "\u00B2")
    for (i in seq_along(vec)) {
      if (vec[i] != 0) label <- paste0(label, factors[i], superscripts[as.character(vec[i])])
    }
    if (label == "") "I" else label
  }
  classify_effect <- function(eff) {
    nonzeros <- sum(eff != 0)
    if (nonzeros == 1) "Main"
    else if (nonzeros == 2) "2FI"
    else "Other"
  }

  cat("Minimal Principal Blocks Covering All Main Effects and 2FIs\n")
  cat("Factors:", paste(x$factors, collapse = " "), "\n\n")

  for (i in seq_along(x$chosen_principal_blocks)) {
    block <- x$chosen_principal_blocks[[i]]
    colnames(block) <- x$factors
    cat(paste0("Principal Block ", i, ":\n"))
    print(as.data.frame(block), row.names = FALSE)

    confounded <- apply(x$normalized_effects, 1, function(eff) {
      eff <- as.numeric(eff)
      if (is_confounded(block, eff)) label_effect(eff, x$factors) else NA
    })
    classified <- apply(x$normalized_effects, 1, classify_effect)
    confounded_main <- na.omit(confounded[classified == "Main"])
    confounded_2FI  <- na.omit(confounded[classified == "2FI"])

    cat("Confounded Main Effects:",
        if (length(confounded_main) == 0) "None" else paste(confounded_main, collapse = ", "),
        "\n")
    cat("Confounded 2FI Effects:",
        if (length(confounded_2FI) == 0) "None" else paste(confounded_2FI, collapse = ", "),
        "\n\n")
  }

  if (!is.null(x$efficiency_factors)) {
    cat("Efficiency factors:\n")
    print(x$efficiency_factors)
    cat("\n")
  }

  print_blocks_side_by_side <- function(block_list, factors, blocks_per_row = 5) {
    n_blocks <- length(block_list)
    n_runs <- nrow(block_list[[1]])
    col_width <- 14

    block_rows <- split(seq_len(n_blocks), ceiling(seq_len(n_blocks)/blocks_per_row))

    for (row_blocks in block_rows) {
      for (b in row_blocks) cat(format(paste0("Block ", b), width = col_width), "")
      cat("\n")

      for (b in row_blocks) cat(format(paste(factors, collapse = " "), width = col_width), "")
      cat("\n")

      for (r in 1:n_runs) {
        for (b in row_blocks) {
          cat(format(paste(block_list[[b]][r, ], collapse = " "), width = col_width), "")
        }
        cat("\n")
      }
      cat("\n")
    }
  }

  for (i in seq_along(x$replications)) {
    cat(paste0("Replication ", i, ":\n"))
    print_blocks_side_by_side(x$replications[[i]], x$factors)
  }

  invisible(x)
}