node2vec_holdout <- function(A, repeats = 30, holdout_ties = 17,
                             rank = 4, p = 1, q = 0.5,
                             walks_per_node = 30, walk_length = 20,
                             window = 4, seed = 6886) {
  A <- 1L * (A > 0)
  A[is.na(A)] <- 0L
  A <- 1L * ((A + t(A)) > 0)
  diag(A) <- 0L
  n <- nrow(A)
  pairs <- which(upper.tri(A), arr.ind = TRUE)
  labels <- A[pairs]

  pick_one <- function(x, w = NULL) {
    x[sample.int(length(x), 1L, prob = w)]
  }

  make_walks <- function(train_A) {
    walks <- vector("list", n * walks_per_node)
    out_index <- 1L
    for (start in seq_len(n)) {
      for (rep_index in seq_len(walks_per_node)) {
        walk <- start
        while (length(walk) < walk_length) {
          current <- tail(walk, 1L)
          neighbors <- which(train_A[current, ] > 0)
          if (!length(neighbors)) break
          if (length(walk) == 1L) {
            next_node <- pick_one(neighbors)
          } else {
            previous <- walk[length(walk) - 1L]
            weight <- ifelse(
              neighbors == previous, 1 / p,
              ifelse(train_A[previous, neighbors] > 0, 1, 1 / q)
            )
            next_node <- pick_one(neighbors, weight)
          }
          walk <- c(walk, next_node)
        }
        walks[[out_index]] <- walk
        out_index <- out_index + 1L
      }
    }
    walks
  }

  make_coordinates <- function(train_A) {
    walks <- make_walks(train_A)
    counts <- matrix(0, n, n)
    for (walk in walks) {
      for (position in seq_along(walk)) {
        lower <- max(1L, position - window)
        upper <- min(length(walk), position + window)
        context <- walk[seq.int(lower, upper)]
        context <- context[context != walk[position]]
        if (length(context)) {
          for (context_node in context) {
            counts[walk[position], context_node] <-
              counts[walk[position], context_node] + 1
          }
        }
      }
    }
    row_total <- rowSums(counts)
    column_total <- colSums(counts)
    denominator <- outer(row_total, column_total)
    ppmi <- matrix(0, n, n)
    observed <- counts > 0 & denominator > 0
    ppmi[observed] <- pmax(
      log(counts[observed] * sum(counts) / denominator[observed]),
      0
    )
    decomposition <- svd(ppmi, nu = rank, nv = 0)
    sweep(
      decomposition$u[, seq_len(rank), drop = FALSE],
      2,
      sqrt(decomposition$d[seq_len(rank)]),
      "*"
    )
  }

  make_pair_features <- function(train_A, coordinates) {
    source <- pairs[, 1]
    target <- pairs[, 2]
    degree <- rowSums(train_A)
    common_neighbors <- vapply(
      seq_along(source),
      function(index) {
        sum(train_A[source[index], ] * train_A[target[index], ])
      },
      numeric(1)
    )
    neighbor_union <- degree[source] + degree[target] - common_neighbors
    baseline <- data.frame(
      degree_sum = degree[source] + degree[target],
      degree_product = degree[source] * degree[target],
      common_neighbors = common_neighbors,
      jaccard = ifelse(
        neighbor_union > 0,
        common_neighbors / neighbor_union,
        0
      )
    )
    embedding <- data.frame(
      dot = rowSums(
        coordinates[source, , drop = FALSE] *
          coordinates[target, , drop = FALSE]
      )
    )
    for (dimension in seq_len(ncol(coordinates))) {
      embedding[[paste0("absolute_difference_", dimension)]] <-
        abs(
          coordinates[source, dimension] -
            coordinates[target, dimension]
        )
    }
    list(
      baseline = baseline,
      embedding = embedding,
      combined = cbind(baseline, embedding)
    )
  }

  auc <- function(outcome, score) {
    positive <- score[outcome == 1]
    negative <- score[outcome == 0]
    mean(outer(positive, negative, ">")) +
      0.5 * mean(outer(positive, negative, "=="))
  }

  fit_and_predict <- function(features, train, test) {
    training_data <- data.frame(
      outcome = labels[train],
      features[train, , drop = FALSE]
    )
    fit <- suppressWarnings(
      glm(outcome ~ ., data = training_data, family = binomial())
    )
    as.numeric(
      predict(
        fit,
        newdata = features[test, , drop = FALSE],
        type = "response"
      )
    )
  }

  one_split <- function(split) {
    set.seed(seed + split)
    tied <- which(labels == 1)
    untied <- which(labels == 0)
    test_ties <- sample(tied, holdout_ties)
    test_nonties <- sample(untied, holdout_ties)
    test <- c(test_ties, test_nonties)
    train <- setdiff(seq_along(labels), test)
    train_A <- A
    train_A[pairs[test_ties, , drop = FALSE]] <- 0L
    train_A[cbind(pairs[test_ties, 2], pairs[test_ties, 1])] <- 0L
    coordinates <- make_coordinates(train_A)
    features <- make_pair_features(train_A, coordinates)
    c(
      baseline = auc(
        labels[test],
        fit_and_predict(features$baseline, train, test)
      ),
      walk_coordinates = auc(
        labels[test],
        fit_and_predict(features$embedding, train, test)
      ),
      combined = auc(
        labels[test],
        fit_and_predict(features$combined, train, test)
      )
    )
  }

  result <- t(vapply(seq_len(repeats), one_split, numeric(3)))
  data.frame(split = seq_len(repeats), result, row.names = NULL)
}
