torgerson <- function (delta, p = 2) { z <- slanczos(-doubleCenter((delta ^ 2) / 2), p) w <- matrix (0, p, p) v <- pmax(z\$values, 0) diag (w) <- sqrt (v) return(z\$vectors %*% w) } basisPrep <- function (n, w = rep (1, n * (n - 1) / 2)) { m <- n * (n - 1) / 2 v <- -symmetricFromTriangle (w, diagonal = FALSE) diag (v) <- -rowSums(v) ev <- eigen (v) eval <- ev\$values[1:(n - 1)] evec <- ev\$vectors[, 1:(n - 1)] z <- evec %*% diag (1 / sqrt (eval)) a <- array (0, c(n - 1, n - 1, m)) k <- 1 for (j in 1:(n-1)) { for (i in (j+1):n) { dif <- z[i,] - z[j,] a [, , k] <- outer (dif, dif) k <- k + 1 } } return (list (z = z, a = a)) } configurationInBasis <- function (x, z) { n <- nrow (x) p <- ncol (x) r <- p * (n - 1) y <- rep (0, r) for (s in 1:p) { k <- (s - 1) * (n - 1) + 1:(n - 1) y[k] <- crossprod (z, x[, s]) / diag (crossprod (z)) } return (y) } columnCenter <- function (x) { return (apply (x, 2, function (z) z - mean (z))) } doubleCenter <- function (x) { n <- nrow (x) j <- diag(n) - (1 / n) return (j %*% x %*% j) } squareDistances <- function (x) { d <- diag (x) return (outer (d, d, "+") - 2 * x) }