COSMOMC: all_l_exact data format

Use of Healpix, camb, CLASS, cosmomc, compilers, etc.
Tor Helge Huse
Posts: 6
Joined: December 15 2004
Affiliation: Astrofysisk Institutt, Oslo

COSMOMC: all_l_exact data format

Post by Tor Helge Huse » April 18 2005

I'm trying to use the all_l_exact data format in cosmomc with some simulated data, but I keep running into computational errors when doing so. As I can't find any documentation on what the all_l_exact data format should be, I'll ask here: What does the all_l_exact data format actually look like?

regards
Tor Helge Huse

PS: I want to first use CAMB to generate some data and then introduce some noise through my own routine, so CAMB data is my starting point.

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

Re: COSMOMC: all_l_exact data format

Post by Antony Lewis » April 18 2005

The all_l_exact format uses the exact full-sky likelihood function given in e.g. Eq. 5 of astro-ph/0502469.

A .dataset file for an all_l_exact simulation would look something like

Code: Select all

name = Planck
has_pol = T
all_l_exact=T
all_l_file = data/Planck.dat
all_l_lmax=2000
Then you just need to make the Planck.dat file with the mock data. The columns in this file are

l C_TT (C_TE C_EE [C_BB]) N_T (N_P) f_sky^eff

The brackets are included if has_pol = T, and C_BB if you have compiled cosmomc with num_cls=4.

C_TT (and the polarized ones similarly) are raw C_l + noise (no l(l+1) or 2\pi factors). i.e. they increase on scales below the beam size due to the larger noise there.

Here's some simple code to make a Planck.dat file from a C_l file output by CAMB (in \mu\rm{K}^2):

Code: Select all

!Simulate CMB C_l data
!Simple version assuming isotropic noise and full sky
!no lensing, no tensors (check input Cl file has only 4 columns!)
!all parameters hard wired in code
!Data points are max likelihood values
!AL: Sept 2004

program CMB_cl_sim
 use AMLutils
 implicit none

!parameters
 character(LEN=*), parameter :: cl_infile = 'MAPcls.dat'
 character(LEN=*), parameter :: out_name = 'planck.dat'
 logical :: DoPol = .true., DoTemp = .true.
 real :: NoiseVar = 3d-4   !Pessimistic Planck 
       !in MicroK^2
 real, parameter :: fwhm_arcmin = 7. ! 13.  !WMAP ~ 13 (0.22 deg min) 
 real, parameter :: ENoiseFac = 2 !factor more noise on the polarization
 integer, parameter :: lmax = 2250


!Other stuff
 integer, parameter ::  C_T = 1, C_E = 2, C_B = 3, C_C = 4
 
 real :: fwhm_deg = fwhm_arcmin/60 
 real :: Cl(lmax,4), NoiseCl(lmax)
 real amp, xlc, sigma2
 integer l,ll
 real T, E, TE,B, scal


print *,'lmax = ',lmax
print *,'fwhm_arcmin = ', fwhm_arcmin
print *,'Noise = ',NoiseVar

   xlc= 180*sqrt(8.*log(2.))/3.14159
   sigma2 = (fwhm_deg/xlc)**2
   do l=2, lmax
     NoiseCl(l) = NoiseVar*exp(l*(l+1)*sigma2)
   end do


Cl=0
B=0
open(unit=1,file=Cl_infile,form='formatted',status='old')
do ll=2,lmax
  read (1,*) l, T, E, TE
  if (l/=ll) stop 'cl must start at l=2'
  scal=twopi/(l*(l+1))
  Cl(l,1) = T*scal
  Cl(l,2) = E*scal
  Cl(l,3) = B*scal
  Cl(l,4) = TE*scal
end do
close(1)

call CreateTxtFile(out_name,34)

 do l=2, lmax

      write (34,'(1I5,6E15.5)') l,(Cl(l,1)+NoiseCl(l)) , &
                                (Cl(l,C_C)), &
                                (Cl(l,C_E)+ENoiseFac*NoiseCl(l)), &
                            !      (Cl(l,C_B)+ENoiseFac*Noise%Cl(l,1))*1e6, &
                                    NoiseCl(l)*1e6,ENoiseFac*NoiseCl(l), 1.
 end do
 close(34)
