
library(mnormt)
library(igraph)
library(cccd)
library(RANN)

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

#Power of the K-NN test for the plot in Figure 8(a) (normal location).

d=10
m=1000
n=800
s=c(1, 5, 20, 50, 100, 300, 600, 900)

mu0=rep(0, d)
sigma0=diag(rep(1, d))

iterations=100
iter=20

power1=matrix(nrow=iter, ncol=length(s)) 
rnull=matrix(nrow=iterations, ncol=length(s))
rp1=vector(length=iterations)

h=seq(0, 3, length.out=iter)


for(j in 1:length(s))
{
  for(i in 1:iterations)
  {
    
    points1<-rmnorm(m, mu0, sigma0)
    points2<-rmnorm(n, mu0, sigma0)
    
    pts<-rbind(points1,points2)
    aa1<-nn2(pts, points1, s[j]+1)$nn.idx[,-1]
    aa2<-nn2(pts, points2, s[j]+1)$nn.idx[,-1]
    
    
    rnull[i, j]<-length(which(c(aa1)>m))+length(which(c(aa2)<=m))
    
  }
}



for(k in 1:length(s))
{   
  for(j in 1:iter)
  {
    delta<-h[j]*rep(1, d)/((m+n)^(1/2))  
    mu0=rep(0, d)
    mu1=delta
    
    
    for(i in 1:iterations)
    {
      
      
      points1<-rmnorm(m, mu0, sigma0)
      points2<-rmnorm(n, mu1, sigma0)
      
      pts<-rbind(points1,points2)
      
      aa1<-nn2(pts, points1, s[k]+1)$nn.idx[,-1]
      aa2<-nn2(pts, points2, s[k]+1)$nn.idx[,-1]
      
      rp1[i]<-length(which(c(aa1)>m))+length(which(c(aa2)<=m))
      
    }
    
    power1[j, k]<-length(which(rp1< quantile(rnull[,k], 0.05)))/iterations
    print(power1[j, k])
  }
}

write.table(as.matrix(power1, nrow=iter, ncol=length(s)), file="NN_Normal.txt") 


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

#Power of the K-NN test for the plot in Figure 8(b) (normal scale).


library(mnormt)
library(igraph)
library(cccd)
library(RANN)

d=10
m=1000
n=800
s=c(1, 5, 20, 50, 100, 300, 600, 900)

mu0=rep(0, d)
sigma0=diag(rep(1, d))

iterations=100
iter=20

power1=matrix(nrow=iter, ncol=length(s)) 
rnull=matrix(nrow=iterations, ncol=length(s))
rp1=vector(length=iterations)

h=seq(0, 3, length.out=iter)


for(j in 1:length(s))
{
  for(i in 1:iterations)
  {
    
    points1<-exp(rmnorm(m, mu0, sigma0))
    points2<-exp(rmnorm(n, mu0, sigma0))
    
    pts<-rbind(points1,points2)
    aa1<-nn2(pts, points1, s[j]+1)$nn.idx[,-1]
    aa2<-nn2(pts, points2, s[j]+1)$nn.idx[,-1]
    
    
    rnull[i, j]<-length(which(c(aa1)>m))+length(which(c(aa2)<=m))
    
  }
}


for(k in 1:length(s))
{   
  for(j in 1:iter)
  {
     
    delta1<-h[j]/(m+n)^(1/2)  
    sigma1=sigma0+delta1*diag(rep(1, d))
    
    for(i in 1:iterations)
    {
      
      
      points1<-exp(rmnorm(m, mu0, sigma0))
      points2<-exp(rmnorm(n, mu0, sigma1))
      
      pts<-rbind(points1,points2)
      
      aa1<-nn2(pts, points1, s[k]+1)$nn.idx[,-1]
      aa2<-nn2(pts, points2, s[k]+1)$nn.idx[,-1]
      
      rp1[i]<-length(which(c(aa1)>m))+length(which(c(aa2)<=m))
      
    }
    
    power1[j, k]<-length(which(rp1< quantile(rnull[,k], 0.05)))/iterations
    print(power1[j, k])
  }
}

write.table(as.matrix(power1, nrow=iter, ncol=length(s)), file="NN_Scale.txt") 



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


