##-----------------------------------------------------------------------------------------------------------------------------##
##----This program is for the paper "Hayter, Lin & Liu (2014), Establishing practical equivalence between three treatments"----##
##--------For the integral in formula (2), the outer integral is calculated by the trapezoidal method in interval [0,U]--------##
##------------------the inner integral is approximated by the trapezoidal method in interval [l(s),u(s)]-----------------------##
##-----------------------------------------------------------------------------------------------------------------------------##

g_newR<-function(gs,M,hstep,S,Nout,c,delta,n,d)
{
  gs[1]<-0
  for (it in 2:Nout) {
    L<- -c*S[it]/sqrt(2)+(delta[2]-delta[1])/(sqrt(2)*d[1])
    U<-  c*S[it]/sqrt(2)+(delta[2]-delta[1])/(sqrt(2)*d[1])
    
    if (L==U) {
      gs[it]<-0
      next
      }
    if ((U-L)/hstep==0) {
      N<-floor((U-L)/hstep)+1
      } else {
        N<-floor((U-L)/hstep)+2
        }

    x<-rep(0,N)
    Step<-(U-L)/(N-1)
    for (i in 1:N) {
      x[i]<-L+(i-1)*Step
      }

    if (M==1) {   #trapezoidal method
      a<- -c*S[it]*d[2]+delta[3]-delta[1]-x[1]/(sqrt(2)*n[1]*d[1])
      b<- -c*S[it]*d[3]+delta[3]-delta[2]+x[1]/(sqrt(2)*n[2]*d[1])
      l<-max(a,b)/sqrt(1/(n[1]+n[2])+1./n[3]);
      a<-c*S[it]*d[2]+delta[3]-delta[1]-x[1]/(sqrt(2)*n[1]*d[1])
      b<-c*S[it]*d[3]+delta[3]-delta[2]+x[1]/(sqrt(2)*n[2]*d[1])
      u<-min(a,b)/sqrt(1/(n[1]+n[2])+1./n[3])
      temp<-dnorm(x[1])*(pnorm(u)-pnorm(l))
      Sum<-temp
      
      for (i in 2:(N-1)) {
        a<- -c*S[it]*d[2]+delta[3]-delta[1]-x[i]/(sqrt(2)*n[1]*d[1])
        b<- -c*S[it]*d[3]+delta[3]-delta[2]+x[i]/(sqrt(2)*n[2]*d[1])
        l<-max(a,b)/sqrt(1/(n[1]+n[2])+1./n[3])
        a<-c*S[it]*d[2]+delta[3]-delta[1]-x[i]/(sqrt(2)*n[1]*d[1])
        b<-c*S[it]*d[3]+delta[3]-delta[2]+x[i]/(sqrt(2)*n[2]*d[1])
        u<-min(a,b)/sqrt(1/(n[1]+n[2])+1./n[3])
        temp<-dnorm(x[i])*(pnorm(u)-pnorm(l))
        Sum<-Sum+2.0*temp
		    }
		    
      a<- -c*S[it]*d[2]+delta[3]-delta[1]-x[N]/(sqrt(2)*n[1]*d[1])
      b<- -c*S[it]*d[3]+delta[3]-delta[2]+x[N]/(sqrt(2)*n[2]*d[1])
      l<-max(a,b)/sqrt(1/(n[1]+n[2])+1./n[3]);
      a<-c*S[it]*d[2]+delta[3]-delta[1]-x[N]/(sqrt(2)*n[1]*d[1])
      b<-c*S[it]*d[3]+delta[3]-delta[2]+x[N]/(sqrt(2)*n[2]*d[1])
      u<-min(a,b)/sqrt(1/(n[1]+n[2])+1./n[3])
      temp<-dnorm(x[N])*(pnorm(u)-pnorm(l))
      Sum<-Sum+temp
      Sum<-Sum*Step/2
      gs[it]<-Sum
      } else {
        print("Method for inner integral error!\n")
        }
    }
    gs
}

sqrt_chisq<-function(x,v)
{
  if (x<=0) return(0)
#  temp<-v^(v/2)*x^(v-1)*exp(-v*x^2/2)/2^(v/2-1)/exp(lgamma(v/2))
  temp<-v/2*log(v)+(v-1)*log(x)-v*x^2/2.-(v/2.-1)*log(2)-log(gamma(v/2))
  temp<-exp(temp)
  temp
}

prob_R<-function(M,U,N,c,delta,n,v,hstep,d)
{
  S<-rep(0,N)
  Step<-U/(N-1)
  for (i in 1:N)
    S[i]<-(i-1)*Step

  if (M==1) {  #trapezoidal method
    gs<-rep(0,N)
    gs<-g_newR(gs,M,hstep,S,N,c,delta,n,d)

    Sum<-0
    for (i in 2:(N-1))
      Sum<-Sum+2*gs[i]*sqrt_chisq(S[i],v)
    Sum<-Sum+gs[N]*sqrt_chisq(S[N],v)
    Sum<-Sum*Step/2
    } else {
      print("Method for outer integral error!\n")
      }

  Sum
}


############# Starting the main routine ##################
#Example 1
n<-c(6,12,21) 
X<-c(9.9,10.5,10.9)
S<-1   
v<-0
for (i in 1:3)
  v<-v+n[i]
v<-v-3
d<-rep(0,3)
d[1]<-1./sqrt(2)*sqrt(1./n[1]+1./n[2])
d[2]<-1./sqrt(2)*sqrt(1./n[1]+1./n[3])
d[3]<-1./sqrt(2)*sqrt(1./n[2]+1./n[3])
R<-max(abs(X[1]-X[2])/(S*d[1]),abs(X[1]-X[3])/(S*d[2]),abs(X[2]-X[3])/(S*d[3]))

alpha<-0.05
M<-1
U<-8
Nout<-200
L<-5
hstep<-2*L/(200-1)

First<-1
delta<-rep(0,3)
last_c<-rep(0,3)
delta_lfc<-rep(0,3)
temp<-NA

Delta<-seq(from=0,to=4,by=0.01)
NDelta<-length(Delta)
c_result<-rep(0,NDelta)
for (it in 1:NDelta) {
  old_c<-100
  for (k in 1:3) {
    if (k==1) {
      delta[2]<-0
      delta[3]<-Delta[it]
      } else if (k==2) {
        delta[1]<-0
        delta[3]<-Delta[it]
        } else {
          delta[1]<-0
          delta[2]<-Delta[it]
          }

    if (First!=1) {
      c<-last_c[k]
      } else {
        c<-0.42 #Here c is the initial value and some value shoule be given with the prior information
        }

    a<-seq(from=0,to=Delta[it],by=0.01)
    Na<-length(a)
    success<-0
    while (!success) {
      for (i in 1:Na) {
        delta[k]<-a[i]
        temp<-prob_R(M,U,Nout,c,delta,n,v,hstep,d)-alpha
        if (temp>0) break
        }
      if (temp>0) {
        success<-1
        } else {
          c<-c+0.001
          }
      }

    lfc<-1
    if (c-0.001<old_c) {
      for (i in 1:3) delta_lfc[i]<-delta[i]
      } else if (c-0.001==old_c) {
        lfc<-lfc+k-1
        }

    old_c<-min(old_c,c-0.001)
    last_c[k]<-c-0.001
    }
  First<-First+1      
  c_result[it]=old_c  
  cat("When Delta=",Delta[it],",   c=",c_result[it],"\n") 
  }