end program CMB_cl_sim
If you want to use a proper simulation rather than taking the data to be the exact model, you should use 1/(2l+1) \sum_m |a_{lm}|^2 in the output, where a_{lm} is a Gaussian generated with variance C_l + N_l, or general a_{lm} from the theory C_l only and then add your own noise.
Last edited by Antony Lewis on April 19 2005, edited 1 time in total.

Tor Helge Huse
Posts: 6
Joined: December 15 2004
Affiliation: Astrofysisk Institutt, Oslo

COSMOMC: all_l_exact data format

Post by Tor Helge Huse » April 18 2005

Thank you for the information.
I'm trying to include tensors and when I do I get a very high Cl_BB term (~300 for l=20) for the first 100 or so l values, then dropping. Is this as it should be? The N_T term is also very large (same order of magnitude as Cl_BB).
When running with the generated data the likelihood values get very large (xE+9) and it does not seem to get anywhere.

Just to clarify: should COBE_normalize be set to true in CAMB?

regards

Tor Helge

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

Re: COSMOMC: all_l_exact data format

Post by Antony Lewis » April 18 2005

I have not tried it with tensors. Obviosuly the relative size of the noise depends on your experiment parameter, though for Planck it is correct that for sensible tensor amplitudes the signal will be noise dominated.

I wouldn't recommend using COBE_normalize, though you can as long as you make sure the final C_l is in \mu\rm{K}^2.

The noise free BB C_l should be about 0.001 the size of the temperature ones, so maybe you have a bug (not using the right column of the totCls output file or something: BB is column 4 of an output file with BB in).

Tor Helge Huse
Posts: 6
Joined: December 15 2004
Affiliation: Astrofysisk Institutt, Oslo

COSMOMC: all_l_exact data format

Post by Tor Helge Huse » April 19 2005

Code: Select all

 do l=2, lmax

      write (34,'(1I5,6E15.5)') l,(Cl(l,1)+NoiseCl(l)) , &
                                (Cl(l,C_C)), &
                                (Cl(l,C_E)+ENoiseFac*NoiseCl(l)), &
                            !      (Cl(l,C_B)+ENoiseFac*Noise%Cl(l,1))*1e6, &
                                    NoiseCl(l)*1e6,ENoiseFac*NoiseCl(l), 1.
 end do
 close(34)
end program CMB_cl_sim
I'm a bit confused when it comes to the above. Why is the C_l(B) term and NoiseCl(l) term multiplied by 1e6? This blows the numbers up a lot, and especially the noise term, and as far as I can that term is used unchanged in cosmomc.

And the line

Code: Select all

 !      (Cl(l,C_B)+ENoiseFac*Noise%Cl(l,1))*1e6, &
is for including tensors right? Is the term

Code: Select all

Noise%Cl(l,1)
correct? I recieve a compilation error as it is not defined anywhere. Should it be

Code: Select all

NoiseCl(l)
instead?

Tor Helge

Edit:
I wouldn't recommend using COBE_normalize, though you can as long as you make sure the final Cl is in μK2.
Now running with a fresh install of CAMB (april05), nothing changed except l_max, k_max and n_s (0.96), n_t(-0.04) and amp_ratio(0.248).
The noise free BB Cl should be about 0.001 the size of the temperature ones, so maybe you have a bug (not using the right column of the totCls output file or something: BB is column 4 of an output file with BB in).
The C_l(BB) is approx. 0.00001*C_l(TT) in the output from CAMB (column 4 vs column 2, l=48).
After running the data through your routine C_l(BB) is approx. 0.0001*C_l(TT)(l=48).

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

Re: COSMOMC: all_l_exact data format

Post by Antony Lewis » April 19 2005

Yes, sorry, I missed deleting those 1e6 factors when converting from my code (which originally used Healpix-normalized output in m\rm{K}^2). Simiarly with the incompatible name. Obviously you need to replace CreateTxtFile as well unless you are using my utils library.

Should be

Code: Select all

!Simulate CMB C_l data
!Simple version assuming isotropic noise and full sky
!no lensing, no tensors (check input Cl file has only 4 columns!)
!all parameters hard wired in code
!Data points are max likelihood values
!AL: Sept 2004

program CMB_cl_sim
 use AMLutils
 implicit none

!parameters
 character(LEN=*), parameter :: cl_infile = 'MAPcls.dat'
 character(LEN=*), parameter :: out_name = 'planck.dat'
 logical :: DoPol = .true., DoTemp = .true.
 real :: NoiseVar = 3d-4   !Pessimistic Planck 
       !in MicroK^2
 real, parameter :: fwhm_arcmin = 7. ! 13.  !WMAP ~ 13 (0.22 deg min) 
 real, parameter :: ENoiseFac = 2 !factor more noise on the polarization
 integer, parameter :: lmax = 2250


