!$Author$
!$Date$
!$Revision$
!$HeadURL$

! Functions for signed log-likelihood ratio statistic from:
! Tian, L. and J. Wu. 2006. Confidence Intervals for the Mean of 
! Lognormal Data with Excess Zeros. Biometrical Journal 48:149-156

subroutine r0_z0(n,x,fvec,ni,iv,nr,rv,iflag)
    integer n,ni,iv(ni),nr
    integer, intent(inout) :: iflag
    double precision x(n),fvec(n),rv(nr)
    double precision :: r0

    ! variables passed to r0
    ! x(1) = psi
    ! iv(1) = nzero
    ! iv(2) = ngtz
    ! rv(1) = psihat
    ! rv(2) = muhat
    ! rv(3) = etahat
    ! rv(4) = t_sum
    ! rv(5) = s_sum

    fvec(1) = r0( x(1), iv(1), iv(2), rv(1), rv(2), rv(3), rv(4), rv(5), iflag ) - rv(6)

end subroutine

double precision function r0(psi, nzero, ngtz, psihat, muhat, etahat, t_sum, s_sum, iflag)
    double precision, intent(in) :: psi, psihat, muhat, etahat, t_sum, s_sum
    integer, intent(in) :: nzero, ngtz
    integer, intent(inout) :: iflag
    double precision :: muhatpsi, etahatpsi

    double precision :: log_likelihood, r0
    double precision tmp0, B, one

    double precision :: constraint, muint, etaint

    ! variables required to use hybrd.f
    integer :: info, nfev, iv(2)
    double precision :: x(2), fvec(2), rv(3), diag(2), fjac(2,2), r(3)
    double precision :: qtf(2),wa1(2),wa2(2),wa3(2),wa4(2)

    external fnu

    ! check solution constraint to get intial guess in proper region
    ! initial guess
    muhatpsi  = muhat
    etahatpsi = etahat
    constraint = psi - muhatpsi - etahatpsi/2
    !write(*,*) "psi, muhatpsi, etahatpsi, constraint", psi, muhatpsi, etahatpsi, constraint
    ! find the point of intersection of the perpendicular to constraint boundary
    muint = 0.4*(2*psi + 0.5*muhatpsi - etahatpsi)
    etaint = 0.4*(psi - muhatpsi + 2*etahatpsi)
    !write(*,*) "intersection point - muint, etaint", muint, etaint
    ! check value against constraint boundary
    if( constraint .ge. 0.0 ) then
       ! point is in wrong region, out of bounds, refect point into solution region
       muhatpsi = 2*muint - muhatpsi
       etahatpsi = 2*etaint - etahatpsi
       constraint = psi - muhatpsi - etahatpsi/2
       !write(*,*) "nfev, psi, muhatpsi, etahatpsi, constraint", 0, psi, muhatpsi, etahatpsi, constraint
    end if
    nfev = 0
    do while( nfev .lt. 10 )
       ! evaluate functions
       x(1) =  muhatpsi
       x(2) =  etahatpsi
       iv(1) = nzero
       iv(2) = ngtz
       rv(1) = psi 
       rv(2) = s_sum
       rv(3) = t_sum
       !write(*,*) "nfev, psi, muhatpsi, etahatpsi, constraint", nfev, psi, muhatpsi, etahatpsi, constraint
       call fnu(2,x,fvec,2,iv,3,rv,info)
       nfev = nfev + 1
       if( (fvec(1) .lt. 0.0d0) .or. (fvec(2) .lt. 0.0d0) ) then
          muhatpsi = 0.5*(muint + muhatpsi)
          etahatpsi = 0.5*(etaint + etahatpsi)
          constraint = psi - muhatpsi - etahatpsi/2
          !write(*,*) "nfev, psi, muhatpsi, etahatpsi, constraint", nfev, psi, muhatpsi, etahatpsi, constraint
       else
          exit
       end if
    end do 
    
       
    ! for value of psi, find muhatpsi, etahatpsi
    x(1) =  muhatpsi
    x(2) =  etahatpsi
    iv(1) = nzero
    iv(2) = ngtz
    rv(1) = psi 
    rv(2) = s_sum
    rv(3) = t_sum
    ! call for signed log-likilihood ratio
    call hybrd(fnu,2,x,fvec,2,iv,3,rv,1.0d-6,100,1,1, &
               1.0d-9,diag,1,1.0d2,0,info,nfev,fjac, &
               2,r,3,qtf,wa1,wa2,wa3,wa4)
    muhatpsi  = x(1) 
    etahatpsi = x(2) 

    if( info .ne. 1 ) then
       iflag = -info
       !write(*,*) "info, muhatpsi, etahatpsi ", info, muhatpsi, etahatpsi
       write(*,*) "# ERROR, confidence interval calculation did not converge, ", nzero+ngtz, " values."
       r0 = 0.0d0
       return
    end if

    tmp0 = 0.0d0
    constraint = psi - muhatpsi - etahatpsi/2
    if( (constraint .ge. 0.0d0) .or. (etahatpsi .le. 0.0d0) ) then
       iflag = -1
       !write(*,*) "muhat ", muhat
       !write(*,*) "etahat ", etahat
       !write(*,*) "nzero ", nzero
       !write(*,*) "ngtz ", ngtz
       !write(*,*) "psi ", psi
       !write(*,*) "s_sum ", s_sum
       !write(*,*) "t_sum ", t_sum
       !write(*,*) "muhatpsi ", muhatpsi
       !write(*,*) "etahatpsi ", etahatpsi
       !write(*,*) "constraint", constraint
       write(*,*) "# ERROR, confidence interval calculation, bad likelihood, ", nzero+ngtz, " values."
    else
       ! signed log-likelihood ratio
       B = psihat - psi
       one = 1.0
       IF (B .EQ. 0.0) B = +0.0
       tmp0 = sign(one, B) &
          * sqrt( 2.0 &
          * (log_likelihood(nzero, ngtz, psihat, muhat, etahat, t_sum, s_sum) &
          - log_likelihood(nzero, ngtz, psi, muhatpsi, etahatpsi, t_sum, s_sum)))
    end if
    r0 = tmp0
