#' ---
#' title: JOFC vs fJOFC
#' output:
#'    html_document:
#'      toc: true
#'      css: ~/RFolder/pandoc.css
#'      highlight: pygments
#' ---

# make figure 3a

set.seed(12345)

rt1 <-matrix(0,5,50)
rt2 <-matrix(0,5,50)

n <- 400
m <- c(2,3,4,5,6)
w <- 200
W1 <- matrix(1,n,n) - diag(n)
W2 <- diag(w,n)

for(k in 1:5){
for(j in 1:50){
  # make the weight matrix
  W <- W2
  for(i in 1:(m[k]-1)){W<-cbind(W,W2)}
  WW<-W
  for(i in 1:(m[k]-1)){W<-rbind(W,WW)}
  W<-W-diag(w,n*m[k])
  W <- W+as.matrix(bdiag(rep(list(W1),m[k])))
  
  
  X0 <- matrix(rnorm(n*2,5,1),n,2)
  X<-list()
  for(l in 1:m[k]){ X[[l]]=jitter(X0,factor=1,amount=0)}
  DX <- list()
  for(l in 1:m[k]){ DX[[l]]=as.matrix(dist(X[[l]]))
  DX[[l]]<-normDiss(DX[[l]],1)}
  
  # make the omnibus matrix
  d<-list()
  for(i in 1:m[k]){
    d[[i]]<-DX[[i]]
    for(l in 1:(m[k]-1)){d[[i]]<-cbind(d[[i]],DX[[i]])}
  }
  O<-Reduce(rbind,d)
  O<-O/2+t(O)/2
  
  # initialize at cMDS
  X0<-torgerson(sqrt(O),2)
  
  t<-proc.time();smds1 <- fJOFC(w,DX,2,verbose=FALSE,init=X0);runtime1<-proc.time()-t
  rt1[k,j]<-runtime1[[3]]/(smds1$niter)
  t<-proc.time();smds2 <- smacofSym(O,2,weightmat=W,init=X0,verbose=FALSE);runtime2<-proc.time()-t
  rt2[k,j]<-runtime2[[3]]/(smds2$niter)
  
}}

means<-rep(0,10)
for(i in 1:5){means[i]=mean(rt1[i,])}
for(i in 6:10){means[i]=mean(rt2[i-5,])}
std<-rep(0,10)
for(i in 1:5){std[i]=2*sd(rt1[i,])/sqrt(50)}
for(i in 6:10){std[i]=2*sd(rt2[i-5,])/sqrt(50)}
x<-data.frame(cbind(means,std))
M<-c(m,m)
x<-cbind(x,M)
Method<-c(rep("fJOFC",5),rep("JOFC",5))
x<-cbind(x,Method)
g1<-ggplot(x,aes(colour=Method,x=M, y=means))+geom_errorbar(aes(ymin=means-std, ymax=means+std), width=0.2,size=1) +
       geom_line(size=1) +
       geom_point()+xlab("m") +
      ylab("Average runtime per iterate")+
       ggtitle("Fixed n")+
      theme(text = element_text(size=10))
g1

###################################################################

# make figure 3b

set.seed(12345)

rt1 <-matrix(0,5,50)
rt2 <-matrix(0,5,50)

n <- c(200,400,600,800,1000)
m <- 3
w <- 200

for(k in 1:5){
  for(j in 1:50){
    # make the weight matrix
    W1 <- matrix(1,n[k],n[k]) - diag(n[k])
    W2 <- diag(w,n[k])
    W <- W2
    for(i in 1:(3-1)){W<-cbind(W,W2)}
    WW<-W
    for(i in 1:(3-1)){W<-rbind(W,WW)}
    W<-W-diag(w,n[k]*3)
    W <- W+as.matrix(bdiag(rep(list(W1),3)))
    
    
    X0 <- matrix(rnorm(n[k]*2,5,1),n[k],2)
    X<-list()
    for(l in 1:3){ X[[l]]=jitter(X0,factor=1,amount=0)}
    DX <- list()
    for(l in 1:3){ DX[[l]]=as.matrix(dist(X[[l]]))
    DX[[l]]<-normDiss(DX[[l]],1)}
    
    # make the omnibus matrix
    d<-list()
    for(i in 1:3){
      d[[i]]<-DX[[i]]
      for(l in 1:(3-1)){d[[i]]<-cbind(d[[i]],DX[[i]])}
    }
    O<-Reduce(rbind,d)
    O<-O/2+t(O)/2
    
    # initialize at cMDS
    X0<-torgerson(sqrt(O),2)
    
    t<-proc.time();smds1 <- fJOFC(w,DX,2,verbose=FALSE,init=X0);runtime1<-proc.time()-t
    rt1[k,j]<-runtime1[[3]]/(smds1$niter)
    t<-proc.time();smds2 <- smacofSym(O,2,weightmat=W,init=X0,verbose=FALSE);runtime2<-proc.time()-t
    rt2[k,j]<-runtime2[[3]]/(smds2$niter)
    
  }}

