# This is written in R
# The input file "dat" is a matrix
# The first column of dat are the phenotype values
# Each of the other columns are the coded genotypes at a SNP
# The number of SNPs is supposed to be greater than 3
# The output is the p-value calculated by 10,000 permutations
# Before running the program, the following should be included
# library(mvtnorm)
# library(wavethresh)
# library(EbayesThresh)

function(dat){
np=4
x=dat[,-1]
y=dat[,1]
n=dim(x)[1]
m=dim(x)[2]

if(m<4){
print("The number of SNPs must be greater than 3")
}

mind=(log(m)/log(2))
if(floor(log(m)/log(2))==mind){

Ans=matrix(0,n,m)
for (i in 1:n){
     x1=as.vector(x[i,1:m])
     x2=wd(x1,family="DaubLeAsymm",filter.number=8)
     if(sum(x1)==0) x3=x2
     else{
     x3=ebayesthresh.wavelet(x2,vscale="level")
     }
     x4=wr(x3)
    Ans[i,]=x4
 }
Xt=Ans
U=numeric(m)
V=numeric(m)
    for ( j in 1:m){
    xmean=mean(Xt[,j])
    U=c(U,sum(y*(Xt[,j]-xmean)))
    V=c(V,((n/(4*(n-1)))*sum((Xt[,j]-xmean)^2)))
    }
    sumv=sum(V)
    if(sumv>0) t=((sum(U))^2)/(sumv)
    else t=0
tmax=t
nsim=(10^np)
tx=numeric(nsim)
ic=0
for( i in i:nsim){
U=numeric(m)
V=numeric(m)
y1=sample(y)
    for ( j in 1:m){
    xmean=mean(Xt[,j])
    U=c(U,sum(y1*(Xt[,j]-xmean)))
    V=c(V,(n/(4*(n-1))*sum((Xt[,j]-xmean)^2)))
    }
    sumv=sum(V)
    if(sumv>0) tx[i]=((sum(U))^2)/sumv
    else tx[i]=0
if(tx[i]>tmax) ic=(ic+1)
}
pval=(ic/nsim)
return(pval)
}

else{
m1=floor(log(m)/log(2))
m1=2^m1
Xt=matrix(0,n,m)
Ans=matrix(0,n,m1)
for (i in 1:n){
     x1=as.vector(x[i,1:m1])
     x2=wd(x1,family="DaubLeAsymm",filter.number=8)
     if(sum(x1)==0) x3=x2
     else{
     x3=ebayesthresh.wavelet(x2,vscale="level")
     }
     x4=wr(x3)
    Ans[i,]=x4
 }
Xt[,(1:m1)]=Ans
Xt[,((m1+1):m)]=x[,((m1+1):m)]
U=numeric(m)
V=numeric(m)
    for ( j in 1:m){
    xmean=mean(Xt[,j])
    U=c(U,sum(y*(Xt[,j]-xmean)))
    V=c(V,((n/(4*(n-1)))*sum((Xt[,j]-xmean)^2)))
    }
    sumv=sum(V)
    if(sumv>0) t=((sum(U))^2)/(sumv)
    else t=0
tmax=t
nsim=(10^np)
tx=numeric(nsim)
ic=0
for( i in i:nsim){
U=numeric(m)
V=numeric(m)
y1=sample(y)
    for ( j in 1:m){
    xmean=mean(Xt[,j])
    U=c(U,sum(y1*(Xt[,j]-xmean)))
    V=c(V,(n/(4*(n-1))*sum((Xt[,j]-xmean)^2)))
    }
    sumv=sum(V)
    if(sumv>0) tx[i]=((sum(U))^2)/sumv
    else tx[i]=0
if(tx[i]>tmax) ic=(ic+1)
}
pval=(ic/nsim)
return(pval)
}

}