##Tri distribution##########################
tri<-function(n,p1,p2){
	a<-runif(n)
	b1<-length(a[a<p1])
	b2<-length(a[a>=p1&a<p1+p2])
	b3<-length(a[a>=p1+p2])
	return(c(b1,b2,b3))
}
## multinorm distribution###################
trinorm<-function(n,a,b){
c<-eigen(b)
d<-c$values
e<-c$vectors
f<-e%*%diag(sqrt(d))%*%t(e)
g<-NULL
for(i in 1:n){
h<-rnorm(3,mean=a)
i<-f%*%as.matrix(h)
i<-as.vector(i)
g<-c(g,i)
}
j<-matrix(data=g,nrow=3,ncol=n)
return(t(j))
}

## CA test##################
ca<-function(x,r1,r2,s1,s2,r,s){
n<-r+s
n1<-r1+s1
n2<-r2+s2
p1<-r1/r
p2<-r2/r
q1<-s1/s
q2<-s2/s
u<-(x*(s*r1-r*s1)+s*r2-r*s2)/n
evaru<-(r*s^2*((x^2*p1+p2)-(x*p1+p2)^2)/n^2)+(s*r^2*((x^2*q1+q2)-(x*q1+q2)^2)/n^2)
return(u/sqrt(evaru))
}
## delta of DC#############
delta<-function(x,r,s,ep1,ep2,eprg1,eprg2){
n<-r+s
a<-(n*r/s)^0.5*(x*ep1+ep2-x*eprg1-eprg2)
b<-x^2*eprg1+eprg2-(x*eprg1+eprg2)^2
c<-a/sqrt(b)
return(c)
}

##generate data##############
gd<-function(rm,sm,f,p,m,f0,f1,f2){
pm<-rbeta(m,(1-f)*p/f,(1-f)*(1-p)/f)
pg0<-f*(1-pm)+(1-f)*(1-pm)^2
pg1<-2*(1-f)*pm*(1-pm)
pg2<-f*(pm)+(1-f)*(pm)^2
km<-pg0*f0+pg1*f1+pg2*f2
pg0ca<-pg0*f0/km
pg1ca<-pg1*f1/km
pg2ca<-pg2*f2/km
pg0co<-pg0*(1-f0)/(1-km)
pg1co<-pg1*(1-f1)/(1-km)
pg2co<-pg2*(1-f2)/(1-km)
ca<-rep(0,3)
co<-rep(0,3)
for(i in 1:m){
a<-tri(rm[i],pg0ca[i],pg1ca[i])
b<-tri(sm[i],pg0co[i],pg1co[i])
ca<-ca+a
co<-co+b
}
result<-matrix(c(ca[1],co[1],ca[2],co[2],ca[3],co[3]),nrow=2,ncol=3)
return(result)
}
##generate the percentile of Tmax########
percentile<-function(n,rou1,rou2,rou3){
	a<-trinorm(n,c(0,0,0),matrix(c(1,rou1,rou2,rou1,1,rou3,rou2,rou3,1),nrow=3,ncol=3))
	b<-rep(0,n)
	for(i in 1:n){
		b[i]=max(a[i,1]^2,a[i,2]^2,a[i,3]^2)
	}
	c<-sort(b)
	result<-c[0.95*n+1]
	return(result)
}
	
