smacofAbove <- function (delta, p = 2, xini = torgerson (symmetricFromTriangle (delta, diagonal = FALSE), p), w = rep (1, length (delta)), bnd = rep (0, length (delta)), itmax = 1000, eps = 1e-10, verbose = FALSE) { m <- length (delta) n <- (1 + sqrt (1 + 8 * m)) / 2 wnd <- which (bnd > 0) l <- length (wnd) xd <- as.vector (dist (xini)) lbd <- max (bnd[wnd] / xd[bnd]) xini <- lbd * xini r <- p * (n - 1) h <- basisPrep (n, p, w) 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) e <- matrix (0, l, r) itel <- 1 sold <- Inf ssq <- sum (w * delta ^ 2) repeat { xmid <- matrix (0, n - 1, p) xx <- matrix (xold, n - 1, p) t <- 1 for (k in 1:m) { ak <- h\$a[, , k] ax <- ak %*% xx d[k] <- sqrt (sum (xx * ax)) if (!is.na (match (k, wnd))) { e[t, ] <- as.vector (ax / d[k]) t <- t + 1 } xmid <- xmid + w[k] * (delta[k] / d[k]) * ax } xnew <- qp (q = diag (r), p = -as.vector (xmid), e = e, f = bnd[wnd]) snew <- sum (w * (delta - d) ^ 2) / ssq 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" ), "\n" ) if ((itel == itmax) || ((sold - snew) < eps)) break xold <- xnew sold <- snew 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, itel = itel, stress = snew, constraints = d[wnd] - bnd[wnd] ) ) }