smacofShepard62 <- function (delta, p = 2, xini = torgerson (symmetricFromTriangle (delta, diagonal = FALSE), p), top = 100, itmax = 1000, eps = 1e-10, verbose = FALSE) { sdelta <- sort (delta) m <- length (delta) n <- (1 + sqrt (1 + 8 * m)) / 2 r <- p * (n - 1) h <- basisPrep (n, p, rep (1, m)) xold <- rep (0, r) for (s in 1:p) { k <- (s - 1) * (n - 1) + 1:(n - 1) xold[k] <- crossprod (h\$z, xini[, s]) / diag (crossprod (h\$z)) } d <- rep (0, m) itel <- 1 repeat { bmat <- matrix (0, n - 1, n - 1) xx <- matrix (xold, n - 1, p) for (k in 1:m) { ak <- h\$a[, , k] d[k] <- sqrt (sum (xx * (ak %*% xx))) } sumd <- sum (d) dhat <- sdelta[rank (d)] sold <- (sum (dhat * d) - sum (delta * d)) / sumd for (k in 1:m) { ak <- h\$a[, , k] bmat <- bmat + ((dhat [k] - delta[k] - sold) / d[k]) * ak } bmat <- bmat / sumd grad <- as.vector (bmat %*% xx) hstp <- optimize (gfunc, c(0, top), x = xold, grad = grad, delta = delta, sdelta = sdelta, a = h\$a) snew <- hstp\$objective xnew <- xold - hstp\$minimum * grad if (verbose) cat( "Iteration: ", formatC (itel, width = 3, format = "d"), "sold: ", formatC ( sold, digits = 8, width = 12, format = "f" ), "snew: ", formatC ( snew, digits = 8, width = 12, format = "f" ), "stepsize: ", formatC ( hstp\$minimum, digits = 8, width = 12, format = "f" ), "\n" ) if ((itel == itmax) || ((sold - snew) < eps)) break xold <- drop (xnew) itel <- itel + 1 } xconf <- matrix (0, n, p) for (s in 1:p) { k <- (s - 1) * (n - 1) + 1:(n - 1) xconf[, s] <- h\$z %*% xnew[k] } return (list ( delta = delta, dist = d, x = xconf, xvec = xnew, grad = grad, itel = itel, stress = snew )) } gfunc <- function (y, x, grad, delta, sdelta, a) { m <- dim (a)[3] n <- dim (a)[1] p <- length (x) / n z <- x - y * grad xx <- matrix (z, n, p) d <- rep (0, m) for (k in 1:m) { ak <- a[, , k] d[k] <- sqrt (sum (xx * (ak %*% xx))) } dhat <- sdelta[rank (d)] val <- (sum (dhat * d) - sum (delta * d)) / sum (d) return (val) }