approx.pvalue<-function(x,t,n,delta,B=10){
  s<-score(n,theta=abs(delta),alt=-1)
  list(P.approx=TOST(s)$P[(t^2+t+2)/2+x])
}

# Just type "approx.pvalue(x,t,n,delta)"


etm.pvalue<-function(x,t,n,delta,B=10){
  j<-(t^2+t+2)/2+x
  s<-score(n,theta=abs(delta),alt=-1)
  es<-ESTEP(s,theta=abs(delta))
  ets<-TOST(es)
  etms<-MSTEP(ets,theta=abs(delta),J=j,B=B)
  list(mle=(2*x-t)/n,T=s$T[j],P.approx=TOST(s)$P[j],etmP=etms$P)
}

# Just type "etm.pvalue(x,t,n,delta)"


ESTEP<-function(stat,theta=0,J=NULL,dec.places=7,prin=F){
# FUNCTION CALCULATES ESTIMATED P-VALUES FROM PROVIDED INITIAL P-VALUES
#
# ARGUMENTS:
#         stat - list with components x, t and P
#        theta - null (boundary) value of theta
#            J - index of datasets to consider (default is all)
#   dec.places - default 5
#
# VALUE
# List with components x, t - data
#                     P.old - initial
#                      phat - estimated nuisance parameter
#                     P - estimated P-values (mis-named but to be 
#                             consistent with functions MSTEP)
#
  P<-NULL
  n<-max(stat$t)
  if(missing(J)){J=1:((n+1)*(n+2)/2)}  
  phat<-profile.mle(stat$x,stat$t,n,theta,tol=0)
  obj<-function(phi,theta,xvals,tvals,n){
    if(phi>0){eta<-0.5*(theta+phi)/phi}
    if(phi==0){eta<-0.5}
    if(eta>1){eta<-1}
    sum(dbinom(xvals,tvals,eta)*dbinom(tvals,n,phi))
  }
  for(j in J){
    if(prin){print(j)}
#    if(round(j/100)*100==j){print(j)}
    take<-(stat$P<=stat$P[j])
    xvals<-stat$x[take]
    tvals<-stat$t[take]
    P<-c(P,obj(phat[j],theta,xvals,tvals,n))
  }
  list(x=stat$x[J],t=stat$t[J],P.old=stat$P[J],phat=phat[J],
     P=round(10^dec.places*P)/10^dec.places)
}

MSTEP<-function(stat,theta=0,J=null,B=10,dec.places=7,prin=F,plt=F){
# FUNCTION CALCULATES MAXIMISED P-VALUES FROM PROVIDED INITIAL P-VALUES
#
# ARGUMENTS:
#         stat - list with components x, t and P
#        theta - null value of theta (default zero)
#            J - index of datasets to consider (default is all)
#            B - search effort on nuisance parameter
#   dec.places - default 5
#       phimax - maximum value of nuisance parameter
#
# VALUE
# List with components x, t - data
#                     P.old - initial
#                         P - maximised P-values 
#
  grid<-function(B,theta=.5){
      b<-round(B/3)+1
      x<-(0:b)/b
      y<-.5+.5*x^theta
      grid.vals<-round(10000*c(1-y,x,y))/10000
      sort(unique(grid.vals))
    }
# 
  P<-NULL
  obj<-function(phi,theta,xvals,tvals,n){
    if(phi>0){eta<-0.5*(theta+phi)/phi}
    if(phi==0){eta<-0.5}
    if(eta>1){eta<-1}
    sum(dbinom(xvals,tvals,eta)*dbinom(tvals,n,phi))
  }
  n<-max(stat$t)
  if(missing(J)){J=1:((n+1)*(n+2)/2)}  
  PHI<-abs(theta)+(1-abs(theta))*grid(B,theta=1)
#
  for(j in J){
    if(prin){print(j)}
    if(round(j/50)*50==j){print(j)}
    xvals<-stat$x[stat$P<=stat$P[j]]
    tvals<-stat$t[stat$P<=stat$P[j]]
    PVALS<-PMAX<-pmax<-NULL
    if(plt){pmax<-obj(mean(PHI[1]),theta,xvals,tvals,n)}
    for(i in 2:length(PHI)){
      if(plt){pmax<-c(pmax,obj(mean(PHI[i]),theta,xvals,tvals,n))}
      PVALS<-c(PVALS,obj(mean(PHI[(i-1):i]),theta,xvals,tvals,n))
      PMAX<-c(PMAX,optimise(obj,PHI[(i-1):i],maximum=TRUE,theta=theta,xvals=xvals,tvals=tvals,n=n)$objective)   
    }
     P<-c(P,max(PVALS,PMAX,obj(PHI[1],theta,xvals,tvals,n),obj(PHI[length(PHI)],theta,xvals,tvals,n)))
  }
  list(x=stat$x[J],t=stat$t[J],P.old=stat$P[J],P=round(10^dec.places*P)/10^dec.places,profile=list(x=PHI,y=pmax))
}

score<-function(n,theta=0,alt=0,tol=0.00000001,J=NULL){
# All possible values of Tango/Nam statistic
  x<-t<-NULL
  for(k in 0:n){
    x<-c(x,0:k)
    t<-c(t,rep(k,k+1))
  }
#
  if(missing(J)){J<-(1:((n+1)*(n+2)/2))}
  x<-x[J]
  t<-t[J]
  phihat<-t/n
  thhat<-(2*x-t)/n
  pmle<-profile.mle(x,t,n,theta,tol=0)
  U<-sqrt(n)*(thhat-theta)/sqrt(pmle-theta^2+tol)
  if(alt==(+1)){P<-1-pnorm(U)}
  if(alt==(-1)){P<-pnorm(U)}
  if(alt!=1&alt!=(-1)){P<-2*(1-pnorm(abs(M)))}
  list(x=x,t=t,T=U,P=P,pmle=pmle)
}


profile.mle<-function(x,t,n,theta,tol=0.0000001){
  thhat<-(2*x-t+tol)/(n+2*tol)
  phihat<-(t+tol)/(n+2*tol)
  b<-(-1)*theta*thhat-phihat
  c<-thhat*theta-(1-phihat)*theta^2
  value<-(-b/2+sqrt(b^2/4-c))
  value[value<theta]<-theta
  value[value>1]<-1
  value
}

swap<-function(stat){
#   Function converts full collection of results for testing
#   theta<delta into set of results for testing theta>-delta.
#   This avoid recomputing maximisations for the greater than
#   hypothesis since the results follow from symmetry.
#
    n<-max(stat$x)
    x<-stat$x
    t<-stat$t
    P<-stat$P
    P[P>1]<-1
    T<-stat$T
    if(is.null(T)){T<-qnorm(P)}
    for(k in 0:n){
      l<-(k^2+k+2)/2
      u<-(k+2)*(k+1)/2
       T[l:u]<-(-T[u:l])
       P[l:u]<-P[u:l]
    }
    list(x=x,t=t,T=T,P=P)
}

TOST<-function(stat){
#   Takes a test statistics for testing theta<delta and converts to 
#   the TOST test based on the maximum P-value from this and the
#   swapped test.
    tmp<-apply(cbind(swap(stat)$P,stat$P),1,max)
    list(x=stat$x,t=stat$t,P=tmp)
}

grid<-function(B,theta=.5){
      b_round(B/3)+1
      x_(0:b)/b
      y_.5+.5*x^theta
      grid.vals_round(10000*c(1-y,x,y))/10000
      sort(unique(grid.vals))
    }