ei <- function (i, n) { return (ifelse(i == (1:n), 1, 0)) } aij <- function (i, j, n) { u <- ei(i, n) - ei (j, n) return (outer (u, u)) } kdelta <- function (i, j) { ifelse (i == j, 1 , 0) } mprint <- function (x, d = 10, w = 15) { print (noquote (formatC ( x, di = d, wi = w, fo = "f" ))) } mnorm <- function (x) { return (sqrt (sum (x ^ 2))) } anorm <- function (x) { return (max (abs (x))) } basis <- function (i, j, n) { s <- sqrt (2) / 2 a <- matrix (0, n, n) if (i == j) a[i, i] <- 1 else { a[i, j] <- a[j, i] <- s } return (a) } center <- 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) } squareDist <- function (x) { d <- diag (x) return (outer (d, d, "+") - 2 * x) } lowerTrapezoidal <- function (x) { n <- nrow (x) p <- ncol (x) if (p == 1) return (x) for (i in 1:(p - 1)) for (j in (i + 1):p) { a <- diag (p) y <- x[i, c(i, j)] y <- y / sqrt (sum (y ^ 2)) a[i, c (i, j)] <- c(1, -1) * y a[j, c (j, i)] <- y x <- x %*% a } return (x) } symmetricFromTriangle <- function (x, lower = TRUE, diagonal = TRUE) { k <- length (x) if (diagonal) n <- (sqrt (1 + 8 * k) - 1) / 2 else n <- (sqrt (1 + 8 * k) + 1) / 2 if (n != as.integer (n)) stop ("input error") nn <- 1:n if (diagonal && lower) m <- outer (nn, nn, ">=") if (diagonal && (!lower)) m <- outer (nn, nn, "<=") if ((!diagonal) && lower) m <- outer (nn, nn, ">") if ((!diagonal) && (!lower)) m <- outer (nn, nn, "<") b <- matrix (0, n, n) b[m] <- x b <- b + t(b) if (diagonal) diag (b) <- diag(b) / 2 return (b) } triangleFromSymmetric <- function (x, lower = TRUE, diagonal = TRUE) { n <- ncol (x) nn <- 1:n if (diagonal && lower) m <- outer (nn, nn, ">=") if (diagonal && (!lower)) m <- outer (nn, nn, "<=") if ((!diagonal) && lower) m <- outer (nn, nn, ">") if ((!diagonal) && (!lower)) m <- outer (nn, nn, "<") return (x[m]) } columnBasis <- function (n) { x <- matrix (0, n, n - 1) x[outer (1:n, 1:(n - 1), "<=")] <- -1 x[outer (1:n, 1:(n - 1), function (i, j) i - j == 1)] <- 1:(n - 1) return (apply (x, 2, function (y) y / sqrt (sum (y ^ 2)))) } matrixBasis <- function (n, p) { x <- matrix (0, n * p, p * (n - 1)) for (j in 1:p) { x [((j - 1) * n) + (1:n), ((j - 1) * (n - 1)) + (1:(n - 1))] <- columnBasis (n) } return (x) }