smacofBelow <- function (delta, p = 2, xini = torgerson (symmetricFromTriangle (delta, diagonal = FALSE), p), w = rep (1, length (delta)), bnd = delta, itmax = 1000, eps = 1e-10, verbose = FALSE) { m <- length (delta) n <- (1 + sqrt (1 + 8 * m)) / 2 wnd <- which (bnd < Inf) l <- length (wnd) r <- p * (n - 1) yold <- rep (0, l) 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) itel <- 1 sold <- Inf ssq <- sum (w * delta ^ 2) 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))) bmat <- bmat + w[k] * (delta[k] / d[k]) * ak } xmid <- as.vector (bmat %*% xx) ccomp <- c(ssq,-bnd[wnd] ^ 2) / 2 bcomp <- matrix (0, r, l + 1) bcomp[, 1] <- -xmid acomp <- array (0, c(r, r, l + 1)) acomp[, , 1] <- diag (r) t <- 1 for (k in wnd) { ak <- directSum (repList (h\$a[, , k], p)) acomp[, , t + 1] <- ak t <- t + 1 } v <- qpqc (yold, acomp, bcomp, ccomp, verbose = FALSE) snew <- 2 * v\$f / ssq xnew <- v\$xmin ynew <- v\$multipliers cons <- v\$constraints 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 yold <- ynew 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, multipliers = ynew, constraints = cons, itel = itel, stress = sum (w * (delta - d) ^ 2) / ssq ) ) }