d  <- read.table("inp_iris.txt",skip=1)
#d  <- read.table("inp_G54.csv",skip=1,sep=",")
dd <- d[,1:4] 
nn <- nrow(dd)
mm <- ncol(dd)
kk <- 4

#Standardization
#dd <- sweep(dd,2,apply(dd,2,mean),FUN="-")
#dd <- sweep(dd,2,apply(dd,2,sd),FUN="/")

cl <- as.integer(runif(nn,min=1,max=kk+1))
xx <- as.matrix(dd)
a <-matrix(0,nrow=kk,ncol=mm)
for(i in 1:nn){
    for(k in 1:kk){
        a[k,1:mm] <- apply(subset(xx,cl==k),2,mean)
    }
}
disx <-matrix(0,nrow=nn,ncol=kk)
for(i in 1:nn){
    for(k in 1:kk){
        disx[i,k]<-sum((xx[i,1:mm]-a[k,1:mm])^2)
    }
}
#*************************************
for(j in 1:20){
    for(i in 1:nn){
        for(k in 1:kk){
            if(min(disx[i,seq(1,kk,by=1)])==disx[i,k])cl[i] <- k
        }
    }
    for(i in 1:nn){
        for(k in 1:kk){
            a[k,1:mm] <- apply(subset(xx,cl==k),2,mean)
        }
    }
    for(i in 1:nn){
        for(k in 1:kk){
            disx[i,k]<-sum((xx[i,1:mm]-a[k,1:mm])^2)
        }
    }
}
#*************************************
cl
#***************************************************
ii  <- matrix(1,nrow=nn,ncol=nn)
mat <- diag(nn)-1/nn*ii
yy  <- xx %*% t(xx)
dis <- (diag(nn)*yy) %*% ii - 2*yy + ii %*% (diag(nn)*yy) #Distance matrix
ary <- -0.5*mat %*% dis %*% mat #Young-Householder's transformation
yy  <- eigen(ary)
emat<- matrix(0,nrow=nn,ncol=nn)
for(i in 1:nn){
    if(0<=yy$values[i])emat[i,i]<-sqrt(yy$values[i])
    if(yy$values[i]<0)emat[i,i]<-0
}
plt <- yy$vectors %*% emat
panel.text<-function(x, y, labels, cex, font, ...){
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0,1,0,1) )
  txt <- c("MDS")
  cex<-2
  font<-"Helvetica"
  text(0.5, 0.5, txt, cex = cex)
} 
pairs(
  plt[,1:4],
  diag.panel=NULL,
  text.panel=panel.text,
  pch=21,
  bg=c("red","yellow","blue","cyan","green","black")[unclass(cl)]
  )