###########################################
### Computes exact confidence Interval  ###
### for the difference of              ####
### two proportions based on           ###
### inverting one two-sided score test  ###
### of Mee or MN.                     ###
#####################################
score.z <- function(y,n,delta.null,alternative = "two.sided") {
  N=sum(n)
  C=sum(y)
  L3=N
  L2=(n[1]+2*n[2])*delta.null-N-C
  L1=(n[2]*delta.null-N-2*y[2])*delta.null+C
  L0=y[2]*delta.null*(1-delta.null)
  c=L2^3/(3*L3)^3 - L1*L2/(6*L3^2) + L0/(2*L3)
  b=ifelse(c>=0,1,-1)*sqrt(L2^2/(3*L3)^2 - L1/(3*L3))
  d <- min(c/b^3,1)
  d <- max(d, -1)
  a=(3.14159265358979+acos(d))/3
  p2=2*b*cos(a)-L2/(3*L3)
  p2 = max(p2, 0)
  p2 = min(p2,1)
  p1=p2+delta.null
  p1 = max(p1, 0)
  p1 = min(p1,1)
  var0 <- p1*(1 - p1)/n[1] + p2*(1 - p2)/n[2]
  z2=(y[1]/n[1]-y[2]/n[2]-delta.null)^2/var0 #Mee or MN doesn't matter in exact test
  if (var0 <= 0) z2 = 0
  return(round(z2,7))
}

exactP.pi2 <- function(pi2, delta0, S.MN, S.obs, midP, n, ygrid) {
 pitrue=apply(ygrid,1,function(r) exp(sum(dbinom(r,size=n,prob=c(delta0+pi2,pi2),log=TRUE))))
 if(!midP) return(sum(pitrue*(S.MN>=S.obs)))
 else return(sum(pitrue*(0.5*(S.MN>=S.obs)+0.5*(S.MN>S.obs))))
  #be very careful about numerical accuracy: Eg: for n=c(6,6), y=c(5,2) and y=c(1,4) need to give same vale of test statistic, but not the case in R. One is considered larger.
}

exactP <- function(delta0, midP=FALSE, y, n, ygrid) {
  S.MN <- apply(ygrid, 1, score.z, n=n, delta.null=delta0)
  S.obs <- score.z(y,n,delta.null=delta0)
  if (abs(delta0) >= 0.99999) return(0)  #trick to get optimize to give upper/lower bound of +/- 1 if L-shaped
  maxP <- optimize(exactP.pi2, lower=max(0,-delta0), upper=min(1,1-delta0), delta0=delta0, S.MN=S.MN,
                               S.obs=S.obs, midP=midP, n=n, ygrid=ygrid, maximum=TRUE)$objective
  return(round(maxP-0.05,7))
}

exactCI <- function(y, n, conf.level=0.95, midP=FALSE) {
 ygrid <- expand.grid(0:n[1],0:n[2])
 p=c((y[1]+1)/(n[1]+2),(y[2]+1)/(n[2]+2))
 point.est= p[1]-p[2]
 c=qnorm(1-(1-conf.level)/2)
 CI.ini = point.est + cbind(-1,1)*c*sqrt(sum(p*(1-p)/n))
 LB <- uniroot(exactP, interval=c(max(-0.99998,CI.ini[1]),point.est), midP=midP, y=y, n=n, ygrid=ygrid, extendInt="up")$root
 UB <- uniroot(exactP, interval=c(point.est,min(0.99998,CI.ini[2])), midP=midP, y=y, n=n, ygrid=ygrid, extendInt="down")$root
 return(round(c(LB,UB),4))
}

n=c(6,6)
y=c(5,1)
exactCI(y,n)
exactCI(y,n,midP=TRUE)