#' ---
#' title: Experiment with Zebrafish brain data
#' output:
#'    html_document:
#'      toc: true
#'      css: ~/RFolder/pandoc.css
#'      highlight: pygments
#' ---

load("tsg-11-30-469.Rbin")

set.seed(12345)

# create the dissimilarity matrices
D<-list()
for(i in 1:20){
  A<-as.matrix(tsg21[[i]])
  g<-graph.adjacency(A,mode=c("undirected"))
  D[[i]]<-1-similarity(g,method=c("jaccard"))
}

w=10

t<-proc.time();smds1 <- fJOFC(10,D,2,verbose=FALSE);runtime1<-proc.time()-t

# make figure 8
w<-c(1,10,25,50,75,100)
DD<-list()
for(i in 1:6)
{
  smds <- fJOFC(w[i],D,2,verbose=FALSE)
  X<-list()
  for(j in 1:20){
    X[[j]]<-smds$conf[(469*(j-1)+1):(469*j),]
  }
  DD[[i]]<-matrix(0,20,20)
  for(j in 1:20){for(k in 1:20){DD[[i]][j,k]<-norm(X[[j]]-X[[k]],"f")}}
}
par(mfrow = c(2,3),
    oma = c(5,4,0,0) + 0.3,
    mar = c(1,0,2,1) + 0.3)
x<-DD[[1]]
rownames(x)<-c(30:11)
colnames(x)<-c(11:30)

image(x[,20:1],xaxt='n',yaxt='n',main="w=1")
axis( 1, at=seq(0,1,length.out=ncol( x ) ),cex.axis=0.7, labels= colnames( x ), las= 2 )
axis( 2, at=seq(0,1,length.out=nrow( x ) ),cex.axis=0.7, labels= rownames( x ), las= 2)
image(DD[[2]][,20:1],xaxt='n',yaxt='n',main="w=10")
image(DD[[3]][,20:1],xaxt='n',yaxt='n',main="w=25")
image(DD[[4]][,20:1],xaxt='n',yaxt='n',main="w=50")
image(DD[[5]][,20:1],xaxt='n',yaxt='n',main="w=75")
image(DD[[6]][,20:1],xaxt='n',yaxt='n',main="w=100")

# make Figure 9

X<-list()
for(i in 1:20){
  X[[i]]<-smds1$conf[(469*(i-1)+1):(469*i),]
}
par(mfrow = c(4,5),
    oma = c(5,4,0,0) + 0.1,
    mar = c(0,0,1,1) + 0.1)
plot(X[[1]],ylab="",xlab="",main="tau=11",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[2]],ylab="",xlab="",main="tau=12",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[3]],ylab="",xlab="",main="tau=13",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[4]],ylab="",xlab="",main="tau=14",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[5]],ylab="",xlab="",main="tau=15",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[6]],ylab="",xlab="",main="tau=16",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[7]],ylab="",xlab="",main="tau=17",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[8]],ylab="",xlab="",main="tau=18",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[9]],ylab="",xlab="",main="tau=19",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[10]],ylab="",xlab="",main="tau=20",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[11]],ylab="",xlab="",main="tau=21",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[12]],ylab="",xlab="",main="tau=22",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[13]],ylab="",xlab="",main="tau=23",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[14]],ylab="",xlab="",main="tau=24",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[15]],ylab="",xlab="",main="tau=25",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[16]],ylab="",xlab="",main="tau=26",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[17]],ylab="",xlab="",main="tau=27",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[18]],ylab="",xlab="",main="tau=28",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[19]],ylab="",xlab="",main="tau=29",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
plot(X[[20]],ylab="",xlab="",main="tau=30",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')

# Figure 10

par(mfrow = c(2,2),
    oma = c(5,4,0,0) + 0.1,
    mar = c(0,0,1,1) + 0.1)

plot(X[[11]],ylab="",xlab="",main="tau=21--22",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
points(X[[12]])
for(i in 1:469){segments( X[[11]][i,1],X[[11]][i,2],X[[12]][i,1],X[[12]][i,2],col='blue')}

plot(X[[12]],ylab="",xlab="",main="tau=22--23",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
points(X[[13]],col='red')
for(i in 1:469){segments( X[[12]][i,1],X[[12]][i,2],X[[13]][i,1],X[[13]][i,2],col='blue')}

plot(X[[13]],ylab="",xlab="",main="tau=23--24",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n',col='red')
points(X[[14]])
for(i in 1:469){segments( X[[13]][i,1],X[[13]][i,2],X[[14]][i,1],X[[14]][i,2],col='blue')}

plot(X[[14]],ylab="",xlab="",main="tau=24--25",xlim=c(-5,5),ylim=c(-5,5),xaxt='n',yaxt='n')
points(X[[15]])
for(i in 1:469){segments( X[[14]][i,1],X[[14]][i,2],X[[15]][i,1],X[[15]][i,2],col='blue')}


## run JOFC for time conparison
# make the weight matrix
# n=469
# w=10
# W1 <- matrix(1,n,n) - diag(n)
# W2 <- diag(w,n)
# W <- W2
# for(i in 1:19){W<-cbind(W,W2)}
# WW<-W
# for(i in 1:19){W<-rbind(W,WW)}
# W<-W-diag(w,469*20)
# W <- W+as.matrix(bdiag(rep(list(W1),20)))
# 
# # make the omnibus matrix
# d<-list()
# for(i in 1:20){
#   d[[i]]<-D[[i]]
#   for(l in 1:19){d[[i]]<-cbind(d[[i]],D[[i]])}
# }
# O<-Reduce(rbind,d)
# O<-O/2+t(O)/2
# 
# t<-proc.time();smds2 <- smacofSym(O,2,weightmat=W,verbose=FALSE);runtime2<-proc.time()-t