#This R code allows to reproduce the results of the paper I. Epifanio, G. Vinué, S. Alemany, Archetypal analysis: Contributions for estimating boundary cases in multivariate accommodation problem, Computers & Industrial Engineering, Volume 64, Issue 3, March 2013, Pages 757-765, ISSN 0360-8352, 10.1016/j.cie.2012.12.011.
(http://www.sciencedirect.com/science/article/pii/S0360835212003063)

#First, we read and preprocess the database USAF 1967 (Zehner et al.).
m <- read.table('med1967.dat', sep = " ")
#Variable selection:
sel <- c(48, 40 ,39, 33, 34, 36)
#Changing to inches: 
mpulg <- m[,sel] / (10 * 2.54)
#Standarizing (normalizing):
smpulg <- scale(mpulg, center = mean(mpulg[,]), scale = sd(mpulg[,])) 

#Next, we remove the more extreme 5% data using the Mahalanobis distance and we check that calculating the depth of each data, we get very similar results.
Sx <- cov(smpulg)
D2 <- mahalanobis(smpulg, colMeans(smpulg), Sx)
#Number of individuals not considered:
sum(D2 > qchisq(0.95, df = dim(mpulg)[2]))
#Individuals that are not considered:
qchi <- which(D2 > qchisq(0.95, df = dim(mpulg)[2]))
library(depth)
dt = c()
for(i in 1 : nrow (smpulg)){ 
 dt[i] <- depth(smpulg[i,], smpulg) 
}
#Number of individuals not considered:
sum(dt == min(dt))
#Individuals that are not considered
qd <- which(dt == min(dt))
#Agreement between Mahalanobis and depth:
intersect(qd, qchi)

#Results shown in Section 3 are obtained as follows.
#Individuals considered for covering 95%:
qlchi <- which(D2 <= qchisq(0.95, df = dim(mpulg)[2])) 
#Database with the 6 variables and the selected individuals:
lsmpulg <- smpulg[qlchi,] 
library("archetypes")
#For reproducing results, seed for randomness:
set.seed(2010) 
#Run archetypes algorithm repeatedly from 1 to 15 archetypes:
lass15 <- stepArchetypes(data = lsmpulg, k = 1:15, verbose = FALSE, nrep = 3) 
#Plot from 1 to 15 archetypes:
screeplot(lass15) 
#3 archetypes:
a3 <- bestModel(lass15[[3]])
parameters(a3)
#7 archetypes:
a7 <- bestModel(lass15[[7]])
parameters(a7) 
#Plotting the percentiles of each archetype:
barplot(a3, smpulg, percentiles = T, which = "beside") 
barplot(a7, smpulg, percentiles = T, which = "beside")

#We get the results in Section 3.3 related with PCA.
pznueva <- prcomp(smpulg, scale = T, retx = T) 
summary(pznueva)
#PCA scores for 3 archetypes:
p3 <- predict(pznueva, parameters(a3)) 
#PCA scores for 7 archetypes:
p7 <- predict(pznueva,parameters(a7)) #pca scores for 7 archetypes
#Representing the scores:
xyplot.pca(p3[,1:2], pznueva$x[,1:2], data.col = gray(0.7), atypes.col = 1, atypes.pch = 15)
xyplot.pca(p7[,1:2], pznueva$x[,1:2], data.col = gray(0.7), atypes.col = 1, atypes.pch = 15)

#At last, we show how the empirical percentiles and the nearest individuals to each one of archetypes can be obtained.
#Function for computing empirical percentiles:
.perc <- function(x, data, digits = 0) {
  Fn <- ecdf(data)
  round(Fn(x) * 100, digits = digits)
}
#Percentiles for 3 archetypes:
.perc(parameters(a3),smpulg) 
#Percentiles for 7 archetypes:
.perc(parameters(a7), smpulg) 

#Which is the nearest individual to archetypes?.
#Example for three archetypes:
i = 3
ai <- bestModel(lass15[[i]])
ras <- rbind(parameters(ai),smpulg)
dras <- dist(ras, method = "euclidean", diag =F, upper = T, p = 2)
mdras <- as.matrix(dras)
diag(mdras) = 100000000000
wh1 <- which.min(mdras[1,]) - (i) 
min(mdras[1,])
wh2 <- which.min(mdras[2,]) - (i)
min(mdras[2,])
wh3 <- which.min(mdras[3,]) - (i)
min(mdras[3,])
pa <- parameters(ai)
dist(rbind(pa[1,],smpulg[wh1,]))
smpulg[wh1,]

#In addition, we can turn the standarized values to the original variables.
parameters(a7)
p <- parameters(a7)
m <- mean(mpulg[,])
s <- sd(mpulg[,])
d <- p
for(i in 1 : 6){
 d[,i] = p[,i] * s[i] + m[i]
}
t(d)

#The code of function xyplot.pca is the following:
xyplot.pca <- function(x, y, data.col = 1, data.pch = 19, data.bg = NULL, atypes.col = 2, atypes.pch = 19, ahull.show = F, ahull.col = atypes.col,
chull = NULL, chull.col = gray(0.7), chull.pch = 19, adata.show = FALSE, adata.col = 3, adata.pch = 13, link.col = data.col, link.lty = 1, ...){
 zs <- x
 data <- y
 
 plot(data, col = data.col, pch = data.pch, bg = data.bg, ...)
 points(zs, col = atypes.col, pch = atypes.pch, ...)

 if(!is.null(chull)){
  points(data[chull, ], col = chull.col, pch = chull.pch, ...)
  lines(data[c(chull, chull[1]), ], col = chull.col, ...)
  }

 if(ahull.show)
  lines(ahull(zs), col = ahull.col)
 if(adata.show){
  adata <- fitted(zs)
  link.col <- rep(link.col, length = nrow(adata))
  link.lty <- rep(link.lty, length = nrow(adata))
  points(adata, col = adata.col, pch = adata.pch, ...)
  for (i in seq_len(nrow(data)))
   lines(rbind(data[i, ], adata[i, ]), col = link.col[i],
   lty = link.lty[i], ...)
  }
 invisible(NULL)
}