#Power of the K-NN test for the plot in Figure 5(b) (lognormal location with dependent covariance matrix)


d=10
m=1000
n=800
s=c(1, 5, 20, 50, 100, 300, 600, 900)

mu0=rep(0, d)
sigma0=rep(1, d)%*%t(rep(1, d))+diag(rep(1,d))


iterations=100
iter=20

power1=matrix(nrow=iter, ncol=length(s)) 
rnull=matrix(nrow=iterations, ncol=length(s))
rp1=vector(length=iterations)

h=seq(0, 3, length.out=iter)


for(j in 1:length(s))
{
for(i in 1:iterations)
{
  
  points1<-exp(rmnorm(m, mu0, sigma0))
  points2<-exp(rmnorm(n, mu0, sigma0))
  
  pts<-rbind(points1,points2)
  aa1<-nn2(pts, points1, s[j]+1)$nn.idx[,-1]
  aa2<-nn2(pts, points2, s[j]+1)$nn.idx[,-1]
  
  
  rnull[i, j]<-length(which(c(aa1)>m))+length(which(c(aa2)<=m))
  
}
}



for(k in 1:length(s))
{   
for(j in 1:iter)
{
  delta<-h[j]*rep(1, d)/((m+n)^(1/2))  
  mu0=rep(0, d)
  mu1=delta
  
  
  for(i in 1:iterations)
  {
    
    
    points1<-exp(rmnorm(m, mu0, sigma0))
    points2<-exp(rmnorm(n, mu1, sigma0))
    
    pts<-rbind(points1,points2)
    
    aa1<-nn2(pts, points1, s[k]+1)$nn.idx[,-1]
    aa2<-nn2(pts, points2, s[k]+1)$nn.idx[,-1]
    
    rp1[i]<-length(which(c(aa1)>m))+length(which(c(aa2)<=m))
    
  }
  
  power1[j, k]<-length(which(rp1< quantile(rnull[,k], 0.05)))/iterations
  print(power1[j, k])
}
}

write.table(as.matrix(power1, nrow=iter, ncol=length(s)), file="NN_Lognormal_Dependent.txt") 




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


#Power of the K-NN test for the plot in Figure 5(a) (lognormal location with independent covariance matrix)



library(mnormt)
library(igraph)
library(cccd)
library(RANN)

d=10
m=1000
n=800
s=c(1, 5, 20, 50, 100, 300, 600, 900)

mu0=rep(0, d)
sigma0=diag(rep(1, d))

iterations=100
iter=20

power1=matrix(nrow=iter, ncol=length(s)) 
rnull=matrix(nrow=iterations, ncol=length(s))
rp1=vector(length=iterations)

h=seq(0, 3, length.out=iter)


for(j in 1:length(s))
{
  for(i in 1:iterations)
  {
    
    points1<-exp(rmnorm(m, mu0, sigma0))
    points2<-exp(rmnorm(n, mu0, sigma0))
    
    pts<-rbind(points1,points2)
    aa1<-nn2(pts, points1, s[j]+1)$nn.idx[,-1]
    aa2<-nn2(pts, points2, s[j]+1)$nn.idx[,-1]
    
    
    rnull[i, j]<-length(which(c(aa1)>m))+length(which(c(aa2)<=m))
    
  }
}



for(k in 1:length(s))
{   
  for(j in 1:iter)
  {
    delta<-h[j]*rep(1, d)/((m+n)^(1/2))  
    mu0=rep(0, d)
    mu1=delta
    
    
    for(i in 1:iterations)
    {
      
      
      points1<-exp(rmnorm(m, mu0, sigma0))
      points2<-exp(rmnorm(n, mu1, sigma0))
      
      pts<-rbind(points1,points2)
      
      aa1<-nn2(pts, points1, s[k]+1)$nn.idx[,-1]
      aa2<-nn2(pts, points2, s[k]+1)$nn.idx[,-1]
      
      rp1[i]<-length(which(c(aa1)>m))+length(which(c(aa2)<=m))
      
    }
    
    power1[j, k]<-length(which(rp1< quantile(rnull[,k], 0.05)))/iterations
    print(power1[j, k])
  }
}

write.table(as.matrix(power1, nrow=iter, ncol=length(s)), file="NN_Lognormal_Independent.txt") 