!Other stuff
 integer, parameter ::  C_T = 1, C_E = 2, C_B = 3, C_C = 4
 
 real :: fwhm_deg = fwhm_arcmin/60 
 real :: Cl(lmax,4), NoiseCl(lmax)
 real amp, xlc, sigma2
 integer l,ll
 real T, E, TE,B, scal


print *,'lmax = ',lmax
print *,'fwhm_arcmin = ', fwhm_arcmin
print *,'Noise = ',NoiseVar

   xlc= 180*sqrt(8.*log(2.))/3.14159
   sigma2 = (fwhm_deg/xlc)**2
   do l=2, lmax
     NoiseCl(l) = NoiseVar*exp(l*(l+1)*sigma2)
   end do


Cl=0
B=0
open(unit=1,file=Cl_infile,form='formatted',status='old')
do ll=2,lmax
  read (1,*) l, T, E, TE
  if (l/=ll) stop 'cl must start at l=2'
  scal=twopi/(l*(l+1))
  Cl(l,1) = T*scal
  Cl(l,2) = E*scal
  Cl(l,3) = B*scal
  Cl(l,4) = TE*scal
end do
close(1)

call CreateTxtFile(out_name,34)

 do l=2, lmax

      write (34,'(1I5,6E15.5)') l,(Cl(l,1)+NoiseCl(l)) , &
                                (Cl(l,C_C)), &
                                (Cl(l,C_E)+ENoiseFac*NoiseCl(l)), &
                            !      (Cl(l,C_B)+ENoiseFac*NoiseCl(l)), &
                                    NoiseCl(l),ENoiseFac*NoiseCl(l), 1.
 end do
 close(34)
end program CMB_cl_sim

Tor Helge Huse
Posts: 6
Joined: December 15 2004
Affiliation: Astrofysisk Institutt, Oslo

COSMOMC: all_l_exact data format

Post by Tor Helge Huse » April 19 2005

Ok, thank you. That makes a whole lot more sense. Testing it now (fingers crossed).

Tor Helge

Samuel Leach
Posts: 16
Joined: October 15 2004
Affiliation: SISSA, Trieste
Contact:

COSMOMC: all_l_exact data format

Post by Samuel Leach » May 11 2005

Dear Antony,

I've been following your discussion here and have been trying out your all_l_exact data format with CosmoMC which of course works rather superbly. Just a couple of related questions here that I'd like you to help me think through:

- I can see that using the exact C_l as the data is a Gaussian realisation, and a rather special one which can be useful for testing the mechanics of the data analysis pipeline for intrinsic bias in parameter reconstruction and in the determination of derived parameters, for instance. Do you know whether this will always lead to the right assessment for the posterior volume ? Of course I should just implement the Gaussian realisation version of the code to experiment a bit.

- What about your Planck noise model ? Can we set up a small database on CosmoCoffee which contains the current broad expectation for Planck and other future experiments, in terms of noise variance and beam width?

- What would be the standard way to cite all the useful information, codes etc from CosmoCofeee ? Did you discuss this together already? These aren't exactly Private Communications.

Thanks always for your time and energies Antony.

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

Re: COSMOMC: all_l_exact data format

Post by Antony Lewis » May 11 2005

Samuel Leach wrote:- I can see that using the exact C_l as the data is a Gaussian realisation, and a rather special one which can be useful for testing the mechanics of the data analysis pipeline for intrinsic bias in parameter reconstruction and in the determination of derived parameters, for instance. Do you know whether this will always lead to the right assessment for the posterior volume ? Of course I should just implement the Gaussian realisation version of the code to experiment a bit.
This is an interesting question, and I don't really know the answer. It probably depends on what you think the "right" answer is for forecasting work when you don't know what the exact realisation will be.

My experience for the CMB is that the error bars come out very similar to those you get from a typical realisation (with the advantage that the peak should be exactly at the input parameter values). However, I've not been able to prove any relation to any exact Bayesian forcasting protocol -- I'd be interested to know any thoughts people have.

