## Erratum for:
## De Leeuw J, Mair P (2009). Multidimensional Scaling Using Majorization:
## SMACOF in R. Journal of Statistical Software, 31(3), 1-30.
## doi:10.18637/jss.v031.i03

## The original replication code (v31i09.R) is not compatible anymore with
## recent versions of "smacof". Also the URL for the ciexyz31_1.txt data
## is not available anymore. Hence a copy is provided in the supplementary
## materials now.


library("smacof")
packageVersion("smacof")   ## 1.8.15 or higher
data("ekman")
ekman.d <- sim2diss(ekman, method = 1)
res.basic <- mds(ekman.d, type = "ordinal")


###################################################
res.sphere <- smacofSphere(ekman.d, algorithm = "primal", type = "ordinal")


###################################################
plot(res.basic, main = "Configuration Basic SMACOF", asp = 1)
plot(res.sphere, main = "Configuration Sphere SMACOF", asp = 1)


###################################################
wl <- c(434, 445, 465, 472, 490, 504, 537, 555, 584, 600, 610, 628, 651, 674)
CIE31 <- read.table("ciexyz31_1.txt", header = FALSE, sep = ",", colClasses = rep("numeric", 4), nrows = 471)
names(CIE31) <- c("WL", "xbar", "ybar", "zbar")
ekmanWL <- subset(CIE31, WL %in% wl)
ekmanWL[, -1] <- ekmanWL[, -1]/ekmanWL[, "ybar"]


###################################################
ekman.mds.wld <- res.basic$confdiss[cumsum(c(1, 13:2))]
ekmanLuv <- convertColor(ekmanWL[, -1], "XYZ", "Luv")
Luv.dist <- dist(ekmanLuv)
Luv.wld <- Luv.dist[cumsum(c(1, 13:2))]


###################################################
plot(wl[-14], ekman.mds.wld, type = "b", ylim = c(0.05, 0.9), xlab = "wavelength", ylab = "distance", main = "Ekman MDS vs. Luv")
lines(wl[-14], Luv.wld/400, col = "red", type = "b")
legend(600, 0.9, c("Ekman", "Luv"), lty = 1, col = c("black", "red"))


###################################################
data("breakfast")
res.rect <- unfolding(breakfast)


###################################################
plot(res.rect, joint = TRUE, xlim = c(-10, 10), asp = 1)


###################################################
data("bread")
res.uc <- smacofIndDiff(bread)
res.uc


###################################################
res.id <- smacofIndDiff(bread, constraint = "identity")
res.id


###################################################
plot(res.uc, main = "Group Configurations Unconstrained", xlim = c(-1.2,1), ylim = c(-1.2,1), asp = 1)
plot(res.id, main = "Group Configurations Identity", xlim = c(-1.2,1), ylim = c(-1.2,1), asp = 1)


###################################################
data("kinshipdelta")
data("kinshipscales")
res.sym <- mds(kinshipdelta)
res.lin <- smacofConstraint(kinshipdelta, constraint = "linear", external = kinshipscales)


###################################################
plot(res.sym, main = "Configuration Plot SMACOF", xlim = c(-1, 1), ylim = c(-1, 1.5), asp = 1)
plot(res.lin, main = "Configuration Plot SMACOF Constraint",xlim = c(-1, 1), ylim = c(-1, 1.5), asp = 1)
plot(res.sym, plot.type = "stressplot", main = "Stress Decomposition SMACOF", xlim = c(0, 16), ylim = c(3, 12), asp = 1)
plot(res.lin, plot.type = "stressplot", main = "Stress Decomposition SMACOF Constraint", xlim = c(0, 16), ylim = c(3, 12), asp = 1)


###################################################
data("trading")
res.sphere <- smacofSphere(trading, ndim = 3)