#' ---
#' title: Error Tolerance
#' output:
#'    html_document:
#'      toc: true
#'      css: ~/RFolder/pandoc.css
#'      highlight: pygments
#' ---

# code to generate Figure 2
source("packagesatstart_JCGS.R")

set.seed(12345)

# rel1 = relative error for fJOFC in the clean setting
# rel2 = relative error for fJOFC in the anomaly setting

rel1<-list()
rel2<-list()

# e = indices of anomalous vertices
e<-list()
e[[1]]<-c(1:10, 401:411, 801:810, 1201:1210)
e[[2]]<-c(1:10, 401:411, 801:810, 1201:1210,1601:1610)
e[[3]]<-c(1:10, 401:411, 801:810, 1201:1210,1601:1610,2001:2010)

for(i in 1:3){rel1[[i]]<-matrix(0,4,25)
rel2[[i]]<-matrix(0,4,25)}

# m = number of modalities
m <- c(4, 5, 6)

for(k in 1:3){
for(j in 1:25){
  X0 <- matrix(rnorm(400*2,5,1),400,2)
  X<-list()
  for(l in 1:m[k]){ X[[l]]=jitter(X0,factor=1,amount=0)}
  Y0<-matrix(rnorm(20,8,2),10,2)
  Y0<-rbind(Y0,X0[11:400,])
  Y<-list()
  for(l in 1:(m[k]-1)){ Y[[l]]=jitter(X0,factor=1,amount=0)}
  Y[[ m[k] ]]=jitter(Y0,factor=1,amount=0)
  
  DX <- list()
  for(l in 1:m[k]){ DX[[l]]=as.matrix(dist(X[[l]]))}
  DY <- list()
  for(l in 1:m[k]){ DY[[l]]=as.matrix(dist(Y[[l]]))}
  
  smds1 <- fJOFC(200,DX,2,verbose=FALSE)
  xfinal<-smds1$conf
  
  smds1i2 <- fJOFC(200,DX,2,verbose=FALSE,itmax=2)
  xit2<-smds1i2$conf
  rel1[[k]][1,j]<-norm(xit2-xfinal,"f")/norm(xfinal,"f")
  
  smds1i3 <- fJOFC(200,DX,2,verbose=FALSE,itmax=3)
  xit3<-smds1i3$conf
  rel1[[k]][2,j]<-norm(xit3-xfinal,"f")/norm(xfinal,"f")
  
  smds1i4 <- fJOFC(200,DX,2,verbose=FALSE,itmax=4)
  xit4<-smds1i4$conf
  rel1[[k]][3,j]<-norm(xit4-xfinal,"f")/norm(xfinal,"f")
  
  smds1i5 <- fJOFC(200,DX,2,verbose=FALSE,itmax=5)
  xit5<-smds1i5$conf
  rel1[[k]][4,j]<-norm(xit5-xfinal,"f")/norm(xfinal,"f")
  
  smds2 <- fJOFC(200,DY,2,verbose=FALSE)
  yfinal<-smds2$conf
  
  smds2i2 <- fJOFC(200,DY,2,verbose=FALSE,itmax=2)
  yit2<-smds2i2$conf
  rel2[[k]][1,j]<-norm(yit2[-e[[k]],]-yfinal[-e[[k]],],"f")/norm(yfinal[-e[[k]],],"f")
  
  smds2i5 <- fJOFC(200,DY,2,verbose=FALSE,itmax=5)
  yit5<-smds2i5$conf
  rel2[[k]][2,j]<-norm(yit5[-e[[k]],]-yfinal[-e[[k]],],"f")/norm(yfinal[-e[[k]],],"f")
  
  smds2i10 <- fJOFC(200,DY,2,verbose=FALSE,itmax=10)
  yit10<-smds2i10$conf
  rel2[[k]][3,j]<-norm(yit10[-e[[k]],]-yfinal[-e[[k]],],"f")/norm(yfinal[-e[[k]],],"f")
  
  smds2i25 <- fJOFC(200,DY,2,verbose=FALSE,itmax=25)
  yit25<-smds2i25$conf
  rel2[[k]][4,j]<-norm(yit25[-e[[k]],]-yfinal[-e[[k]],],"f")/norm(yfinal[-e[[k]],],"f")
  

}
}


means1<-matrix(0,12,1)
for(i in 1:3){for(j in 1:4){means1[(i-1)*4+j]<-mean(rel1[[i]][j,])}}
std1<-matrix(0,12,1)
for(i in 1:3){for(j in 1:4){std1[(i-1)*4+j]<-sd(rel1[[i]][j,])/5}}
p1<-cbind(means1,std1)
p1<-data.frame(p1)
m<-c(rep('4',4),rep('5',4),rep('6',4))
reps<-c(rep(c(2,3,4,5),3))
p1<-cbind(p1,m)
p1<-cbind(p1,reps)

g1<-ggplot(p1,aes(colour=m,x=reps, y=X1))+geom_errorbar(aes(ymin=X1-2*X2, ymax=X1+2*X2), width=0.2,size=1) +
  geom_line(size=1) +
  geom_point()+xlab("Maximum Iterations") +
  ylab("Relative Error")+
  ggtitle("Matched Setting")+
  theme(text = element_text(size=10))
g1


means2<-matrix(0,12,1)
for(i in 1:3){for(j in 1:4){means2[(i-1)*4+j]<-mean(rel2[[i]][j,])}}
std2<-matrix(0,12,1)
for(i in 1:3){for(j in 1:4){std2[(i-1)*4+j]<-sd(rel2[[i]][j,])/5}}
p2<-cbind(means2,std2)
p2<-data.frame(p2)
m<-c(rep('4',4),rep('5',4),rep('6',4))
reps<-c(rep(c(2,5,10,25),3))
p2<-cbind(p2,m)
p2<-cbind(p2,reps)

g2<-ggplot(p2,aes(colour=m,x=reps, y=X1))+geom_errorbar(aes(ymin=X1-2*X2, ymax=X1+2*X2), width=0.2,size=1) +
  geom_line(size=1) +
  geom_point()+xlab("Maximum Iterations") +
  ylab("Relative Error")+
  ggtitle("Anomaly Setting")+
  theme(text = element_text(size=10))
g2