# 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