means<-rep(0,10)
for(i in 1:5){means[i]=mean(rt1[i,])}
for(i in 6:10){means[i]=mean(rt2[i-5,])}
std<-rep(0,10)
for(i in 1:5){std[i]=2*sd(rt1[i,])/sqrt(50)}
for(i in 6:10){std[i]=2*sd(rt2[i-5,])/sqrt(50)}
x<-data.frame(cbind(means,std))
N<-c(n,n)
x<-cbind(x,N)
Method<-c(rep("fJOFC",5),rep("JOFC",5))
x<-cbind(x,Method)
g1<-ggplot(x,aes(colour=Method,x=N, y=means))+geom_errorbar(aes(ymin=means-std, ymax=means+std), width=0.2,size=1) +
  geom_line(size=1) +
  geom_point()+xlab("n") +
  ylab("Average runtime per iterate")+
  ggtitle("Fixed m")+
  theme(text = element_text(size=10))
g1

################################################################

# oos Experiment for Figure 4

set.seed(12345)
# insample runtime
rt1<-matrix(0,25,5)
# out of sample runtime
rt2<-matrix(0,25,5)

norms<-matrix(0,25,5)

n=c(200,300,400,500,600)
m<-c(10,15,20,25,30)

# left panel, vary m n=200
for(j in 1:25){
  for(k in 1:5){
    X0 <- matrix(rnorm(200*3,5,1),200,3)
    X<-list()
    for(l in 1:m[k]){ X[[l]]=jitter(X0,factor=1,amount=0)}
    DX <- list()
    for(l in 1:m[k]){ DX[[l]]=as.matrix(dist(X[[l]]))}
    # do the insample embedding
    t<-proc.time();smds1 <- fJOFC(100,DX,3,verbose=FALSE);runtime1<-proc.time()-t
    rt1[j,k]<-runtime1[[3]]
    Xhat<-list()
    for(i in 1:m[k]){Xhat[[i]]<-smds1$conf[((i-1)*200+1):(200*i),]}
    
    # out of sample embed
    Y<-list()
    for(l in 1:m[k]){ Y[[l]]=X[[l]][-200,]}
    DY <- list()
    for(l in 1:m[k]){ DY[[l]]=as.matrix(dist(Y[[l]]))}
    t<-proc.time();smds2 <- fJOFC(100,DY,3,verbose=FALSE);runtime1<-proc.time()-t
    Xinsample<-list()
    for(i in 1:m[k]){Xinsample[[i]]<-smds2$conf[((i-1)*199+1):(199*i),]}
    t<-proc.time();smds3<-oosfJOFC(DX,Xinsample,100);runtime2<-proc.time()-t
    rt2[j,k]<-runtime2[[3]]
    xhaty<-list()
    for(i in 1:m[k]){xhaty[[i]]<-rbind(Xinsample[[i]],smds3$Y[[i]])  }
    norms1<-rep(0,m[k])
    for(i in 1:m[k]){W<-procrustes(xhaty[[i]],Xhat[[i]])$W
    norms1[i]<-sqrt(sum( ( (xhaty[[i]]%*%W)[200,]-Xhat[[i]][200,])^2)) }
    norms[j,k]<-sum(norms1)
    }
}

# values in table 3 top row
means<-colMeans(norms)

