nmds <- function(delta,x,Ntop)
{
# delta: N times N dissimilarity matrix
# x: inital configuration
# Ntop: Number of total iterations
N<-dim(delta)[1]
ds <- 0.1/N^3
rdelta <- matrix(rank(delta),ncol=N)
#----------- 1 of internal functions
rnorm <- function(x1)
{
avg <-  apply(x1,2,mean)
x1 <- t(t(x1)- avg)
x1/(sum(x1^2))^0.5
}
#----------- 2 of internal functions
drank <- function(x1)
{matrix(rank(as.matrix(dist(x1))),ncol=N)-rdelta}
#-------------------------
for (i in c(1:Ntop))
{
x <- rnorm(x)
cat(i,":",sum(drank(x)^2),"\n")
dr <- array(0,c(N,N,dim(x)[2]))
for (i in c(1:dim(x)[2]))
{dr[,,i] <- outer(x[,i],x[,i],FUN="-")}
r <- apply(dr^2,1,sum)^0.5
dr <- dr/r
x1 <- matrix(0,ncol=dim(x)[2],nrow=N)
for (i in c(1:dim(x)[2]))
{
diag(dr[,,i])<-0
x1[,i] <- x[,i] -ds*apply(dr[,,i]*drank(x),1,sum)
}
x<- x1
}
return(rnorm(x))
}