module br_mod_dist implicit none ! ********************************************************************* ! * br_mod -- An F90 module for computing the Blackwell-Rao * ! * estimator given signal samples from the posterior * ! * * ! * Written by Hans Kristian Eriksen * ! * * ! * Copyright 2006, all rights reserved * ! * * ! * * ! * NB! The code is provided as is, and *no* guarantees are given * ! * as far as either accuracy or correctness goes. * ! * * ! * If used for published results, please cite these papers: * ! * * ! * - Eriksen et al. 2006, ApJ, submitted, astro-ph/0606088 * ! * - Chu et al. 2005, Phys. Rev. D, 71, 103002 * ! * * ! ********************************************************************* integer, parameter, private :: i4b = selected_int_kind(8) integer, parameter, private :: sp = selected_real_kind(5,30) integer, parameter, private :: dp = selected_real_kind(12,200) integer, parameter, private :: lgt = kind(.true.) integer(i4b), private :: lmin, lmax, numsamples, numchain logical(lgt), private :: first_eval real(dp), private :: offset real(dp), allocatable, dimension(:,:,:), private :: sigmas contains ! Initialization routines subroutine initialize_br_mod(lmin_in, sigmas_in) implicit none integer(i4b), intent(in) :: lmin_in real(sp), dimension(lmin_in:,1:,1:), intent(in) :: sigmas_in lmin = lmin_in lmax = size(sigmas_in(:,1,1)) + lmin - 1 numchain = size(sigmas_in(lmin,:,1)) numsamples = size(sigmas_in(lmin,1,:)) first_eval = .true. allocate(sigmas(lmin:lmax,numchain,numsamples)) sigmas = real(sigmas_in,dp) end subroutine initialize_br_mod subroutine clean_up_br_mod implicit none if (allocated(sigmas)) deallocate(sigmas) end subroutine clean_up_br_mod ! Base computation routine subroutine compute_br_estimator(cls, lnL) implicit none real(dp), dimension(lmin:), intent(in) :: cls real(dp), intent(out) :: lnL integer(i4b) :: i, j, l real(dp) :: subtotal, x if (first_eval) then call compute_largest_term(cls) first_eval = .false. end if ! Compute the Blackwell-Rao estimator lnL = 0.d0 do i = 1, numchain do j = 1, numsamples subtotal = 0.d0 do l = lmin, lmax x = sigmas(l,i,j)/cls(l) subtotal = subtotal + & & 0.5d0 * real(2*l+1,dp) * (-x + log(x)) - log(real(sigmas(l,i,j),dp)) end do lnL = lnL + exp(subtotal-offset) end do end do if (lnL > 1e-20) then lnL = log(lnL) else lnL = log(1e-30) end if end subroutine compute_br_estimator ! Routine for reading the Gibbs sigma samples subroutine read_gibbs_chain(filename, unit, lmax, numchains, numsamples, data) implicit none character(len=*), intent(in) :: filename integer(i4b), intent(in) :: unit integer(i4b), intent(out) :: lmax, numchains, numsamples real(sp), pointer, dimension(:,:,:) :: data integer(i4b) :: l, status, blocksize, readwrite, numspec integer(i4b) :: fpixel, group, numargs logical(lgt) :: anyf real(sp) :: nullval character(len=80) :: comment integer(i4b), dimension(4) :: naxes real(sp), pointer, dimension(:,:,:,:) :: indata status = 0 readwrite = 0 nullval = 0. numargs = 1 ! Open the result file call ftopen(unit,trim(filename),readwrite,blocksize,status) ! Read keywords call ftgkyj(unit,'LMAX', lmax, comment,status) call ftgkyj(unit,'NUMSAMP', numsamples, comment,status) call ftgkyj(unit,'NUMCHAIN', numchains, comment,status) call ftgkyj(unit,'NUMSPEC', numspec, comment,status) allocate(data(0:lmax,numchains,numsamples)) allocate(indata(0:lmax,1,numchains,numargs+numsamples)) ! Read the binned power spectrum array group = 1 fpixel = 1 call ftgpve(unit,group,fpixel,size(indata),nullval,indata,anyf,status) call ftclos(unit,status) data = indata(:,1,:,numargs+1:numargs+numsamples) deallocate(indata) end subroutine read_gibbs_chain ! Utility routine for initializing the offset to be subtracted from each term ! to avoid overflow errors. Only called with the first power spectrum subroutine compute_largest_term(cls) implicit none real(dp), dimension(lmin:), intent(in) :: cls integer(i4b) :: i, j, l real(dp) :: subtotal, x ! Compute the Blackwell-Rao estimator offset = -1.6375e30 do i = 1, numchain do j = 1, numsamples subtotal = 0.d0 do l = lmin, lmax x = sigmas(l,i,j)/cls(l) subtotal = subtotal + & & 0.5d0 * real(2*l+1,dp) * (-x + log(x)) - log(real(sigmas(l,i,j),dp)) end do offset = max(offset,subtotal) end do end do end subroutine compute_largest_term end module br_mod_dist