Mean<-rep(0,10)
for(i in 1:5){Mean[i]<-mean(rt1[,i])}
for(i in 1:5){Mean[i+5]<-mean(rt2[,i])}
Std<-rep(0,10)
for(i in 1:5){Std[i]<-2*sd(rt1[,i])/5}
for(i in 1:5){Std[i+5]<-2*sd(rt2[,i])/5}
M<-c(m,m)
x<-data.frame(cbind(cbind(Mean,Std),M))
Method<-c(rep("fJOFC",5),rep("oos-fJOFC",5))
x<-cbind(x,Method)
g1<-ggplot(x,aes(colour=Method,x=M, y=Mean))+geom_errorbar(aes(ymin=Mean-Std, ymax=Mean+Std), width=0.2,size=1) +
  geom_line(size=1) +
  geom_point()+xlab("m") +
  ylab("Average runtime")+
  ggtitle("Fixed n=200")+
  theme(text = element_text(size=10))
g1

####
# fixed m, vary n

set.seed(12345)
# insample runtime
rt1<-matrix(0,25,5)
# out of sample runtime
rt2<-matrix(0,25,5)

norms<-matrix(0,25,5)

n=c(200,300,400,500,600)
m<-c(10,15,20,25,30)

# right panel, vary n m=10
for(j in 1:25){
  for(k in 1:5){
    X0 <- matrix(rnorm(n[k]*3,5,1),n[k],3)
    X<-list()
    for(l in 1:10){ X[[l]]=jitter(X0,factor=1,amount=0)}
    DX <- list()
    for(l in 1:10){ DX[[l]]=as.matrix(dist(X[[l]]))}
    # do the insample embedding
    t<-proc.time();smds1 <- fJOFC(100,DX,3,verbose=FALSE);runtime1<-proc.time()-t
    rt1[j,k]<-runtime1[[3]]
    Xhat<-list()
    for(i in 1:10){Xhat[[i]]<-smds1$conf[((i-1)*n[k]+1):(n[k]*i),]}
    
    # out of sample embed
    Y<-list()
    for(l in 1:10){ Y[[l]]=X[[l]][-n[k],]}
    DY <- list()
    for(l in 1:10){ DY[[l]]=as.matrix(dist(Y[[l]]))}
    t<-proc.time();smds2 <- fJOFC(100,DY,3,verbose=FALSE);runtime1<-proc.time()-t
    Xinsample<-list()
    for(i in 1:10){Xinsample[[i]]<-smds2$conf[((i-1)*(n[k]-1)+1):((n[k]-1)*i),]}
    t<-proc.time();smds3<-oosfJOFC(DX,Xinsample,100);runtime2<-proc.time()-t
    rt2[j,k]<-runtime2[[3]]
    xhaty<-list()
    for(i in 1:10){xhaty[[i]]<-rbind(Xinsample[[i]],smds3$Y[[i]])  }
    norms1<-rep(0,10)
    for(i in 1:10){W<-procrustes(xhaty[[i]],Xhat[[i]])$W
      norms1[i]<-sqrt(sum( ( (xhaty[[i]]%*%W)[n[k],]-Xhat[[i]][n[k],])^2)) }
    norms[j,k]<-sum(norms1)
  }
}

# values in table 3 bottom row
means<-colMeans(norms)

Mean<-rep(0,10)
for(i in 1:5){Mean[i]<-mean(rt1[,i])}
for(i in 1:5){Mean[i+5]<-mean(rt2[,i])}
Std<-rep(0,10)
for(i in 1:5){Std[i]<-2*sd(rt1[,i])/5}
for(i in 1:5){Std[i+5]<-2*sd(rt2[,i])/5}
N<-c(n,n)
x<-data.frame(cbind(cbind(Mean,Std),N))
Method<-c(rep("fJOFC",5),rep("oos-fJOFC",5))
x<-cbind(x,Method)
g1<-ggplot(x,aes(colour=Method,x=N, y=Mean))+geom_errorbar(aes(ymin=Mean-Std, ymax=Mean+Std), width=0.2,size=1) +
  geom_line(size=1) +
  geom_point()+xlab("n") +
  ylab("Average runtime")+
  ggtitle("Fixed m=10")+
  theme(text = element_text(size=10))
g1