CosmoMC - running generic sampler on 2013 version

[Russell W Johnston]

Hey there

I am trying to run the generic sampler on the latest version of cosmomc.
The last version I ran used 'num_hard' in settings.f90 to set the number of parameters. However, I see this is now redundant.

When I try to run the program I get the following error:

Code: Select all

file_root:FR/chains
 Random seeds:  9367,  2262 rand_inst:   1
 Change the Parameterization type implementation to use Parameter names
 Computing tensors: F
 Doing CMB lensing: F
 Doing non-linear Pk:           0
 lmax              = 6500
 lmax_computed_cl  =   10
 max_eta_k         =    20.0000000000000     
 transfer kmax     =   0.800000011920929     
 Number of C_ls =    4
 parameter ranges not found: 1

I am trying to run with only 5 parameters and using the 'params_generic.ini' as a template. So in that file I have set:

Code: Select all

param1 = 0. -10  10. 0.01 0.01
param2 = 0. -10  10. 0.01 0.01
param3 = 0. -10  10. 0.01 0.01
param4 = 0. -10  10. 0.01 0.01
param5 = 0. -10  10. 0.01 0.01

use_CMB = F
use_HST = F
use_mpk = F
use_clusters = F
use_BBN = F
use_Age_Tophat_Prior = F
use_SN = F
use_min_zre = 0

also, in calclike.f90 i have:

Code: Select all

 function GenericLikelihoodFunction(Params)
    type(ParamSet)  Params
    real(mcp) :: GenericLikelihoodFunction
    real(mcp), allocatable, save :: covInv(:,:)
    real(mcp) X(num_params_used)

    alpha1 = Params%P(1)
    alpha2 = Params%P(2)
    beta1  = Params%P(3)
    beta2  = Params%P(4)
    beta3  = Params%P(5)

   GenericLikelihoodFunction = sfr_like()

   call MpiStop('GenericLikelihoodFunction: need to write this function!')
end function GenericLikelihoodFunction

So basically, its not clear to me how you tell cosmomc how many paramaters you are using in generic mode and thus why I get error message:

parameter ranges not found: 1

I presume its to do with .paranames file but not sure how you link it. Any help would be appreciated!

cheers

[Antony Lewis]

Yes, I must admit the generic interface has not been tested since the major planck update.

Try this: In GeneralTypes.f90 change

Code: Select all

        if (generic_mcmc) then
            Names%nnames=0
            if (Feedback>0) write (*,*) 'Change the Parameterization type implementation to use Parameter names'
        else

to

Code: Select all

        if (generic_mcmc) then
            num_theory_params= Ini_Read_Int_File(Ini, 'num_theory_params')
            Names%nnames=0
            if (Feedback>0) write (*,*) 'Change the Parameterization type implementation to use Parameter names'
        else

then add

Code: Select all

num_theory_params=5

to your input .ini file? Alternative (I think) you could code up a TParameterization subclass if you want to use named parameters.

[Jason Dossett]

I also got the following to work.

In driver.F90 change

Code: Select all

    call SetTheoryParameterization(DefIni, NameMapping)
    call DataLikelihoods%AddNuisanceParameters(NameMapping)
    call CMB_Initialize(Params%Info)

to

Code: Select all

    if(generic_mcmc)then
        num_theory_params= Ini_Read_Int('num_theory_params') 
    else
        call SetTheoryParameterization(DefIni, NameMapping)
        call DataLikelihoods%AddNuisanceParameters(NameMapping)
        call CMB_Initialize(Params%Info)
    end if

As with Antony's fix you will need to add the num_theory_params line to your input file.

[Russell W Johnston]

Yes, I must admit the generic interface has not been tested since the major planck update.

Try this: In GeneralTypes.f90 change

Code: Select all

        if (generic_mcmc) then
            Names%nnames=0
            if (Feedback>0) write (*,*) 'Change the Parameterization type implementation to use Parameter names'
        else

to

Code: Select all

        if (generic_mcmc) then
            num_theory_params= Ini_Read_Int_File(Ini, 'num_theory_params')
            Names%nnames=0
            if (Feedback>0) write (*,*) 'Change the Parameterization type implementation to use Parameter names'
        else

then add

Code: Select all

num_theory_params=5

to your input .ini file? Alternative (I think) you could code up a TParameterization subclass if you want to use named parameters.

Hi Antony

Thanks for getting back so quickly and for your help. The slight modification to the code seems to work like a charm. I'm not trying to re-invent the wheel, but for those interested in using cosmoMC purely as a generic sampler and are a bit unfamiliar, I would do the following steps:

• 1. Use the 'params_generic.ini' as a template. As Antony suggests, put somewhere in this file
  

  Code: Select all

  num_theory_params=number of parameters you want

  Set each parameter like e.g. :
  
  

  Code: Select all

  param1 = 0. -10  10. 0.01 0.01 
  param2 = 0. -10  10. 0.01 0.01 
  param3 = 0. -10  10. 0.01 0.01 
  param4 = 0. -10  10. 0.01 0.01 
  param5 = 0. -10  10. 0.01 0.01 
  

  I also set these to false, but not sure if it really matters:

  Code: Select all

  use_CMB = F 
  use_HST = F 

  2. In 'calclike.f90' I edited in the following way:
  

  Code: Select all

  function GenericLikelihoodFunction(Params)
      type(ParamSet)  Params
      real(mcp) :: GenericLikelihoodFunction
      real(mcp), allocatable, save :: covInv(:,:)
      real(mcp) X(num_params_used)
          
         alpha1 = Params%P(1)
         alpha2 = Params%P(2)
         beta1  = Params%P(3)
         beta2  = Params%P(4)
         beta3  = Params%P(5)
  
        GenericLikelihoodFunction = my_likelihood_function()
  end function GenericLikelihoodFunction

  (I have a module where I define all my parameter names)
  
  3. In settings.f90 set

  Code: Select all

  generic_mcmc=.true.

  4. In GeneralTypes.f90 change
  

  Code: Select all

          if (generic_mcmc) then
              Names%nnames=0
              if (Feedback>0) write (*,*) 'Change the Parameterization type implementation to use Parameter names'
          else

  to
  

  Code: Select all

          if (generic_mcmc) then
              num_theory_params= Ini_Read_Int_File(Ini, 'num_theory_params')
              Names%nnames=0
              if (Feedback>0) write (*,*) 'Change the Parameterization type implementation to use Parameter names'
          else

That's it I think. Thanks again Antony!!!

Cheers
Russ