smacofLossRC <- function (d, w, delta) { m <- length (d) h <- .C( "smacofLossC", as.double (d), as.double (w), as.double (delta), as.integer (m), stress = as.double (0) ) return (h\$stress) } smacofDistRC <- function (x, p) { n <- round (length (x) / p) m <- n * (n - 1) / 2 h <- .C("smacofDistC", as.double (x), as.integer (n), as.integer(p), dist = as.double (rep(0, m))) return (h\$dist) } smacofInitialRC <- function (delta, p) { m <- length (delta) n <- round ((1 + sqrt (1 + 8 * m)) / 2) r <- n * (n + 1) / 2 h <- .C( "smacofInitialC", as.double(delta), as.integer(n), as.integer(p), as.double (rep(0, n)), as.double (rep(0, r)), as.double (rep(0, n * n)), as.double (rep (0, n)), x = as.double (rep (0, n * p)) ) return (h\$x) } smacofBmatRC <- function (d, w, delta) { m <- length (w) n <- round ((1 + sqrt (1 + 8 * m)) / 2) r <- n * (n + 1) / 2 h <- .C( "smacofBmatC", as.double(d), as.double(w), as.double(delta), as.integer (n), bmat = as.double(rep(0, r)) ) return (h\$bmat) } smacofVmatRC <- function (w) { m <- length (w) n <- round ((1 + sqrt (1 + 8 * m)) / 2) r <- n * (n + 1) / 2 h <- .C("smacofVmatC", as.double(w), as.integer (n), vmat = as.double(rep(0, r))) return (h\$vmat) } smacofGuttmanRC <- function (x, b, vinv) { m <- length (b) n <- round ((sqrt (1 + 8 * m) - 1) / 2) r <- length (x) p <- round (r / n) h <- .C( "smacofGuttmanC", as.double(x), as.double (b), as.double (vinv), as.integer(n), as.integer(p), as.double (rep(0, r)), y = as.double (rep(0, r)) ) return (h\$y) } smacofGradientRC <- function (x, b, v) { m <- length (b) n <- round ((sqrt (1 + 8 * m) - 1) / 2) r <- length (x) p <- round (r / n) h <- .C( "smacofGradientC", as.double(x), as.double (b), as.double (v), as.integer(n), as.integer(p), y = as.double (rep(0, r)) ) return (h\$y) } smacofHmatRC <- function (x, b, v, d, w, delta) { m <- length (w) n <- round ((1 + sqrt (1 + 8 * m)) / 2) q <- n * (n + 1) / 2 r <- length (x) p <- round (r / n) h <- .C( "smacofHmatC", as.double(x), as.double (b), as.double (v), as.double (w), as.double (delta), as.double (d), as.integer(n), as.integer(p), as.double (rep (0, q)), hmat = as.double (rep (0, r * (r + 1) / 2)) ) return (h\$hmat) } smacofHessianRC <- function (x, b, v, d, w, delta) { m <- length (w) n <- round ((1 + sqrt (1 + 8 * m)) / 2) q <- n * (n + 1) / 2 r <- length (x) p <- round (r / n) h <- .C( "smacofHessianC", as.double(x), as.double (b), as.double (v), as.double (w), as.double (delta), as.double (d), as.integer(n), as.integer(p), as.double (rep (0, q)), hmat = as.double (rep (0, r * (r + 1) / 2)) ) return (h\$hmat) } smacofUpdateRC <- function (x, w, delta, b, vinv) { m <- length (b) n <- round ((sqrt (1 + 8 * m) - 1) / 2) r <- length (x) p <- round (r / n) q <- n * (n - 1) / 2 h <- .C( "smacofUpdateC", as.double (x), as.double (w), as.double (delta), as.double (vinv), as.integer (n), as.integer (p), dnew = as.double (rep(0, q)), bnew = as.double (b), work = as.double (rep (0, n * p)), snew = as.double (0), xnew = as.double (rep (0, n * p)) ) return (h) } smacofRC <- function (w, delta, p, xold = smacofInitialRC(delta, p), itmax = 1000, eps = 1e-10, verbose = FALSE) { itel <- 1 dold <- smacofDistRC (xold, p) sold <- smacofLossRC (dold, w, delta) bold <- smacofBmatRC (dold, w, delta) vinv <- mpowerRC (smacofVmatRC (w), -1) repeat { xnew <- smacofGuttmanRC (xold, bold, vinv) eiff <- max (abs (xold - xnew)) dnew <- smacofDistRC (xnew, p) bnew <- smacofBmatRC (dnew, w, delta) snew <- smacofLossRC (dnew, w, delta) if (verbose) { cat( "itel ", formatC(itel, width = 4, format = "d"), "eiff ", formatC( eiff, width = 15, digits = 10, format = "f" ), "sold ", formatC( sold, width = 15, digits = 10, format = "f" ), "snew ", formatC( snew, width = 15, digits = 10, format = "f" ), "\n" ) } if ((eiff < eps) || (itel == itmax)) { break } itel <- itel + 1 xold <- xnew bold <- bnew dold <- dnew sold <- snew } return (list ( x = xnew, d = dnew, b = bnew, s = snew, itel = itel )) } smacofRCU <- function (w, delta, p, xold = smacofInitialRC(delta, p), itmax = 1000, eps = 1e-10, verbose = FALSE) { itel <- 1 dold <- smacofDistRC (xold, p) sold <- smacofLossRC (dold, w, delta) bold <- smacofBmatRC (dold, w, delta) vinv <- mpowerRC (smacofVmatRC (w), -1) repeat { h <- smacofUpdateRC(xold, w, delta, bold, vinv) xnew <- h\$xnew dnew <- h\$dnew snew <- h\$snew bnew <- h\$bnew eiff <- max (abs (xold - xnew)) if (verbose) { cat( "itel ", formatC(itel, width = 4, format = "d"), "eiff ", formatC( eiff, width = 15, digits = 10, format = "f" ), "sold ", formatC( sold, width = 15, digits = 10, format = "f" ), "snew ", formatC( snew, width = 15, digits = 10, format = "f" ), "\n" ) } if ((eiff < eps) || (itel == itmax)) { break } itel <- itel + 1 xold <- xnew bold <- bnew dold <- dnew sold <- snew } return (list ( x = xnew, d = dnew, b = bnew, s = snew, itel = itel )) } smacofEllipseRC <- function (s, t, w, delta, p, eps, relative = FALSE) { m <- length (w) n <- round ((1 + sqrt (1 + 8 * m)) / 2) r <- n * p hsma <- smacofRCU(w, delta, p) eps <- ifelse (relative, 2 * hsma\$s * eps, 2 * eps) hvec <- smacofHessianRC(hsma\$x, hsma\$b, smacofVmatRC (w), hsma\$d, w, delta) xmat <- matrix (hsma\$x, n, p) plot ( xmat, col = "RED", pch = 21, bg = "RED", xlab = paste("dim ", formatC( s, digits = 1, format = "d" )), ylab = paste("dim ", formatC( t, digits = 1, format = "d" )) ) hmat <- matrix (trimatRC(hvec), r, r) z <- cbind (sin (seq(-pi, pi, length = 100)), cos (seq(-pi, pi, length = 100))) * sqrt(eps) for (i in 1:n) { ii <- c((s - 1) * n + i, (t - 1) * n + i) y <- xmat[i, c(s, t)] amat <- c(hmat[ii[1], ii[1]], hmat[ii[2], ii[1]], hmat[ii[2], ii[2]]) hjac <- jacobi (amat) zi <- z %*% diag (1 / sqrt (hjac\$values)) zi <- tcrossprod(zi, hjac\$vectors) zi <- zi + matrix (y, 100, 2, byrow = TRUE) lines (list(x = zi[, 1], y = zi[, 2])) } return ( list ( x = xmat, dist = hsma\$d, bmat = hsma\$b, hessian = hmat, stress = hsma\$s, itel = hsma\$itel ) ) }