#Code for examples in H-plots for displaying non-metric dissimilarity matrices
#Code for h-plots, cmdscale, isoMDS, Sammon, Laub et al.
#Isomap, multimaps, and other manifold learning
#techniques in matlab


#####################
#Congruence coefficient
####################


"congru" <- function(x,y) {
# x, y are the dissimilarity matrices
# congruence coefficient

vx=x[lower.tri(x)]


vy=y[lower.tri(y)]

cxy=sum(vx*vy)/(sqrt(sum(vx^2))*sqrt(sum(vy^2)))

return(cxy)
}


library(MASS)



#####################
#Example 1 
####################

#dissimilarity matrix

n=4
d=matrix(1.5,n,n)
diag(d)=0

d[1,2]=1
d[2,1]= 1

d[1,3]=5
d[3,1]=5
d[4,1]=7
d[1,4]=7


d[2,3]=10
d[3,2]=10
d[4,2]=12
d[2,4]=12


cit=c('MA',  'VL',  'MO',  'SP')

#hplot, with 2 dimensions
p=prcomp(d)
v=p$rotation
pr=p$sdev
pv=cbind(pr[1]*v[1:n,1],pr[2]*v[1:n,2])

plot(pv,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
text(pv[,1],pv[,2],cit,cex=2)
as.matrix(dist(pv))
congru(d,as.matrix(dist(pv)))

#h-plot in one dimension
as.matrix(dist(pv[,1]))
congru(d,as.matrix(dist(pv[,1])))

#cmdscale

ms= cmdscale(d,k=1) #with k=2 the second eigenvalue is negative

plot(ms[,1],c(0,0,0,0),type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
text(ms[,1],c(0,0,0,0),cit,cex=2)

as.matrix(dist(ms))
congru(d,as.matrix(dist(ms)))

#isomds, in 2-dimensions
y = cmdscale(d,k=2, add=T) 

ms=isoMDS(d,y$points)

plot(ms$points)
plot(ms$points,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
text(ms$points[,1],ms$points[,2],cit,cex=2)

#isomds in 1-dim
ms=isoMDS (d, y = cmdscale(d, 1),k=1)
plot(ms$points[,1],c(0,0,0,0),type='n',asp=1,xlab='',ylab='')
text(ms$points[,1],c(0,0,0,0),cit)
as.matrix(dist(ms$points)
congru(d,as.matrix(dist(ms$points)))

#Sammon in 2-d

ms=sammon(d,y$points) 


plot(ms$points,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
text(ms$points[,1],ms$points[,2],cit,cex=2)

#Sammon in 1-d

ms=sammon(d, y = cmdscale(d, 1),k=1)
plot(ms$points[,1],c(0,0,0,0),type='n',asp=1,xlab='',ylab='')
text(ms$points[,1],c(0,0,0,0),cit,cex=2)
as.matrix(dist(ms$points)
congru(d,as.matrix(dist(ms$points)))

#Isomap in matlab 

ii=read.table('aeroisomap.dat',sep=" ") 

congru(d,as.matrix(dist(ii)))

plot(ii,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
text(ii[,1],ii[,2],cit,cex=2)



#####################
#Example 4 
####################


d1=matrix(c(0.00, 2.36, 2.59, 1.78, 4.74, 4.82 ,4.98, 4.72,2.36, 0.00, 2.39, 1.60 ,4.98, 
5.06, 5.22, 4.96,2.59, 2.39, 0.00, 2.09, 5.29, 5.37, 5.53 ,5.27,1.78, 1.60 ,2.09, 0.00,
 5.08, 5.16, 5.32, 5.06,4.74, 4.98, 5.29, 5.08 ,0.00, 1.20, 1.82, 1.62,4.82, 5.06, 5.37, 5.16,
  1.20, 0.00, 2.98, 1.78,4.98, 5.22, 5.53 ,5.32, 1.82, 2.98, 0.00, 2.02,
4.72, 4.96, 5.27 ,5.06, 1.62, 1.78, 2.02, 0.00),8)


d2=matrix(c(0.00, 4.15, 2.03, 4.14, 1.26, 4.33, 0.690, 4.85,
4.15 ,0.00 ,4.70, 0.570, 4.37 ,1.82, 4.24, 2.02,
2.03, 4.70 ,0.00 ,4.69 ,1.85, 4.88, 1.68 ,5.40,
4.14, 0.570, 4.69 ,0.00 ,4.36, 1.83, 4.23, 2.67,
1.26, 4.37 ,1.85 ,4.36, 0.00, 4.55, 0.730, 5.07,
4.33 ,1.82 ,4.88 ,1.83 ,4.55 ,0.00 ,4.42, 2.14,
0.690 ,4.24 ,1.68, 4.23, 0.730, 4.42, 0.00, 4.94,
4.85, 2.02, 5.40, 2.67, 5.07, 2.14, 4.94 ,0.00),8)

d=d1-d2

n=8
e=rep(1,n)

Q=(diag(n) - e%*%t(e)/n)

C=-0.5* (diag(n) - e%*%t(e)/n)%*%d%*%(diag(n) - e%*%t(e)/n)

E=eigen(C)

X=E$vectors%*%diag(sqrt(abs(E$values)))


pl=cbind(sqrt(E$values[1])*E$vectors[1:n,1],sqrt(abs(E$values[2]))*E$vectors[1:n,2])
plot(pl,type='n',asp=1,xlab='',ylab='',cex=1.5)
text(pl[,1],pl[,2],1:n,cex=2)

windows()
pll=cbind(sqrt(abs(E$values[n-1]))*E$vectors[1:n,n-1],sqrt(abs(E$values[n]))*E$vectors[1:n,n])
plot(pll,type='n',asp=1,xlab='',ylab='',cex=1.5)
text(pll[,1],pll[,2],1:n,cex=2)

#h-plot

n=8
p=princomp(d)
v=p$loadings
pr=p$sdev
pv=cbind(pr[1]*v[1:n,1],pr[2]*v[1:n,2])


#plot(pv,type='n',asp=1,xlab='Dim 1',ylab='Dim 2')
plot(pv,type='n',asp=1,xlab='',ylab='',cex=1.5)
text(pv[,1],pv[,2],1:n,cex=2)


#multimaps in matlab


iiw=read.table('laubmultiw.dat',sep=" ") 

#map 2
ii=read.table('laubmulti2.dat',sep=" ") 
plot(ii,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
points(ii[,1],ii[,2],pch=1,cex=10*iiw[,2])
text(ii[,1],ii[,2],1:n,cex=2)

#####example 1, Laub et al. ####

#d is d in example 1

n=4
e=rep(1,n)

Q=(diag(n) - e%*%t(e)/n)

C=-0.5* (diag(n) - e%*%t(e)/n)%*%d%*%(diag(n) - e%*%t(e)/n)

E=eigen(C)

X=E$vectors%*%diag(sqrt(abs(E$values)))


pl=cbind(sqrt(E$values[1])*E$vectors[1:n,1],sqrt(abs(E$values[2]))*E$vectors[1:n,2])
  
plot(pl,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
text(pl[,1],pl[,2],cit,cex=2)

windows()
pll=cbind(sqrt(abs(E$values[n-1]))*E$vectors[1:n,n-1],sqrt(abs(E$values[n]))*E$vectors[1:n,n])
plot(pll,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
text(pll[,1],pll[,2],cit,cex=2)

#multimaps in matlab

ii=read.table('aeromulti1.dat',sep=" ") 
iiw=read.table('aeromultiw.dat',sep=" ") 

#map 1
plot(ii,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
points(ii[,1],ii[,2],pch=1,cex=10*iiw[,1])
text(ii[,1],ii[,2],cit,cex=2)

#map 2
ii=read.table('aeromulti2.dat',sep=" ") 
plot(ii,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
points(ii[,1],ii[,2],pch=1,cex=10*iiw[,2])
text(ii[,1],ii[,2],cit,cex=2)

#####################
#Point patterns 
####################


#####################
#Experiment 1: suppl. material
####################
load("exp1disiG3lg.bin")

d<-disiG3lg #dissimilarity matrix computed in Ayala et al. CSDA2006

#h-plot
p=princomp(d)
v=p$loadings
pr=p$sdev
pv=cbind(pr[1]*v[1:300,1],pr[2]*v[1:300,2])



plot(pv,type='n',xlab="Dim 1", ylab="Dim 2")
points(pv[1:100,],col=4,pch=22) #blue
points(pv[101:200,],col=2,pch=21) #red
points(pv[201:300,],col=3,pch=24) #green

congru(d,as.matrix(dist(pv)))

#####################
#Corneal endothelia 
####################


load("dissimilaritiesendothelia.bin")


d<-disiETlg #dissimilarity matrix computed in Ayala et al. CSDA2006

library(cluster)
h<-as.dist(d)


k<-pam(h,k=2,diss=T)

p=princomp(d)
v=p$loadings
pr=p$sdev
pv=cbind(pr[1]*v[1:153,1],pr[2]*v[1:153,2])

congru(d,as.matrix(dist(pv)))


plot(pv[,1],pv[,2],type='n',xlab="Dim 1", ylab="Dim 2") 
points(pv[k$clustering==1,1],pv[k$clustering==1,2],pch=21,col=1) 
points(pv[k$clustering==2,1],pv[k$clustering==2,2],pch=24,col=2) 

#ANOVA


library(vegan)

#cmdscale
ms= cmdscale(d)
 

anod <- anosim(dist(ms), as.factor(k$clustering))

ado<- adonis(dist(ms)~ as.factor(k$clustering))

mr <- mrpp(dist(ms), as.factor(k$clustering))

#isomds
ms=isoMDS(d)


anod <- anosim(dist(ms$points), as.factor(k$clustering))

ado<- adonis(dist(ms$points)~ as.factor(k$clustering))

mr <- mrpp(dist(ms$points), as.factor(k$clustering))


#sammon
ms=sammon(d) 

anod <- anosim(dist(ms$points), as.factor(k$clustering))

ado<- adonis(dist(ms$points)~ as.factor(k$clustering))

mr <- mrpp(dist(ms$points), as.factor(k$clustering))

#h-plot

p=princomp(d)
v=p$loadings
pr=p$sdev
pv=cbind(pr[1]*v[1:153,1],pr[2]*v[1:153,2])



anod <- anosim(dist(pv), as.factor(k$clustering))

ado<- adonis(dist(pv)~ as.factor(k$clustering))

mr <- mrpp(dist(pv), as.factor(k$clustering))

#isomap

 ii=read.table('endisomape1.dat',sep=" ") # isomapR1 with epsilon 1: 146 points
 
 #all except that 7 seven points eliminanted in isomap with epsilon=1
 o=c(25,27,29,44,84,136,153)
 pu=setdiff(1:153,o)
 
ms=ii;


anod <- anosim(dist(ms), as.factor(k$clustering[pu]))

ado<- adonis(dist(ms)~ as.factor(k$clustering[pu]))

mr <- mrpp(dist(ms), as.factor(k$clustering[pu]))



#rank dissimilarities and h-plot

ds= matrix(rank(d),153,153)

p=princomp(ds)
v=p$loadings
pr=p$sdev

pv=cbind(pr[1]*v[1:153,1],pr[2]*v[1:153,2])


anod <- anosim(dist(pv), as.factor(k$clustering))

ado<- adonis(dist(pv)~ as.factor(k$clustering))

mr <- mrpp(dist(pv), as.factor(k$clustering))

#####################
#Example 3 with a = 2: suppl. material 
####################


la=c('A', 'B', 'C', 'D', 'E')


n=5
a=2
d=matrix(1,5,5)

d[1,1]=0
d[2,2]=0
d[3,3]=0
d[4,4]=0
d[5,5]=0
d[1,2]=a
d[2,1]=a
d[5,4]=4 
d[4,5]=d[5,4]
d[5,3]=2
d[3,5]=d[5,3]

#h-plot

p=princomp(d)
v=p$loadings
pr=p$sdev
pv=cbind(pr[1]*v[1:n,1],pr[2]*v[1:n,2])


plot(pv,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
text(pv[,1],pv[,2],la,cex=3)              

congru(d,as.matrix(dist(pv)))

#cmdscale
ms= cmdscale(d)


plot(ms,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
text(ms[,1],ms[,2],la,cex=3)

congru(d,as.matrix(dist(ms)))


#isomds
ms=isoMDS(d)

congru(d,as.matrix(dist(ms$points)))


plot(ms$points,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
text(ms$points[,1],ms$points[,2],la,cex=3)

#sammon
ms=sammon(d) 

congru(d,as.matrix(dist(ms$points)))

plot(ms$points,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
text(ms$points[,1],ms$points[,2],la,cex=3)

#isomap in matlab

#computing congruence coefficient for isomap
ii=read.table('brand2isomap.dat',sep=" ") 

congru(d,as.matrix(dist(ii)))

plot(ii,type='n',asp=1,xlab='',ylab='',cex.axis=1.5)
text(ii[,1],ii[,2],la,cex=3)



#####################
#Example 4: suppl. material 
####################

m=read.table('morse code confusions asym.dat') #from http://people.few.eur.nl/groenen/mmds/datasets/


 
#first those that are presented in first place, i.e. di. from columns to 1 to 36, while d.i from column 37 to 72 
da=cbind(t(m),m) #dimension: (36,36*2)

# di., where code i-th is first presented
#d.i, where code i-th is second presented

  
#changing from similarity to dissimilarity
dda=1-da


la=c(LETTERS ,1:10)
la[36]="0"

p=prcomp(dda)
v=p$rotation
pr=p$sdev
pv=cbind(pr[1]*v[1:72,1],pr[2]*v[1:72,2])

plot(pv,type='n',asp=1,xlab='Dim 1',ylab='Dim 2')
text(pv[1:36,1],pv[1:36,2],la)
text(pv[37:72,1],pv[37:72,2],la,col=2)



