


source("FstepArchetypesMod.R")
source("FstepArchetypoids.R")
source("FFunctions_to_calculate_real_archetypes_with_swap2.R")

library(archetypes)
library(fda)
library(R.matlab)


 Wfd=readMat('Wfd.mat')

Wfd=Wfd$Wfd

X=t(Wfd$coef) #coefficients

nbasis = 11;
norder = 6;
wbasis = create.bspline.basis(c(-2.8,2.8), nbasis);

PM=eval.penalty(wbasis, rng=c(-2.8,2.8))

WfdR=fd(coef=Wfd$coef, basisobj=wbasis)


K=4
set.seed(2018)

norep=20
lass10 <- FstepArchetypesMod(data=X,k=K,verbose=FALSE,nrep=norep,PM=PM,saveHistory=F) 
ai <- bestModel(lass10[[1]]) #archetypes
 
an4 <-  FstepArchetypoids(K,TRUE,X,lass10,FALSE,K-1,PM)

aw4 <-  FstepArchetypoids(K,FALSE,X,lass10,FALSE,K-1,PM) #which_alpha archetypoids

ini_arch=c()
 for (j in 1:K){
 ini_arch[j]=which.max(ai$betas[j,])
 }
 
abeta <- FrealArchetypes(X,huge = 200, K, ini_arch,PM) #which beta archetypoid

# 2 60 28 18



#Represent in matlab


#explained

1 -abeta$rss/(norm2.Fnormfn(t(X),PM)/dim(X)[1])
# 0.968896


#################
#alphas archetypoids

n=60
huge=200
 archs_ini <- c(2,18,28,60)
 data=X
  
  n <- ncol(t(data))
  x_gvv <- rbind(t(data), rep(huge, n))
  
  zs <- x_gvv[, archs_ini] 
  zs <- as.matrix(zs)
  
  alphas <- matrix(0, nrow = K, ncol = n)
  for (j in 1 : n){
    alphas[, j] = coef(nnls(zs, x_gvv[,j]))
  }


al=alphas
 

barplot(al,names.arg=1:60,cex.names=.6)

#barplot(al,col=c(2,8,4,7),names.arg=1:60,cex.names=.6)

stars(t(al),draw.segments=T,labels=1:60)