##correlation################
cc<-function(ep0,ep1,ep2){
a1<-ep0*(ep1+2*ep2)
b1<-sqrt(ep0*(1-ep0))*sqrt(ep0*(ep1+2*ep2)+ep2*(ep1+2*ep0))
a2<-ep0*ep2
b2<-sqrt(ep0*(1-ep0))*sqrt(ep2*(1-ep2))
a3<-ep2*(ep2+2*ep0)
b3<-sqrt(ep2*(1-ep2))*sqrt(ep0*(ep1+2*ep2)+ep2*(ep1+2*ep0))
c1<-a1/b1
c2<-a2/b2
c3<-a3/b3
return(c(c1,c2,c3))		
}
## transitional parameter##################
ca1no<-NULL
ca2no<-NULL
ca3no<-NULL
qmax<-NULL
ca1nonl<-NULL
ca2nonl<-NULL
ca3nonl<-NULL
ep1nl<-NULL
ep2nl<-NULL
eprg0<-NULL
eprg1<-NULL
eprg2<-NULL
gamma1<-NULL
gamma2<-NULL
gamma3<-NULL
delta1<-NULL
delta2<-NULL
delta3<-NULL
z1gc<-NULL
z2gc<-NULL
z3gc<-NULL
z1dc<-NULL
z2dc<-NULL
z3dc<-NULL
tmaxgc<-NULL
tmertgc<-NULL
tpearsongc<-NULL
tmaxdc<-NULL
tmertdc<-NULL
tpearsondc<-NULL
tmaxno<-NULL
tmertno<-NULL
tpearsonno<-NULL
result<-NULL
##main programme########################
m1<-10000
## m1 is the replication times
m<-3
## m is the subpopulation
rm<-c(200,100,200)
## rm is the case of each subpopulation
sm<-c(100,300,100)
## sm is the control of each subpopulation
f<-0.005
## f is the inbreeding coefficient
p<-0.1
## p is the allel frequency
f0<-0.02
f1<-0.02
f2<-0.02
## fi = P(case|Gi) for i = 0, 1, 2 with f0 > 0.
r<-sum(rm)
## r is the number of case
s<-sum(sm)
## s is the number of control
n<-r+s
for(i in 1:m1){
	matrix1<-gd(rm,sm,f,p,m,f0,f1,f2)
        if((matrix1[1,1]<0.5)||(matrix1[1,2]<0.5)||(matrix1[1,3]<0.5)||(matrix1[2,1]<0.5)||(matrix1[2,2]<0.5)||(matrix1[2,3]<0.5)){
        matrix1<-matrix(c(matrix1[1,1]+0.5,matrix1[2,1]+0.5,matrix1[1,2]+0.5,matrix1[2,2]+0.5,matrix1[1,3]+0.5,matrix1[2,3]+0.5),nrow=2,ncol=3)
}
	r0<-matrix1[1,1]
	r1<-matrix1[1,2]
	r2<-matrix1[1,3]
	s0<-matrix1[2,1]
	s1<-matrix1[2,2]
	s2<-matrix1[2,3]
	n0<-r0+s0
	n1<-r1+s1
	n2<-r2+s2
	ca1no[i]<-ca(0,r1,r2,s1,s2,r,s)
	## ca1no is the CA trend test when x is 0 without corrected
	ca2no[i]<-ca(0.5,r1,r2,s1,s2,r,s)
	## ca2no is the CA trend test when x is 0.5 without corrected
	ca3no[i]<-ca(1,r1,r2,s1,s2,r,s)
	## ca3no is the CA trend test when x is 1 without corrected
	
	l<-50
	## l is the number of null loci
	pnl<-runif(l,min=p-0.05,max=p+0.05)
	for(j in 1:l){
	    matrix2<-gd(rm,sm,f,pnl[j],m,f0,f0,f0)
            if((matrix2[1,1]<0.5)||(matrix2[1,2]<0.5)||(matrix2[1,3]<0.5)||(matrix2[2,1]<0.5)||(matrix2[2,2]<0.5)||(matrix2[2,3]<0.5)){
        matrix2<-matrix(c(matrix2[1,1]+0.5,matrix2[2,1]+0.5,matrix2[1,2]+0.5,matrix2[2,2]+0.5,matrix2[1,3]+0.5,matrix2[2,3]+0.5),nrow=2,ncol=3)
}
	r0nl<-matrix2[1,1]
	r1nl<-matrix2[1,2]
	r2nl<-matrix2[1,3]
	s0nl<-matrix2[2,1]
        s1nl<-matrix2[2,2]
	s2nl<-matrix2[2,3]
	n0nl<-r0nl+s0nl
	n1nl<-r1nl+s1nl
	n2nl<-r2nl+s2nl
	ca1nonl[j]<-ca(0,r1nl,r2nl,s1nl,s2nl,r,s)
	ca2nonl[j]<-ca(0.5,r1nl,r2nl,s1nl,s2nl,r,s)
	ca3nonl[j]<-ca(1,r1nl,r2nl,s1nl,s2nl,r,s)
        ep1nl[j]<-r1nl/r
	ep2nl[j]<-r2nl/r
        eprg0[j]<-n0nl/n
        eprg1[j]<-n1nl/n
        eprg2[j]<-n2nl/n
	}
	gamma1[i]<-mean(ca1nonl^2)
	gamma2[i]<-mean(ca2nonl^2)
	gamma3[i]<-mean(ca3nonl^2)
	ep1nlm<-mean(ep1nl)
	ep2nlm<-mean(ep2nl)
        eprg0m<-mean(eprg0)
	eprg1m<-mean(eprg1)
	eprg2m<-mean(eprg2)
        rou1<-cc(eprg0m,eprg1m,eprg2m)[1]
	rou2<-cc(eprg0m,eprg1m,eprg2m)[2]
	rou3<-cc(eprg0m,eprg1m,eprg2m)[3]
	m2<-100
	## m2 is the replication times to deceide the percentile of Tmax
	qmax[i]<-percentile(m2,rou1,rou2,rou3) 
	delta1[i]<-delta(0,r,s,ep1nlm,ep2nlm,eprg1m,eprg2m)
	delta2[i]<-delta(0.5,r,s,ep1nlm,ep2nlm,eprg1m,eprg2m)
	delta3[i]<-delta(1,r,s,ep1nlm,ep2nlm,eprg1m,eprg2m)
	z1gc[i]<-sign(ca1no[i])*sqrt(ca1no[i]^2/gamma1[i])
	z2gc[i]<-sign(ca2no[i])*sqrt(ca2no[i]^2/gamma2[i])
	z3gc[i]<-sign(ca3no[i])*sqrt(ca3no[i]^2/gamma3[i])
	z1dc[i]<-ca1no[i]-delta1[i]
	z2dc[i]<-ca2no[i]-delta2[i]
	z3dc[i]<-ca3no[i]-delta3[i]
	tmaxgc[i]<-max(z1gc[i]^2,z2gc[i]^2,z3gc[i]^2)
	tmertgc[i]<-(z1gc[i]+z3gc[i])^2/(2*(1+rou2))
	tpearsongc[i]<-(z1gc[i]^2+z3gc[i]^2-2*rou2*z1gc[i]*z3gc[i])/(1-rou2^2)
	tmaxdc[i]<-max(z1dc[i]^2,z2dc[i]^2,z3dc[i]^2)
	tmertdc[i]<-(z1dc[i]+z3dc[i])^2/(2*(1+rou2))
	tpearsondc[i]<-(z1dc[i]^2+z3dc[i]^2-2*rou2*z1dc[i]*z3dc[i])/(1-rou2^2)
	tmaxno[i]<-max(ca1no[i]^2,ca2no[i]^2,ca3no[i]^2)
	tmertno[i]<-(ca1no[i]+ca3no[i])^2/(2*(1+rou2))
	tpearsonno[i]<-(ca1no[i]^2+ca3no[i]^2-2*rou2*ca1no[i]*ca3no[i])/(1-rou2^2)	
}
result[1]<-mean(ca1no^2>qchisq(0.95,1))
## result[1] is the Type I error(power) of T0 with out corrected
result[2]<-mean(ca2no^2>qchisq(0.95,1))
## result[2] is the Type I error(power) of T0.5 with out corrected
result[3]<-mean(ca3no^2>qchisq(0.95,1))
## result[3] is the Type I error(power) of T1 with out corrected
result[4]<-mean(tmertno>qchisq(0.95,1))
## result[4] is the Type I error(power) of MERT without corrected
result[5]<-mean(tpearsonno>qchisq(0.95,2))
##result[5] is the Type I error(power) of Pearson without corrected
result[6]<-mean(tmaxno>qmax)
## result[6] is the Type I error(power) of MAX without corrected
result[7]<-mean(z1dc^2>qchisq(0.95,1))
## result[7] is the Type I error(power) of T0 with DC corrected
result[8]<-mean(z2dc^2>qchisq(0.95,1))
## result[8] is the Type I error(power) of T0.5 with DC corrected
result[9]<-mean(z3dc^2>qchisq(0.95,1))
## result[9] is the Type I error(power) of T1 with DC corrected
result[10]<-mean(tmertdc>qchisq(0.95,1))
## result[10] is the Type I error(power) of MERT with DC corrected
result[11]<-mean(tpearsondc>qchisq(0.95,2))
## result[11] is the Type I error(power) of Pearson with DC corrected
result[12]<-mean(tmaxdc>qmax)
## result[12] is the Type I error(power) of MAX with DC corrected
result[13]<-mean(z1gc^2>qchisq(0.95,1))
## result[13] is the Type I error(power) of T0 with GC corrected
result[14]<-mean(z2gc^2>qchisq(0.95,1))
## result[14] is the Type I error(power) of T0.5 with GC corrected
result[15]<-mean(z3gc^2>qchisq(0.95,1))
## result[15] is the Type I error(power) of T1 with GC corrected
result[16]<-mean(tmertgc>qchisq(0.95,1))
## result[16] is the Type I error(power) of MERT with GC corrected
result[17]<-mean(tpearsongc>qchisq(0.95,2))
##result[17] is the Type I error(power) of Pearson with GC corrected
result[18]<-mean(tmaxgc>qmax)
## result[18] is the Type I error(power) of MAX with GC corrected
print(result)