#[0]--------------------------
if(!require(cluster)){install.packages("cluster"); library(cluster)}
# input data----------------------------------
df <- read.csv("DataBrandy.csv",row.names=1,header=TRUE)
nf<-dim(df)[[2]]; df<-scale(df)
#[1] Gap statistics---------
rmax<-15  # Mannually determined
opKmean <- clusGap(df, FUNcluster=kmeans,K.max = rmax, B=1000, verbose = interactive())
#opKmean <- clusGap(df, FUNcluster=clara,K.max = rmax, B=1000, verbose = interactive())
plot(opKmean,ylim=c(0,0.8))
ns<-dim(opKmean$Tab)[1]
#[2] diff ------------------
sabun<-c(rep(0,rmax-1))
opK<-opKmean$Tab[,3]
for (i in 1:(rmax-1)) {              
  sabun[i] <- opK[i+1] -opK[i]          
}
barplot(sabun,ylim=c(-0.1,0.2),names.arg=c(2:ns),
        xlab="k",ylab="ąGap", main="k vs. ąGap")
abline(h=0)
#-[3] visualization of clusters------------
maxc<-2  # Manually determined
df.km <- kmeans(df,maxc)
clusplot(df,df.km$cluster, color=TRUE, shade=TRUE, labels=2, lines=0)
dflabel <- as.factor(substring(rownames(df),1,2))
table <- table(dflabel,df.km$cluster); table
#[4]  Output-------------------------------------------------------------
fdata<-cbind(df.km$cluster,df)
colnames(fdata)[[1]]<-c("GpEst"); filename<-c("DataKmeans.csv")
write.csv(fdata,file=filename)