ftn_catt=function(data, score){
	nn=apply(data,2,sum)
	n=sum(nn)

	Rbar=sum(nn*score)/n
	s2=sum(nn*(score-Rbar)^2)
	phi=sum(data[1,])/n

	catt=sum(data[1,]*(score-Rbar))/sqrt(phi*(1-phi)*s2)
	p.value = 2*(1-pnorm(abs(catt)))


	out=list(estimate = catt,p.value = p.value)
	return(out) 

} # end ftn_catt

ftn_pmin2 = function(y){
	f1_inv = qchisq(y,lower.tail=F,1)
	integrand = function(z){exp(-z/2)*asin(2*f1_inv/z-1)}

	pmin2 = 0.5*exp(-f1_inv/2) + 0.5*y - 1/(2*pi)*integrate(integrand,f1_inv,-2*log(y))$value
	return(pmin2)
} # end ftn_pmin2


ftn_CLRT = function(r0,r1,r2,s0,s1,s2){

	r = r0+r1+r2
	s = s0+s1+s2

	n0 = r0+s0
	n1 = r1+s1
	n2 = r2+s2

	n = n0+n1+n2

	f0 = r0/n0
	f1 = r1/n1
	f2 = r2/n2

	lambda1 = f1/f0
	lambda2 = f2/f0


	ftn_l2 = function(f0,f1,f2){
	r0*log(f0) + s0*log(1-f0)+r1*log(f1) + s1*log(1-f1)+r2*log(f2) + s2*log(1-f2)
	}

	f0_null = r/n
	l0 = ftn_l2(f0_null,f0_null,f0_null)


	if(lambda1 >= 1 & lambda2 >= lambda1 & lambda2 > 1){
	CLRT = 2*(ftn_l2(f0,f1,f2) - l0)
	} else if(lambda2 >0 & lambda1 >= lambda2 & lambda2 < 1){
	CLRT = 2*(ftn_l2(f0,f1,f2) - l0)
	} else{
	f0_r = (r0+r1)/(n0+n1)
	f1_r = (r0+r1)/(n0+n1)
	f2_r = r2/n2

	f0_d = r0/n0
	f1_d = (r1+r2)/(n1+n2)
	f2_d = (r1+r2)/(n1+n2)

	CLRT_R = 2*(ftn_l2(f0_r,f1_r,f2_r) -l0)
	CLRT_D = 2*(ftn_l2(f0_d,f1_d,f2_d) -l0)

	CLRT = max(CLRT_R,CLRT_D)
	}

	return(CLRT)

}  # end ftn_CLRT

ftn_pCLRT = function(BB,r0,r1,r2,s0,s1,s2,CLRT){

	r = r0+r1+r2
	s = s0+s1+s2

	n0 = r0+s0
	n1 = r1+s1
	n2 = r2+s2

	n = n0+n1+n2

	count = 0
	for(bb in 1:BB){
	b_r = rmultinom(1,r,c(n0/n,n1/n,n2/n))
	b_s = rmultinom(1,s,c(n0/n,n1/n,n2/n))

	b_r0 = b_r[1];	b_r1 = b_r[2];	b_r2 = b_r[3]
	b_s0 = b_s[1];	b_s1 = b_s[2];	b_s2 = b_s[3]
	
	b_CLRT = ftn_CLRT(b_r0,b_r1,b_r2,b_s0,b_s1,b_s2)

	count = count + as.numeric(b_CLRT > CLRT)
	} # end bb

	return(count/BB)
} # end ftn_pCLRT


ftn_MAX = function(r0,r1,r2,s0,s1,s2){
	xx = rbind(c(r0,r1,r2),c(s0,s1,s2))

#	if(min(xx) == 0){xx = xx+0.5}

	rec = ftn_catt(xx,c(0,0,1))$estimate
	add = ftn_catt(xx,c(0,0.5,1))$estimate
	dom = ftn_catt(xx,c(0,1,1))$estimate
	
	MAX = max(abs(rec),abs(add),abs(dom))
	return(MAX)
}

ftn_p_MAX_CLRT = function(BB,r0,r1,r2,s0,s1,s2,MAX,CLRT){

	r = r0+r1+r2
	s = s0+s1+s2

	n0 = r0+s0
	n1 = r1+s1
	n2 = r2+s2

	n = n0+n1+n2

	count_CLRT = 0
	count_MAX = 0

	for(bb in 1:BB){

#	if(bb%%100000 == 0){cat("bb=",bb/BB*100,"%\n")}
	b_r = rmultinom(1,r,c(n0/n,n1/n,n2/n))
	b_s = rmultinom(1,s,c(n0/n,n1/n,n2/n))

	b_r0 = b_r[1];	b_r1 = b_r[2];	b_r2 = b_r[3]
	b_s0 = b_s[1];	b_s1 = b_s[2];	b_s2 = b_s[3]
	
	b_CLRT = ftn_CLRT(b_r0,b_r1,b_r2,b_s0,b_s1,b_s2)
	count_CLRT = count_CLRT + as.numeric(b_CLRT > CLRT)

	b_MAX = ftn_MAX(b_r0,b_r1,b_r2,b_s0,b_s1,b_s2)
	count_MAX = count_MAX + as.numeric(b_MAX > MAX)

	} # end bb
	
	out = list(pMAX = count_MAX/BB,pCLRT = count_CLRT/BB)
	return(out)

} # end ftn_p_MAX_CLRT



ftn_p_MAX= function(BB,r0,r1,r2,s0,s1,s2,MAX){

	r = r0+r1+r2
	s = s0+s1+s2

	n0 = r0+s0
	n1 = r1+s1
	n2 = r2+s2

	n = n0+n1+n2

	count_MAX = 0

	for(bb in 1:BB){

	if(bb%%100000 == 0){cat("bb=",bb/BB*100,"%\n")}
	b_r = rmultinom(1,r,c(n0/n,n1/n,n2/n))
	b_s = rmultinom(1,s,c(n0/n,n1/n,n2/n))

	b_r0 = b_r[1];	b_r1 = b_r[2];	b_r2 = b_r[3]
	b_s0 = b_s[1];	b_s1 = b_s[2];	b_s2 = b_s[3]
	
	b_MAX = ftn_MAX(b_r0,b_r1,b_r2,b_s0,b_s1,b_s2)
	count_MAX = count_MAX + as.numeric(b_MAX > MAX)

	} # end bb
	
	return(count_MAX/BB)

} # end ftn_p_MAX