
#Figure 1


#waveform

casose <-c('we1.txt')

i=1

we=read.table(casose[i])
X=we[,1:101]

t=seq(1,21,0.2)
y=t(X)

#########
#Sample: 5 of each group


par(mfrow=c(3,1),mar = c(2, 3, 1, 1)) 
 matplot(t, y[,1:5], type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,ylim=c(-2.5,6.5)) #data
matplot(t, y[,151:160], type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,ylim=c(-2.5,6.5)) #data
matplot(t, y[,301:305], type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,ylim=c(-2.5,6.5)) #data


#original functions source generating data


 matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data

h1=pmax(6 ? abs(t ?11), 0)
h2=pmax(6 ? abs(t-4 ?11), 0)
h3=pmax(6 ? abs(t+4 ?11), 0)

#matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data
lines(t, h1, type = "l",lty=1,xlab="",ylab="",col=2,lwd=3) #h1
lines(t, h2, type = "l",lty=2,xlab="",ylab="",col=3,lwd=3) #h2
lines(t, h3, type = "l",lty=3,xlab="",ylab="",col=4,lwd=3) #h3



################# AA
library(archetypes)
norep=20
K=3
 set.seed(1234)
# lassd <- stepLArchetypoids3(data=X, k=3, norep) # no robust

#out=lassd[[1]][[1]]


 set.seed(1234)
lass <- LstepArchetypesMod(data = X, k=K, nrep = norep,verbose = FALSE)
ai <- bestModel(lass[[1]]) #archetypes



a3a=ai$archetypes #archetypes

#a3a=t(y[,lassd[[1]][[1]]])



 matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data



#matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data
lines(t, a3a[1,], type = "l",lty=1,xlab="",ylab="",col=3,lwd=3) #arch
lines(t, a3a[2,], type = "l",lty=2,xlab="",ylab="",col=4,lwd=3) #arch
lines(t, a3a[3,], type = "l",lty=3,xlab="",ylab="",col=2,lwd=3) #arch

####### PCA

p=princomp(X)


 matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data

#PC1


#matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data
lines(t, apply(X,2,mean)+10*p$loadings[,1], type = "l",lty=1,xlab="",ylab="",col=2,lwd=3) #pc+ c1
lines(t, apply(X,2,mean)-10*p$loadings[,1], type = "l",lty=2,xlab="",ylab="",col=6,lwd=3) #pc-c1


#PC2

 matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data


lines(t, apply(X,2,mean)+10*p$loadings[,2], type = "l",lty=3,xlab="",ylab="",col=4,lwd=3) #pc+ c1
lines(t, apply(X,2,mean)-10*p$loadings[,2], type = "l",lty=4,xlab="",ylab="",col=5,lwd=3) #pc-c1



#cluster

set.seed(1234)
km=kmeans(X, 3, iter.max = 10, nstart = norep)

 matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data



#matplot(t, y, type = "l",lty=1,xlab="",ylab="",cex.lab=1.5,cex.axis=1.5,col=8) #data
lines(t, km$centers[1,], type = "l",lty=1,xlab="",ylab="",col=2,lwd=3) #center
lines(t, km$centers[2,], type = "l",lty=2,xlab="",ylab="",col=3,lwd=3) #center
lines(t, km$centers[3,], type = "l",lty=3,xlab="",ylab="",col=4,lwd=3) #center