end function

subroutine fnu(n,x,fvec,ni,iv,nr,rv,iflag)
    integer, intent(in) :: n,ni,iv(ni),nr
    integer, intent(inout) :: iflag
    double precision, intent(in) :: rv(nr)
    double precision, intent(inout) :: x(n)
    double precision, intent(out) :: fvec(n)
    double precision B, one

    integer :: nzero, ngtz
    double precision :: constraint, muhatpsi, etahatpsi, psi, s_sum, t_sum 

    ! variables passed to fnu
    muhatpsi  = x(1) 
    etahatpsi = x(2) 
    nzero  = iv(1)
    ngtz   = iv(2)
    psi    = rv(1)
    s_sum  = rv(2)
    t_sum  = rv(3)
                   
    constraint = psi - muhatpsi - etahatpsi/2
    ! check for out of range value
    if( constraint > 700.0 ) then
        ! return flag indicating failure of routine, please terminate
        iflag = -1
    else
        !write(*,*) 'fnu: constraint: ', constraint
        fvec(1) = (nzero * exp(constraint)) &
                / (1 - exp(constraint)) &
                - ngtz &
                + s_sum / etahatpsi &
                - ngtz * muhatpsi / etahatpsi
        fvec(2) = (nzero * exp(constraint)) &
                / 2 / (1 - exp(constraint)) &
                - ngtz / 2.0 &
                - ngtz / 2.0 / etahatpsi &
                + t_sum / 2.0 / etahatpsi / etahatpsi &
                - muhatpsi * s_sum / etahatpsi / etahatpsi &
                + ngtz * muhatpsi * muhatpsi / 2.0 / etahatpsi / etahatpsi

        ! check constraint
        !constraint = psi - muhatpsi - etahatpsi/2
        if( constraint .gt. 0 ) then
          one = 1.0
          B = fvec(1)
          IF (B .EQ. 0.0) B = +0.0
           fvec(1) = fvec(1) + sign(one,B)*constraint
          B = fvec(2)
          IF (B .EQ. 0.0) B = +0.0
           fvec(2) = fvec(2) + sign(one,B)*constraint
        end if
    end if

end subroutine

double precision function log_likelihood(nzero, ngtz, psi, mu, eta, t_sum, s_sum)
    integer, intent(in) :: nzero, ngtz
    double precision, intent(in) :: psi, mu, eta, t_sum, s_sum

    real pi
    parameter (pi = 3.1415926535)

    log_likelihood = nzero * log(1 - exp(psi - mu - eta/2)) &
                   + ngtz * (psi - mu - eta/2) &
                   - ngtz * log(2*pi) / 2 &
                   - ngtz * log(eta) / 2 &
                   - t_sum / 2 / eta &
                   + mu * s_sum / eta &
                   - ngtz * mu * mu / 2 / eta
end function