library(BiasedUrn)
generate_genotype=function(sample, M, K, L)
{
	nk=K^2 # nk is the total number of grids
	c=sample/nk
	nnode=rep(c,nk)
	node=matrix(0,nk,4*c)
	cnode=matrix(0,2*sample,3)
	nn=0
	for (i in 1:nk)
	{
		for(j in 1:c)	{nn=nn+1;	node[i,j]=nn}
	}

	nb=rep(0,nk)
	nbb=matrix(0,nk,4)	

	for(m in 1:nk)
	{
		i=floor((m-1)/K+1)
		j=m-(i-1)*K
		nn=0
		if(i>1) {nn=nn+1;nbb[m,nn]=(i-2)*K+j}
		if(i<K) {nn=nn+1;nbb[m,nn]=(i)*K+j}
		if(j>1) {nn=nn+1;nbb[m,nn]=(i-1)*K+j-1}
		if(j<K) {nn=nn+1;nbb[m,nn]=(i-1)*K+j+1}
		nb[m]=nn
	}
	
	tc=tm=rep(0,nk)
	nn=0
	t=0
	nnn=0
	m0=1
	while( (nn<=(sample-3)) & (nnn<= (4*sample)) )
	{
		nnn=nnn+1
		tc=nnode*(nnode-1)/2/m0
		tcc=sum(tc)
		tm=M*nnode*nb/2/m0
		tmm=sum(tm)
		t=t+exp(tcc+tmm)
		a=tcc/(tcc+tmm)
		if(runif(1)<a) 
		{
			nn=nn+1
			cnode[sample+nn,3]=t
			m=sample(c(1:nk),1,replace=F,prob=tc/tcc)
			mm=sample(c(1:nnode[m]),2,replace=F)
			nnode[m]=nnode[m]-1
			cnode[sample+nn,1]=node[m,mm[1]]
			cnode[sample+nn,2]=node[m,mm[2]]
			m1=min(mm)
			m2=max(mm)
			node[m,m1]=nn+sample
			for(i in m2:(nnode[m]+1))  node[m,i]=node[m,(i+1)]	
		}
		else
		{
			m=sample(c(1:nk),1,replace=F,prob=tm/tmm)
			mm=sample(c(1:nnode[m]),1,replace=F)
			nnode[m]=nnode[m]-1
			mnode=node[m,mm]
			for(i in mm:(nnode[m]+1))  node[m,i]=node[m,(i+1)]
			nb0=sample(c(1:nb[m]),1,replace=F)
			m0=nbb[m,nb0]
			nnode[m0]=nnode[m0]+1
			node[m0,nnode[m0]]=mnode
		}
	}
pnode=which(cnode[,3]>0)
pt=cnode[pnode,3]
n0=sample(pnode,2*L,replace=T,prob=pt/sum(pt))

geno=matrix(0,sample,L)
for(i in 1:L)
	{
	geno1=geno2=rep(0,sample)	
	geno[,i]=count_node(sample, cnode, n0[2*i-1], geno1)+count_node(sample, cnode, n0[2*i], geno2)
	}
geno
}


