N=c(3,3,2)
L=length(N)
K=3


library(base)

#package for analysing Markov chains
library(coda)
library(boa)

#input the Markov basis for the null model
source("basis_no3way.txt")

#input the data tables, in our case 10 chose 2 tables corresponding to all possible pairs of SNPs among the 10 SNPs with the lowest marginal p-values.
load("ListOfTables.RDATA")

#we ran 3 Markov chains with 40000 iterations each and used a burn-in of 10000 iterations each
n.chains=3
n.iter=40000
burnin=10000

#we save the chi-square statistic of each iteration
sims=array(NA,c(n.iter,n.chains,length(List.tables),1))
dimnames(sims)=list(NULL,NULL,NULL,c("S"))

#out will contain the chi-square statistic of the given data table, its p-value before Bonferroni correction, and the value 1 for significant and 0 for not significant using a 5% cutoff (after Bonferroni correction).
out=rep(2,3*length(List.tables))
dim(out)=c(3,length(List.tables))
for (k in 1:length(List.tables)){

	#if table is not of dimension 3x3x2, we fill in zero cell counts to make it full dimensional
	D=List.tables[[k]]
	if (dim(D)[1] != 3 | dim(D)[2] != 3) {
		D1=rep(0,18)
		dim(D1)=c(3,3,2)
		N1=names(D[1,,1])
		N2=names(D[,1,1])
		if (dim(D)[1]==1){N1="0"}
		if (dim(D)[2]==1){N2="0"}
		d=c("0", "1", "2")
		N1=c(N1,rep(4,3-length(N1)))
		N2=c(N2,rep(4,3-length(N2)))
		a=which(N1==d)
		b=which(N2==d)
		D1[a,b,1:2]=D
		D=D1
	}
	
	#fit table to the no 3-way interaction model by iterative proportional fitting
	c=combn(1:L,K-1)
	C=list()
	for (i in 1:length(c[1,])){
		C=c(C,list(c[,i]))
	}

	pf=loglin(D,C,fit=TRUE)
	square=(D-pf$fit)^2
	indices=which(square==0)
	pf.nozeros=pf$fit
	for (i in indices){
  		pf.nozeros[i]=1
	}
	matrix=square/pf.nozeros
	S=sum(matrix)

	#chi-square statistic of data table
	S.new=S

	fun=function(x){
		sum(log(seq(1,x)))
	}

	#mcmc: we sample an element of the Markov basis, add or subtract it to the table from the previous 	iteration and accept the new table (given all entries are non-negative) with probability r1.
	for (m in 1:n.chains){
		for (s in 1:n.iter){
			rf=sample(length(F),1)
			t=sample(c(1,-1),1)
			D.new=D+t*F[[rf]]
			if (all(D.new >=0)){
				b1=log(runif(1))
				D.log=D
				D.new.log=D.new
				D.log[D==0]=1
				D.new.log[D.new==0]=1
				D.llog=as.matrix(D.log[D.log!=D.new.log])
				D.new.llog=as.matrix(D.new.log[D.new.log!=D.log])
				r1=sum(apply(D.llog,c(1,2),fun))-sum(apply(D.new.llog,c(1,2),fun))
				if (b1<r1){
					D=D.new
					S.new=sum((D-pf$fit)^2/pf.nozeros)
				}
			}
			sims[s,m,k,]=c(S.new)
		}
	}
	B=sapply(1:n.chains,function(i)assign(paste("A",i,sep=""),sims[,i,k,],envir=.GlobalEnv))

	
	A1=as.mcmc(A1)
	A2=as.mcmc(A2)
	A3=as.mcmc(A3)

	x1=mcmc.list(A1,A2,A3)
	y1=window(x1,burnin,n.iter,1)

	v=seq(0.0001, 1, by=.0001)
	v1=summary(y1,quantiles=v)$quantiles
	if (S>v1[length(v)]){
		p.value=0
	}
	if (S<=v1[length(v)]){
		d=abs(v1-S)
		e=which(d==min(d))
		p.value=1-e[1]/length(d)
	}
	if (p.value<0.05/length(List.tables)){
		out[,k]=c(S,p.value,1)
		if (S<=summary(y1,quantiles=1-0.05/length(List.tables))$quantiles){
			out[,k]=c(S,p.value,0)
		}
	}
	if (p.value>=0.05/length(List.tables)){
		out[,k]=c(S,p.value,0)
	}	
}