Writing a Likelihood using CMB cross-correlation function

Use of Healpix, camb, CLASS, cosmomc, compilers, etc.
Post Reply
Esteban Chalbaud
Posts: 9
Joined: June 01 2019
Affiliation: Master's student
Contact:

Writing a Likelihood using CMB cross-correlation function

Post by Esteban Chalbaud » October 22 2019

Hello,

I am trying to write a likelihood for the first time and use it inside CosmoMC, so the idea is to use equation (4.11) in the following paper: 10.1088/1475-7516/2011/07/027 , but i am confused how to sample over the parameter beta which is inside the likelihood. What i have done so far is:

1) modify in CosmologyType.f90 the number of the parameter max_inipower_params from default value (10) into 11
integer, parameter :: max_inipower_params = 10

2) I declare in the same file specifically at the function: Type, extends(TTheoryParams) :: CMBParams
the name of the variable that i want to sample:
real(mcp) boost

3) I added in the paramnames file: params_CMB.paramnames the parameter that i am introducing (beta in this case)

After this i recompile everything, but when i try to run my chain. i get this error:
SetTheoryParameterNumbers: parameter numbers do not match
MpiStop: 0

I found this way of do this in the following presentation: http://icg.port.ac.uk/~jschewts/cantata/CAMB/CosmoMC_lecture.pdf
since i did not find any other good source. In general i am curious in how do you add a parameter that you want to sample inside a likelihood? i think i have to do the same on getdist to plot the probability distribution for that added parameter. If there is any better tutorial in how to do this Let me know.


Thanks in advance.

Antony Lewis
Posts: 1512
Joined: September 23 2004
Affiliation: University of Sussex
Contact:

Re: Writing a Likelihood using CMB cross-correlation function

Post by Antony Lewis » October 22 2019

You don't need to modify any of the theory types to use a parameter in a likelihood, everything should be inside your likelihood class (i.e. revert your current changes). Basically you call likelihood%loadParamNames('myparamsfile.paramnames') to assign your nuisance parameters. The values of the nuisance parameters are then passed to your likelihood%LogLike method. So all you need to do is specify the ranges in your .ini files as you would any other parameter.

For example see how it's done in supernovae_JLA and CMB_BK_Planck.

Esteban Chalbaud
Posts: 9
Joined: June 01 2019
Affiliation: Master's student
Contact:

Re: Writing a Likelihood using CMB cross-correlation function

Post by Esteban Chalbaud » October 25 2019

Ok now is working. For those interested in what i did here is my Likelihood :

________________________________________________________________________
module boost_estimator
use CosmologyTypes
use CosmoTheory
use Likelihood_Cosmology
implicit none
private

type, extends(TCosmoCalcLikelihood) :: boost_estimatorLikelihood
integer(kind=4) :: lmaxtt = 1024
real(mcp), allocatable :: ell(:)
real(mcp), allocatable :: ClsTT(:)
real(mcp), allocatable :: ClsTTcross(:)
real(mcp), allocatable :: ClsTE(:)
real(mcp), allocatable :: ClsTE_err(:)
real(mcp), allocatable :: ClsEE(:)
real(mcp), allocatable :: ClsEE_err(:)
contains

procedure :: boost_estimator_ReadCls
procedure :: LogLike => boost_estimator_LnLike

end type boost_estimatorLikelihood

public boost_estimatorLikelihood, boost_estimatorLikelihood_Add
contains

!Subroutine area ................... here we want to reade the .ini likelihood file from CosmoMC
subroutine boost_estimatorLikelihood_Add(LikeList, ini)
class(TLikelihoodList) LikeList
class(TSettingIni) Ini
Type (boost_estimatorLikelihood), pointer :: this
character(len=256) :: ClTT_filename = ''
character(LEN=:), allocatable :: string

if (Ini%Read_Logical('boost_use_CMB_TT_only',.false.)) then !EMS there should be use_CMB_only=T inside .ini file
allocate(this)

this%LikelihoodType = 'boost_estimator'
this%name = Ini%Read_String('boost_CMB_TT_only_name')
call Ini%Read('boost_CMB_TT_only_lmax',this%lmaxtt)
string = Ini%Read_String('boost_CMB_TT_only_fname')
if (string /= '') call this%boost_estimator_ReadCls(string)
this%needs_background_functions = .true.

