anchor.diff:      ! recompute fxs
       fxs = exp(lfwa(xs))
   
  -    ! recompute integration range
  -    a1 = alpha2
  -    a2 = alpha1
  -    do i = 1,nx
  -       x = log(xs(i))
  -       f = fxs(i)
  -       if (x > alpha1 .and. x < a1 .and. f > small_number) a1 = log(xs(i-1))
  -       if (x > a2 .and. x < alpha2 .and. f > small_number) a2 = log(xs(i+1))
  -    end do
  -    
  -    alpha1 = a1
  -    alpha2 = a2 
  +    alpha1 = log(minval(xs, fxs > small_number))
  +    alpha2 = log(maxval(xs, fxs > small_number))
  +
  +!    ! recompute integration range
  +!    a1 = alpha2
  +!    a2 = alpha1
  +!    do i = 1,nx
  +!       x = log(xs(i))
  +!       f = fxs(i)
  +!       if (x > alpha1 .and. x < a1 .and. f > small_number) a1 = log(xs(i-1))
  +!       if (x > a2 .and. x < alpha2 .and. f > small_number) a2 = log(xs(i+1))
  +!    end do
  +!   
  +!    alpha1 = a1
  +!    alpha2 = a2 
   
       ! recompute the maximum
       dx = (alpha2 - alpha1) / (nx * 1.0_real64)
  
positive.code:   end function log_fpxa
  
    real(real64) function hdir(alpha)
      ! posterior mean entropy given the data and alpha
      ! computed from histogram multiplicities 
      use gamma_funcs, only: digamma
      use constants
      
  
negative.code:     alpha1 = -20.0_real64
      alpha2 = 10.0_real64
  
      ! initialize amax
      amax = 1.0_real64
      log_fpxa_amax = log_fpxa(amax) ! log p(n|a_max)
      fpa_amax = fpa(amax) ! propto p(a)
  
      ! set intervals equally spaced on log scale
      dx = (alpha2 - alpha1) / real(nx -1.0_real64)
      do i = 1,nx
         xs(i) = alpha1 + (i - 1) * dx
      end do
      xs = exp(xs)
      
      fxs = lfwa(xs)
      ! find amax such that the weight in the integral - fwa(alpha) - is maximal
      i = maxloc(fxs,1,fxs < largest)
      amax = xs(i)
      log_fpxa_amax = log_fpxa(amax)
  