library (rcdd) library (numDeriv) library (vertexenum) library (MASS) inverseMDS <- function (x) { n <- nrow (x) m <- ncol (x) x <- apply (x, 2, function (y) y - mean (y)) nm <- n - (m + 1) kk <- cbind (1, x, matrix (rnorm (n * nm), n , nm)) kperp <- as.matrix (qr.Q (qr (kk))[, -(1:(m + 1))]) dd <- Euclid (x) k <- 1 base <- matrix (0, n * (n - 1) / 2, nm * (nm + 1) / 2) for (i in 1:nm) { for (j in 1:i) { oo <- outer (kperp[, i], kperp[, j]) if (j != i) { oo <- oo + t(oo) } base[, k] <- lower_triangle (dd + (1 - oo)) k <- k + 1 print (c(i, j, k)) } } return (base = cbind (lower_triangle (dd), base)) } inversePlus <- function (base, affine = TRUE) { if (affine) { hrep <- makeH ( a1 = d2q (-base), b1 = d2q (rep (0, nrow (base))), a2 = d2q (rep (1, ncol (base))), b2 = d2q (1) ) } else { hrep <- makeH (a1 = d2q (-base), b1 = d2q (rep (0, nrow (base)))) } vrep <- scdd (hrep) hrep <- q2d (hrep) vrep <- q2d (vrep\$output) pr <- tcrossprod (hrep[,-c(1, 2)], vrep[,-c(1, 2)])[-1,] return (list ( base = base, hrep = hrep, vrep = vrep, pr = pr )) } twoPoints <- function (x, y) { dx <- inverseMDS (x) dy <- inverseMDS (y) mx <- ncol (dx) my <- ncol (dy) ex <- eigen (crossprod (cbind (dx, -dy))) iz <- which (abs (ez\$values) < 1e-10) ax <- ez\$vectors[1:mx, iz] ay <- ez\$vectors[mx + (1:my), iz] fx <- dx %*% ax fy <- dy %*% ay return (list (fx = fx, fy = fy)) } second_partials_stress <- function (x, delta, w = nonDiag (nrow (x))) { n <- nrow (x) p <- ncol (x) d <- Euclid (x) fac <- (w * delta) / (d + diag (n)) dd <- d * d v <- vmat (w) deri <- direct_sum (repList (v, p)) xx <- as.vector (x) for (i in 1:(n - 1)) { for (j in (i + 1):n) { aa <- direct_sum (repList (aijn (i, j, n), p)) ax <- drop (aa %*% xx) deri <- deri - fac[i, j] * (aa - outer (ax, ax) / dd[i, j]) } } return (deri) } second_partials_numerical <- function (x, delta, w = nonDiag (nrow (x))) { stress <- function (x, delta, w) { n <- nrow (delta) p <- length (x) / n d <- Euclid (matrix (x, n, p)) res <- delta - d return (sum (w * res * res) / 2) } return (hessian (stress, x, delta = delta, w = w)) } lower_triangle <- function (x) { n <- nrow (x) return (x[outer(1:n, 1:n, ">")]) } fill_symmetric <- function (x) { m <- length (x) n <- (1 + sqrt (1 + 8 * m)) / 2 d <- matrix (0, n, n) d[outer(1:n, 1:n, ">")] <- x return (d + t(d)) } Euclid <- function (x) { c <- tcrossprod (x) d <- diag (c) return (sqrt (abs (outer (d, d, "+") - 2 * c))) } circular <- function (n) { x <- seq (0, 2 * pi, length = n + 1) z <- matrix (0, n + 1, 2) z[, 1] <- sin (x) z[, 2] <- cos (x) return (z[-1,]) } direct_sum <- function (x) { n <- length (x) nr <- sapply (x, nrow) nc <- sapply (x, ncol) s <- matrix (0, sum (nr), sum (nc)) k <- 0 l <- 0 for (j in 1:n) { s[k + (1:nr[j]), l + (1:nc [j])] <- x[[j]] k <- k + nr[j] l <- l + nc[j] } return (s) } aijn <- function (i, j, n) { a <- matrix (0, n, n) a[i, i] <- 1 a[j, j] <- 1 a[i, j] <- -1 a[j, i] <- -1 return (a) } repList <- function (x, n) { z <- list () for (i in 1:n) { z <- c (z, list (x)) } return (z) } nonDiag <- function (n) { return (matrix (1, n, n) - diag (n)) } bmat <- function (x, delta, w = nonDiag (nrow (x))) { d <- Euclid (x) bmat <- -(w * delta) / (d + diag (nrow (x))) diag (bmat) <- -rowSums (bmat) return (bmat) } vmat <- function (w) { vmat <- -w diag (vmat) <- -rowSums (vmat) return (vmat) } torgerson <- function (delta) { n <- nrow (delta) dd <- delta * delta sd <- rowSums (dd) / n ed <- sum (dd) / (n * n) cc <- -(dd - outer (sd, sd, "+") + ed) / 2 return (eigen (dd)) } cleanUp <- function (a, eps = 1e-3) { nv <- nrow (a) ind <- rep (TRUE, nv) for (i in 1:(nv - 1)) { xx <- a[i,] for (j in (i + 1):nv) { if (!ind[j]) next yy <- a[j,] mm <- max (abs (xx - yy)) if (mm < eps) ind[j] <- FALSE } } return (ind) } bruteForce <- function (a, b, eps = 1e-3) { n <- nrow (a) m <- ncol (a) cb <- combn (n, m) nb <- ncol (cb) ind <- rep(TRUE, nb) ht <- numeric() for (i in 1:nb) { gg <- a[cb[, i], ] bg <- b[cb[, i]] qg <- qr(gg) if (qg\$rank < m) { ind[i] <- FALSE next } hh <- solve (qg, bg) hg <- drop (a %*% hh) if (min (b - hg) < -eps) { ind[i] <- FALSE next } ht <- c(ht, hh) } n1 <- ncol (cb) n2 <- sum (ind) ht <- t (matrix (ht, m, n2)) ind <- cleanUp (ht, 1e-3) n3 <- sum (ind) return (list ( x = ht[ind,], n1 = n1, n2 = n2, n3 = n3 )) } rankTest <- function (x, a, b, eps = 1e-3) { h <- drop (a %*% x) ind <- which (abs (h - b) < eps) r <- qr (a[ind,])\$rank f <- min (b - h) > -eps return (list (rank = r, feasibility = f)) } makeDC <- function (x) { y <- -x diag(y) <- -rowSums (y) return (y) } bmat <- function (delta, w, d) { n <- nrow (w) dd <- ifelse (d == 0, 0, 1 / d) return (makeDC (w * delta * dd)) } smacof <- function (delta, w, xini, eps = 1e-6, itmax = 100, pen = 0, verbose = TRUE) { n <- nrow (xini) xold <- xini dold <- as.matrix (dist (xold)) sold <- sum (w * (delta - dold) ^ 2) / 2 itel <- 1 v <- ginv (makeDC (w)) repeat { b <- bmat (delta, w, dold) xnew <- v %*% b %*% xold / (1 + pen) dnew <- as.matrix (dist (xnew)) snew <- sum (w * (delta - dnew) ^ 2) / 2 if (verbose) { cat ( formatC (itel, width = 4, format = "d"), formatC ( sold, digits = 10, width = 13, format = "f" ), formatC ( snew, digits = 10, width = 13, format = "f" ), "\n" ) } if ((itel == itmax) || (sold - snew) < eps) break () itel <- itel + 1 sold <- snew dold <- dnew xold <- xnew } return (list ( x = xnew, d = dnew, s = snew, itel = itel )) }