##############################################
count_node=function(sample, cnode, n0, geno)
{
	if((n0<= sample)&(n0>0)) 
	{
		geno[n0]=1	
	}
	else if(n0>sample)
	{
		geno=count_node(sample, cnode, cnode[n0,1], geno)
		geno=count_node(sample, cnode, cnode[n0,2], geno)
	}
	geno
}
#################################################
trait_null0_qual=function(sample, K, model, RR,mm,sigma,percent)
{
	nk=K^2 # nk is the total number of grids
	c=sample/nk
	map=rep(0,sample)
	
	nn=0
	for (i in 1:nk)
	{
		for(j in 1:c)	{nn=nn+1;	map[nn]=i}
	}
	if(model==0) y=rnorm(sample)
	i0=j0=6
	if(model==1) #
	{
		y=rnorm(sample)
		for(i in i0:(i0+mm)) 
		{
			for(j in j0:(j0+mm))
			{
				nn=(i-1)*K+(j)
				a=which(map==nn)
				y[a]=y[a]+RR
			}
		}
	}
	
	if(model==2)
	{
		y=rnorm(sample)
		r=0.4*RR
		for(i in 1:K)
		{
			for(j in 1:K)
			{
				mu=r*exp(-((i-i0)^2+(j-j0)^2)/sigma**2/2)
				nn=(i-1)*K+(j)
				a=which(map==nn)
				y[a]=y[a]+mu
			}
		}
	}
nn=floor(sample*(1-percent))
aa=sort(y)[nn]
y=(y>aa)
}
#####################################################
uncorrect_single=function(y,x)
{
	n=length(y)
	#y1=rnorm(n)
	pv=1-pchisq((n-1)*cor(y,x)^2,1)
	#pv1=1-pchisq((n-1)*cor(y1,x)^2,1)
	#pvalue=cbind(pv,pv1)
}
######################################################
GC=function(pv)
{
     T=qchisq(1-pv,1)
     lambda=median(T)/0.456
     pvalue=1-pchisq(T/lambda,1)
}
#################################################
pc_linear_single=function(y, x, pc)
{
	n=length(y)
	x=cbind(y,x)
	x=lm(x~pc)$resid
	pv=1-pchisq((n-1)*cor(x[,1],x[,-1])^2,1)
}
######################################
pc_per_single=function(y, x, pc, numper)
{
	n=length(y)
	m=ncol(x)
	#a=median(y)
	#y=(y>a)
	ncase=sum(y)
	T=cor(y,x)**2
	Tper1=Tper2=matrix(0,numper,m)
	model=glm (y ~pc, family= binomial())
	dodds=exp (model$linear.predictors)
	m1=rep(1, n)
	perm=rMFNCHypergeo(numper, m1, ncase, dodds)
	for(i in 1:numper)
	{
		#print(paste("per=",i))
		#y1=sample(y)
		y2=perm[,i]
		#Tper1[i,]=cor(y1,x)**2
		Tper2[i,]=cor(y2,x)**2
	}
	pv1=pv2=rep(0,m)
	for(i in 1:m)
	{
		#pv1[i]=sum(Tper1[,i]>T[i])/numper
		#pv2[i]=sum(Tper2[,i]>T[i])/numper
		#pv1[i]=sum(Tper1[,i]>T[i])/numper+sum(Tper1[,i]==T[i])/numper/2
		pv2[i]=sum(Tper2[,i]>T[i])/numper+sum(Tper2[,i]==T[i])/numper/2
	}
	pv2#=cbind(pv1,pv2)
}
#############################################################
new_choose_h=function(y, x, pc, k)
{
	ee=0.000001
	n=length(y)
	m=ncol(x)
	xx=cbind(y,x)
	for(i in 1:k) pc[,i]=(pc[,i]-min(pc[,i]))/(max(pc[,i])-min(pc[,i]))
	ww=matrix(0,n,n)
	nn=(c(1:m)-0.5)/m
   mm=18	
   bb=2**c(0:mm)
	hh=c(2**c(1:4),rep(1,mm+1)/bb)
   hh=hh**(2/(5+k))
	pvalue=matrix(0,m,length(hh))
   pv=rep(0,length(hh))
	iter=0
	for(h in hh)
	{
	iter=iter+1
	for(i in 1:n) 
	{
		w=rep(1,n)
		for(j in 1:k)
		{
			t=(pc[,j]-pc[i,j])/h
			t=(1-t^2)^2*(t^2<1)
			w=w*t
		}
	ww[i,]=w/sum(w)
	}
	x=xx
	x=x-ww%*%x
	aa=(apply(x[,-1],2,var)+ee)*(var(x[,1]+ee))
	pv0=1-pchisq((n-1)*cov(x[,1],x[,-1])^2/aa,1)
 	pv0=sort(pv0)
	pvalue[,iter]=pv0
    pv[iter]=max(abs(pv0-nn))
	}
ind=which(pv==min(pv))[1]
h=hh[ind]
}
######################################################
pc_nonp_single=function(y, x, pc, k, h)
{
	n=length(y)
	x=cbind(y,x)
	#h=h^(1/k)
	for(i in 1:k) pc[,i]=(pc[,i]-min(pc[,i]))/(max(pc[,i])-min(pc[,i]))
	ww=matrix(0,n,n)
	for(i in 1:n) 
	{
		w=rep(1,n)
		for(j in 1:k)
		{
			t=(pc[,j]-pc[i,j])/h
			t=(1-t^2)^2*(t^2<1)
			w=w*t
		}
	ww[i,]=w/sum(w)
	}
	x=x-ww%*%x
	pv=1-pchisq((n-1)*cor(x[,1],x[,-1])^2,1)
}
###########################################################
trait_null0=function(sample, K, model, RR,mm)
{
	nk=K^2 # nk is the total number of grids
	c=sample/nk
	map=rep(0,sample)
	
	nn=0
	for (i in 1:nk)
	{
		for(j in 1:c)	{nn=nn+1;	map[nn]=i}
	}
	if(model==0)y=rnorm(sample)
	i0=j0=6
	if(model==1) #
	{
		y=rnorm(sample)
		for(i in i0:(i0+mm)) 
		{
			for(j in j0:(j0+mm))
			{
				nn=(i-1)*K+(j)
				a=which(map==nn)
				y[a]=y[a]+RR
			}
		}
	}
	
	if(model==2)
	{
		y=rnorm(sample)
		r=0.4*RR
		for(i in 1:K)
		{
			for(j in 1:K)
			{
				mu=r*exp(-((i-i0)^2+(j-j0)^2)/18)
				nn=(i-1)*K+(j)
				a=which(map==nn)
				y[a]=y[a]+mu
			}
		}
	}
y
}
###################################################
uncorrect_region1=function(y,x,k)
{
	n=length(y)
	m=ncol(x)
	mm=floor(m/k)
   if(k==1) xx=x
	else 
   {
   	q=colMeans(x)/2
		w=1/sqrt(q*(1-q))
		x=t(t(x)*w)
		xx=matrix(0,n,mm)
		for(i in 1:mm)
		{
			m1=(i-1)*k+1
			m2=i*k
			xx[,i]=rowSums(x[,m1:m2])
		}
	}	
	pv=1-pchisq((n-1)*cor(y,xx)^2,1)
}
##################################################
pc_linear_region1=function(y, x, pc, k)
{
	n=length(y)
	m=ncol(x)
	mm=floor(m/k)
    if(k==1) xx=x
	else 
    {
    q=colMeans(x)/2
	w=1/sqrt(q*(1-q))
	x=t(t(x)*w)
	xx=matrix(0,n,mm)
	for(i in 1:mm)
	{
		m1=(i-1)*k+1
		m2=i*k
		xx[,i]=rowSums(x[,m1:m2])
	}
    }
	xx=cbind(y,xx)
	x=lm(xx~pc)$resid
	pv=1-pchisq((n-1)*cor(x[,1],x[,-1])^2,1)
}
################################################
pc_nonp_region1=function(y, x, pc, k, h, kk)
{
	n=length(y)
	m=ncol(x)
	mm=floor(m/kk)
    if(kk==1) xx=x
	else 
    {
    q=colMeans(x)/2
	w=1/sqrt(q*(1-q))
	x=t(t(x)*w)
	xx=matrix(0,n,mm)
	for(i in 1:mm)
	{
		m1=(i-1)*kk+1
		m2=i*kk
		xx[,i]=rowSums(x[,m1:m2])
	}
    }
	x=cbind(y,xx)
	#h=h^(1/k)
	for(i in 1:k) pc[,i]=(pc[,i]-min(pc[,i]))/(max(pc[,i])-min(pc[,i]))
	ww=matrix(0,n,n)
	for(i in 1:n) 
	{
		w=rep(1,n)
		for(j in 1:k)
		{
			t=(pc[,j]-pc[i,j])/h
			t=(1-t^2)^2*(t^2<1)
			w=w*t
		}
	ww[i,]=w/sum(w)
	}
	x=x-ww%*%x
	pv=1-pchisq((n-1)*cor(x[,1],x[,-1])^2,1)
}
#################################################
uncorrect_region=function(y,x)
{
	n=length(y)
	#y1=rnorm(n)
	q=colMeans(x)/2
	w=1/sqrt(q*(1-q))
	x=x%*%w
	pv=1-pchisq((n-1)*cor(y,x)^2,1)
	#pv1=1-pchisq((n-1)*cor(y1,x)^2,1)
	#pvalue=c(pv,pv1)
}
###############################################
GC_lambda=function(pv,lambda)
{
     T=qchisq(1-pv,1)
     pvalue=1-pchisq(T/lambda,1)
}
##################################################
pc_linear_region=function(y, x, pc)
{
	n=length(y)
	q=colMeans(x)/2
	w=1/sqrt(q*(1-q))
	x=x%*%w
	x=cbind(y,x)
	x=lm(x~pc)$resid
	pv=1-pchisq((n-1)*cor(x[,1],x[,-1])^2,1)
}
##############################################
pc_nonp_region=function(y, x, pc, k, h)
{
	n=length(y)
	q=colMeans(x)/2
	w0=1/sqrt(q*(1-q))
	x=x%*%w0
	x=cbind(y,x)
	#h=h^(1/k)
	for(i in 1:k) pc[,i]=(pc[,i]-min(pc[,i]))/(max(pc[,i])-min(pc[,i]))
	ww=matrix(0,n,n)
	for(i in 1:n) 
	{
		w=rep(1,n)
		for(j in 1:k)
		{
			t=(pc[,j]-pc[i,j])/h
			t=(1-t^2)^2*(t^2<1)
			w=w*t
		}
	ww[i,]=w/sum(w)
	}
	x=x-ww%*%x
	pv=1-pchisq((n-1)*cor(x[,1],x[,-1])^2,1)
}
###############################################
quan_to_qual=function(y,percent)
{
  nn=floor(sample*(1-percent))
  aa=sort(y)[nn]
  y=(y>aa)
}
#############################################
uncorrect_lambda=function(y,x)
{
	n=length(y)
    T=(n-1)*cor(y,x)^2
    lambda=median(T)/0.456
    	
}
##############################################
pc_per_region=function(y, x, pc, numper)
{
	n=length(y)
	q=colMeans(x)/2
	w=1/sqrt(q*(1-q))
	x=x%*%w
	x=cbind(y,x)
	x=lm(x~pc)$resid
	pv=1-pchisq((n-1)*cor(x[,1],x[,-1])^2,1)
}
#############################################
gen_power=function(nsample, x, herid, nprot, alpha, ncau)
{
#herid=herit
#alpha=0.01
x1=x
q=colMeans(x1)/2
index=which(q<(alpha+0.5))
n0=length(index)
ii=sample(c(1:n0))
index=index[ii[1:ncau]]
nprot=floor(nprot*ncau)
xcau=x1[,index]
qcau=q[index]
nc=ncau-nprot
varr=2*qcau*(1-qcau)
u=runif(ncau)
u=u/sum(u)
	
beta=sqrt(herid*u/varr)
if(nprot>=1)
	{
	beta[(nc+1):ncau]=-beta[(nc+1):ncau]
	}
y=rowSums(t(beta*t(xcau)))#+rnorm(nsample)
#if(model<3) trait_power(sample, K, model, RR) 
}
####################################################
trait_nulla=function(sample, K, model, RR, m, kk)
{
	nk=K^2 # nk is the total number of grids
	c=sample/nk
	map=rep(0,sample)
	
	nn=0
	for (i in 1:nk)
	{
		for(j in 1:c)	{nn=nn+1;	map[nn]=i}
	}
	if(model == 0)y=rnorm(sample)
	if(model==1) #
	{
		y=rnorm(sample)
		i0=j0=6
		for(i in i0:(i0+m)) 
		{
			for(j in j0:(j0+m))
			{
				nn=(i-1)*K+(j)
				a=which(map==nn)
				y[a]=y[a]+RR
			}
		}
		i0=6;j0=K-6
		for(i in i0:(i0+m)) 
		{
			for(j in j0:(j0+m))
			{
				nn=(i-1)*K+(j)
				a=which(map==nn)
				y[a]=y[a]+RR
			}
		}
		i0=K-6;j0=6
		for(i in i0:(i0+m)) 
		{
			for(j in j0:(j0+m))
			{
				nn=(i-1)*K+(j)
				a=which(map==nn)
				y[a]=y[a]+RR
			}
		}
		i0=K-6;j0=K-6
		for(i in i0:(i0+m)) 
		{
			for(j in j0:(j0+m))
			{
				nn=(i-1)*K+(j)
				a=which(map==nn)
				y[a]=y[a]+RR
			}
		}
	}

	i0=j0=6
	if(model==2)
	{
		y=rnorm(sample)
		r=0.4*RR
		for(i in 1:K)
		{
			for(j in 1:K)
			{
				mu=r*exp(-((i-i0)^2+(j-j0)^2)/kk**2/2)
				nn=(i-1)*K+(j)
				a=which(map==nn)
				y[a]=y[a]+mu
			}
		}
	}
y
}
######################################
gen_kpop=function(sample, L, Fst, k, qq, mu)
{
	F=Fst
	if(k==1) q1=sample(qq,L,replace=TRUE)
   else 
   {
   	q=matrix(0,L,k)
   	q0=rep(0,k)
    
	nn=0
	while (nn < L)
	{		
	
      p=sample(qq,1,replace=TRUE)	
		a=p*(1-F)/F
		b=(1-p)*(1-F)/F
      bb=0
      for(i in 1:k)
      {
			q0[i]=rbeta(1,a,b)
         bb=bb+q0[i]
			#print(paste("i=",i,"q0=",q0[i]))
      }
		bb=bb/k
		if ((bb>0.0005)) 
		{
			nn=nn+1
			q[nn,]=q0
		}
	}
}
##################################
kk=k
nn=sample/kk
m=0
g=matrix(0,sample,L)
y=rep(0,sample)
for(i in 1:kk)
{
   for( j in 1:nn)
   {
      m=m+1
      y[m]=(i-1)*mu
      if(k==1) q0=q1
      else q0=q[,i]
      g[m,]=((runif(L)-q0) < 0)+((runif(L)-q0) < 0)
    }
}
y=y+rnorm(sample)
data=list(y,g)
}

#############################################
gen_kpop_r=function(sample, L, Fst, k, qq, mu, npop)
{
	F=Fst
	if(k==1) q1=sample(qq,L,replace=TRUE)
    else 
    {
    q=matrix(0,L,k)
    q0=rep(0,k)
    
	nn=0
	while (nn < L)
	{		
      p=sample(qq,1,replace=TRUE)	
		a=p*(1-F)/F
		b=(1-p)*(1-F)/F
      bb=0
      for(i in 1:k)
      {
			q0[i]=rbeta(1,a,b)
      	bb=bb+q0[i]
      }
		
		bb=bb/k
		if ((bb>0.0005)) 
		{
			nn=nn+1
			q[nn,]=q0
		}
	}
}
##################################
kk=k
nn=sample/kk
m=0
g=matrix(0,sample,L)
y=rep(0,sample)
for(i in 1:kk)
{
   for( j in 1:nn)
   {
      m=m+1
      if(i>(k-npop))y[m]=mu
      if(k==1) q0=q1
      else q0=q[,i]
      g[m,]=((runif(L)-q0) < 0)+((runif(L)-q0) < 0)
    }
}
y=y+rnorm(sample)
data=list(y,g)
}
##############################################