


library(MASS)
library(archetypes)

huge=200

no.scalefn <- function(x, ...){
  return(x)
}
no.rescalefn <- function(x, zs, ...){
  if ( is.null(zs) )
    return(matrix(NA, nrow = 0, ncol = 0))
  return(zs)
}

ginv.zalphasfn <- function(alphas, x, ...) {
  require(MASS)  
  return(t(ginv(alphas %*% t(alphas)) %*% alphas %*% t(x)))
}

stepArchetypes_mod <- function (data, k, nrep = 3, verbose = TRUE){
  mycall <- match.call()
  as <- list()
  for (i in 1:length(k)) {
    as[[i]] <- list()
    class(as[[i]]) <- "repArchetypes"
    for (j in seq_len(nrep)) {
      if (verbose) 
        cat("\n*** k=", k[i], ", rep=", j, ":\n", sep = "")
      as[[i]][[j]] <- archetypes(data, k = k[i], family = archetypesFamily("original",scalefn = no.scalefn,rescalefn = no.rescalefn))
    }
  }
  return(structure(as, class = "stepArchetypes", call = mycall))
}


set.seed(2013)

m=read.table('mfdgrsqcurvas50_orto2.txt')
X=m[,2:50]
nm=dim(X)

lass10 <- stepArchetypes_mod(data = X, k = 1:10, verbose = FALSE, nrep = 5) 
screeplot(lass10) #Figure 3


i=3
na=i
ai=bestModel(lass10[[i]])

#alphas
gv <- rbind(t(parameters(ai)), rep(200, na))
zs <- as.matrix(gv)

x_gvv <- rbind(t(X), rep(200, na))
zs <- as.matrix(gv)

alphas <- matrix(0, nrow = na, ncol = nm[1])
for (j in 1 : nm[1]){
 alphas[, j] = coef(nnls(zs, x_gvv[,j]))
}
#alphas[,ar]


#write.table(t(alphas),'alfas3se_orto2.txt',col.names=F,row.names=F) 

#alphas[,eti==1]
#alphas[,eti==2]
#alphas[,eti==3]
#alphas[,eti==0]

#aa=read.table('alfas3se_orto2.txt')

aa=t(alphas)


library(clValid)
library(pdfCluster)


library(trimcluster)



set.seed(1234)

km=kmeans(aa, 3, iter.max = 10, nstart = 20)

cl=km$cluster

##################

m=read.table('mfdall_orto2.txt')
X=m[,2:50]
nm=dim(X)



i=3
na=i
ai=bestModel(lass10[[i]])

#alphas
gv <- rbind(t(parameters(ai)), rep(200, na))
zs <- as.matrix(gv)

x_gvv <- rbind(t(X), rep(200, na))
zs <- as.matrix(gv)

alphas <- matrix(0, nrow = na, ncol = nm[1])
for (j in 1 : nm[1]){
 alphas[, j] = coef(nnls(zs, x_gvv[,j]))
}



aa=t(alphas)

dd=c()
eall=c()
for( j in 1:dim(aa)[1]){
 dd=as.matrix(dist(rbind(aa[j,],km$centers)))
eall[j]=which.min(dd[2:4,1])
}


write.table(eall,'eallorto2.txt',col.names=F,row.names=F)

#########################

##############
#k-means
##############

m=read.table('mfdgrsqcurvas50_orto2.txt')
X=m[,2:50]



set.seed(1234)

km=kmeans(X, 3, iter.max = 10, nstart = 20)



##################

m=read.table('mfdall_orto2.txt')
X=m[,2:50]
nm=dim(X)




dd=c()
eall=c()
for( j in 1:dim(X)[1]){
 dd=as.matrix(dist(rbind(X[j,],km$centers)))
eall[j]=which.min(dd[2:4,1])
}


write.table(eall,'eallorto2k.txt',col.names=F,row.names=F)