Note that that actual quadrupole appears to be low, and hence not very typical. Error bars can depend slightly on the value of the quadrupole realisation, which is an example where using best-fit theory value may not be realistic for future experiment error bars (given what we already know about the realisation from WMAP).
Samuel Leach wrote: - What about your Planck noise model ? Can we set up a small database on CosmoCoffee which contains the current broad expectation for Planck and other future experiments, in terms of noise variance and beam width?
Sarah started work on a Wiki giving details of future experiments etc. This would be useful, but of course someone has to put in the work to input all the data. (though the idea with a Wiki of course is that anyone can do it)
Samuel Leach wrote: - What would be the standard way to cite all the useful information, codes etc from CosmoCofeee ? Did you discuss this together already? These aren't exactly Private Communications.
I've referred to plots on CosmoCoffee - just cited the URL. Otherwise I guess mention CosmoCoffee in the Acknowledgements - but we haven't really discussed this.

Pascal Vaudrevange
Posts: 50
Joined: March 26 2006
Affiliation: DESY
Contact:

COSMOMC: all_l_exact data format

Post by Pascal Vaudrevange » May 23 2007

Reviving this old thread.. I would like to create some mock datasets with WMAP3 accuracy. What would the values for noise_var be? Thanks

Pascal Vaudrevange
Posts: 50
Joined: March 26 2006
Affiliation: DESY
Contact:

COSMOMC: all_l_exact data format

Post by Pascal Vaudrevange » May 24 2007

I guess I can answer the question myself. I found this paper by Antony

http://arxiv.org/abs/astro-ph/0502469

which uses N_l^{TT}=0.03\mu K^2 . (Of course, this is only a crude approximation as the beam is not Gaussian..)

Jochen Weller
Posts: 42
Joined: September 24 2004
Affiliation: Ludwig-Maximilians-University Munich
Contact:

COSMOMC: all_l_exact data format

Post by Jochen Weller » March 08 2008

I have a question regarding the input format of the exact Cl mock data.
In cmbdata.F90 you write:

Code: Select all

  !File format:
     !l C_TT (C_TE C_EE [C_BB]) N_T (N_P) fsky_eff 
     !Must have num_cls set to correct value for file
         do l = 2, aset%all_l_lmax
            read (file_unit, *, end=100, err=100) idum,inobs(1:ncls), aset%all_l_noise(l,:), 
            aset%all_l_fsky(l)
           if (idum /= l) stop 'Error reading all_l_file'
            !set BB to pure noise if not in file
          if &#40;aset%has_pol .and. ncls < num_cls&#41; inobs&#40;num_cls&#41; = aset%all_l_noise&#40;l,2&#41;
          aset%all_l_obs&#40;l,&#58;&#41; = inobs&#40;1&#58;num_cls&#41;
        end do

However in the example code to create the mock data you write:

Code: Select all

     do l=2, lmax

      write &#40;34,'&#40;1I5,6E15.5&#41;'&#41; l,&#40;Cl&#40;l,1&#41;+NoiseCl&#40;l&#41;&#41; , &
                                &#40;Cl&#40;l,C_C&#41;&#41;, &
                                &#40;Cl&#40;l,C_E&#41;+ENoiseFac*NoiseCl&#40;l&#41;&#41;, &
                            !      &#40;Cl&#40;l,C_B&#41;+ENoiseFac*NoiseCl&#40;l&#41;&#41;, &
                                    NoiseCl&#40;l&#41;,ENoiseFac*NoiseCl&#40;l&#41;, 1.
       end do 
[q]
Now I am a bit puzzled by why the second column has to be written Cl(l,1)+NoiseCl(l)
I would have expected to have Cl(l,C_Temp) properly normalized to be in there, from the above read statement?[/q]

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

Re: COSMOMC: all_l_exact data format

Post by Antony Lewis » March 08 2008

I'm not sure I get the question, but all C_l here are the bare ones without factors of l(l+1)/2\pi.

[q]C_Temp=1.[/q]

Jochen Weller
Posts: 42
Joined: September 24 2004
Affiliation: Ludwig-Maximilians-University Munich
Contact:

COSMOMC: all_l_exact data format

Post by Jochen Weller » March 08 2008

What I meant it looks like in the 2nd column you are writing Cl+noise, but you just read Cl

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

Re: COSMOMC: all_l_exact data format

Post by Antony Lewis » March 08 2008

Yes, the C_l in the data file are the 'observed' C_l and include noise. The exact likelihood is a function of only C_l + N_l for full sky isotropic noise.

Post Reply