call this%loadParamNames(trim(DataDir)//'boost.paramnames')

call LikeList%Add(this)
end if
end subroutine boost_estimatorLikelihood_Add

!Subroutine area ................... here we want to read the power spectrum

subroutine boost_estimator_ReadCls(this, ClTT_filename)
class(boost_estimatorLikelihood) :: this
character(LEN=:), allocatable :: ClTT_filename
integer j, idx, m
integer :: IO = 0
real(mcp) c1value
real(mcp) c2value
real(mcp) c3value
real(mcp) err
real(mcp) err2
real(mcp) err3

allocate ( this%ClsTT(this%lmaxtt+1) )
allocate ( this%ClsTTcross(this%lmaxtt+1) )
allocate ( this%ClsEE(this%lmaxtt+1) )
allocate ( this%ClsEE_err(this%lmaxtt+1) )
allocate ( this%ClsTE(this%lmaxtt+1) )
allocate ( this%ClsTE_err(this%lmaxtt+1) )

print*,'POWER ',trim(ClTT_filename)

open(unit = 100, file = ClTT_filename, status = 'old', action = 'read')

!do j = 0, this%lmaxtt
!this%ClsTT(j) = -1.d0
!this%ClsTTcross(j) = 0.d0
!this%ClsEE(j) = -1.d0
!this%ClsEE_Err(j) = 0.d0
!this%ClsTE(j) = -1.d0
!this%ClsTE_Err(j) = 0.d0
!end do

do while (IO == 0)
read(100, *, IOSTAT = IO) idx, c1value
if (IO == 0 .and. idx<=this%lmaxtt) then
this%ClsTT(idx) = c1value
!this%ClsTTcross(idx) = err
!print *, 'Power Spectra', idx, this%ClsTT(idx)
end if

end do

!call this%loadParamNames(trim(DataDir)//'boost.paramnames')

end subroutine boost_estimator_ReadCls


!We declare the function that calculates the likelihood a rela valued one
real(mcp) function boost_estimator_LnLike(this, CMB, Theory, DataParams) !the funtion takes parameters from othe modules
class(boost_estimatorLikelihood) :: this
class(CMBParams) CMB
class(TCosmoTheoryPredictions), target :: Theory
real(mcp) :: DataParams(:)
real(mcp) :: boost, theoryval1, theoryval2, noise1, noise2, observedf, coefficientA, coefficientB, theoryf
integer l
integer m
boost_estimator_LnLike = 0.d0
boost=DataParams(1)

associate(Cls=> Theory%Cls(1,1))
do l=0,this%lmaxtt ! loop over multipoles !skip C_0 and C_1
if (l >= 2 .and. this%ClsTT(l) > 0.) then
theoryval1 = CLs%CL(l)/(l*(l+ 1.d0))*2.d0*pi ! already multiplied by l*(l+1)/twopi
theoryval2 = CLs%CL(l + 1)/((l+1.d0)*(l+ 2.d0 ))*2.d0*pi ! already multiplied by l*(l+1)/twopi

noise1 = (14.d-6)**2*exp(l*(l+1.d0)*8.84)
noise2 = (14.d-6)**2*exp((l+1.d0)*(l+2.d0)*8.84)
observedf = this%ClsTTcross(l)

do m=-l,l ! loop over m !skip C_0 and C_1
coefficientA = l*sqrt((l+1.d0)**2-m**2)/sqrt(4*(l + 1.d0)**2-1.d0)
coefficientB = -(l+2)*sqrt((l+1.d0)**2-m**2)/sqrt(4*(l + 1.d0)**2-1.d0)
end do

theoryf = coefficientA*theoryval1 + coefficientB*theoryval2

boost_estimator_LnLike = boost_estimator_LnLike + &
((0.8)**2)*(observedf - (boost)*theoryf)**2 / ((theoryval1 + noise1)*(theoryval1 + noise1))
coefficientA = 0.d0
coefficientB = 0.d0
!print *, 'Power Spectra', boost
end if
end do
end associate
end function boost_estimator_LnLike
!ending the module
end module boost_estimator

________________________________________________________________________

But i still have a the MCMC stops and complains due to some negative Cls, there is something that i am missing? this is the error:
___________________________________________________________________________

Number of MPI processes: 2
file_root:boost_estimator_test
NOTE: num_massive_neutrinos ignored, using specified hierarchy
NOTE: num_massive_neutrinos ignored, using specified hierarchy
Random seeds: 18762, 22676 rand_inst: 2
Random seeds: 18662, 22676 rand_inst: 1
POWER /home/esteban/CosmoMC_test/data/homemade_like/power_cross_test.txt
POWER /home/esteban/CosmoMC_test/data/homemade_like/power_cross_test.txt
Doing non-linear Pk: F
Doing CMB lensing: T
Doing non-linear lensing: T
TT lmax = 2500
EE lmax = 2500
ET lmax = 2500
BB lmax = 2500
PP lmax = 2500
lmax_computed_cl = 2500
Computing tensors: F
max_eta_k = 14000.0000
transfer kmax = 1.00000000
adding parameters for: ECHM
Fast divided into 1 blocks
7 parameters ( 6 slow ( 0 semi-slow), 1 fast ( 0 semi-fast))
starting Monte-Carlo
Chain:0 drag accpt: 1.00000000 fast/slow 2.06666660 slow: 30
Chain:1 drag accpt: 1.00000000 fast/slow 2.20000005 slow: 30
Chain:0 drag accpt: 1.00000000 fast/slow 2.06666660 slow: 60
Chain:1 drag accpt: 1.00000000 fast/slow 2.23333335 slow: 60
Chain:0 drag accpt: 1.00000000 fast/slow 1.93333328 slow: 90
Chain:1 drag accpt: 1.00000000 fast/slow 2.31111121 slow: 90
Chain1, MPI done 'burn', Samples =118, like = 0.00000000
Time: 491.86910210100007 output lines= 116
slow changes 118 power changes 0
MPI_Min_Sample_Update 79 118
Chain2, MPI done 'burn', Samples =119, like = 0.00000000
Time: 500.03799697899990 output lines= 117
slow changes 119 power changes 0
MPI_Min_Sample_Update 79 119
Chain:0 drag accpt: 1.00000000 fast/slow 2.09999990 slow: 120
1 all_burn done
1 DoUpdates
Chain:1 drag accpt: 1.00000000 fast/slow 2.28333330 slow: 120
2 all_burn done
2 DoUpdates
Chain:0 drag accpt: 1.00000000 fast/slow 2.05333328 slow: 150
Chain:1 drag accpt: 1.00000000 fast/slow 2.40000010 slow: 150
Chain:0 drag accpt: 1.00000000 fast/slow 2.07777786 slow: 180
Chain:1 drag accpt: 0.994475126 fast/slow 2.48888898 slow: 180
Chain:1 drag accpt: 0.981308401 fast/slow 2.47619057 slow: 210
Chain:0 drag accpt: 1.00000000 fast/slow 2.00952387 slow: 210
Chain:1 drag accpt: 0.948616624 fast/slow 2.48333335 slow: 240
Chain:0 drag accpt: 1.00000000 fast/slow 1.98333335 slow: 240
Chain:1 drag accpt: 0.950704217 fast/slow 2.45185184 slow: 270
Chain:0 drag accpt: 1.00000000 fast/slow 2.04444456 slow: 270
Chain:1 drag accpt: 0.955413997 fast/slow 2.40000010 slow: 300
Chain:0 drag accpt: 1.00000000 fast/slow 2.09999990 slow: 300
Chain 2 MPI communicating
Chain 1 MPI communicating
Current convergence R-1 = 103.575783 chain steps = 322
slow changes 322 power changes 0
updating proposal density
Chain:1 drag accpt: 0.942857146 fast/slow 2.41212130 slow: 330
Chain:0 drag accpt: 0.985074639 fast/slow 2.10909081 slow: 330
Chain:1 drag accpt: 0.935064912 fast/slow 2.47777772 slow: 360
Chain:0 drag accpt: 0.957446814 fast/slow 2.22222233 slow: 360
Chain:1 drag accpt: 0.926365793 fast/slow 2.51282048 slow: 390
Chain:0 drag accpt: 0.955882370 fast/slow 2.26410246 slow: 390
Chain:1 drag accpt: 0.895522416 fast/slow 2.59761906 slow: 420
Chain:0 drag accpt: 0.952380955 fast/slow 2.31190467 slow: 420
Chain:1 drag accpt: 0.889328063 fast/slow 2.65555549 slow: 450
Chain:0 drag accpt: 0.955413997 fast/slow 2.39333344 slow: 450
Chain:1 drag accpt: 0.893854737 fast/slow 2.69374990 slow: 480
Chain:0 drag accpt: 0.958083808 fast/slow 2.46041656 slow: 480
Chain:1 drag accpt: 0.897887349 fast/slow 2.73333335 slow: 510
Chain:0 drag accpt: 0.956848025 fast/slow 2.48235297 slow: 510
Chain:1 drag accpt: 0.903010011 fast/slow 2.77407408 slow: 540
Chain:0 drag accpt: 0.955752194 fast/slow 2.54814816 slow: 540
Chain:1 drag accpt: 0.907643318 fast/slow 2.82456136 slow: 570
Chain:0 drag accpt: 0.957983196 fast/slow 2.58771920 slow: 570
Chain 2 MPI communicating
Chain 1 MPI communicating
Current convergence R-1 = 13.4338703 chain steps = 602
slow changes 577 power changes 0
updating proposal density
Chain:1 drag accpt: 0.906344414 fast/slow 2.83333325 slow: 600
Chain:0 drag accpt: 0.958466470 fast/slow 2.61333323 slow: 600
Chain:1 drag accpt: 0.907781005 fast/slow 2.85238099 slow: 630
Chain:0 drag accpt: 0.958904088 fast/slow 2.64603186 slow: 630
Chain:1 drag accpt: 0.901639342 fast/slow 2.88484859 slow: 660
Chain:0 drag accpt: 0.953757226 fast/slow 2.67121220 slow: 660
Calculator_CAMB: negative C_l (could edit to silent error here)
MpiStop: 1
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 1 in communicator MPI_COMM_WORLD
with errorcode 128.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.


___________________________________________________________________________

I am attaching the plot that i got using only my likelihood just before the MCMC stopped.

Image

Antony Lewis
Posts: 1512
Joined: September 23 2004
Affiliation: University of Sussex
Contact:

Re: Writing a Likelihood using CMB cross-correlation function

Post by Antony Lewis » October 25 2019

If you are using CMB spectra you need to set the this%cl_lmax() to the required CMB ranges that your likelihood uses.

Post Reply