lesgourg / class_public

Public repository of the Cosmic Linear Anisotropy Solving System (master for the most recent version of the standard code; GW_CLASS to include Cosmic Gravitational Wave Background anisotropies; classnet branch for acceleration with neutral networks; ExoCLASS branch for exotic energy injection; class_matter branch for FFTlog)
230 stars 285 forks source link

Guidance with adding flat logarithmic parameter constraints #376

Closed astrogirl1 closed 3 years ago

astrogirl1 commented 3 years ago

Hi all, I am trying to run MontePython and CLASS together to get an MCMC analysis for the decaying dark matter scenario. I use a '.param' file to supply the values to the program.

I have introduced a new parameter : new_parameter in CLASS and it all works fine. The trouble I am facing is, my new parameter is used as it is in the equations I need, but the constraints I have are in the flat logarithmic form, in the form of : log_10{new_parameter}. I am unsure how to add this log_10{new_parameter} = [values] in my '.param' file and convert to new_parameter for the equations.

I am happy to supply more specific information.

ThomasTram commented 3 years ago

Hi MontePython does not have the option of "logarithmic sampling" of a parameter. But what you can do, (and what people usually do) is to define a new input parameter in CLASS called log10_new_parameter, read it, and then store it in new parameter. That way, CLASS uses new_parameter internally but the sampling can be in log10_new_parameter. If you want to make it nice, you should still be able to take new_parameter as input, and you should also make sure that one cannot accidentally pass both new_parameter and log10_new_parameter. You can adopt the following lines from input.c to accomplish that:

      class_call(parser_read_double(pfc,"A_s",&param1,&flag1,errmsg),
                 errmsg,
                 errmsg);
      class_call(parser_read_double(pfc,"ln10^{10}A_s",&param2,&flag2,errmsg),
                 errmsg,
                 errmsg);
      class_test((flag1 == _TRUE_) && (flag2 == _TRUE_),
                 errmsg,
                 "In input file, you cannot enter both A_s and ln10^{10}A_s, choose one");
      if (flag1 == _TRUE_)
        ppm->A_s = param1;
      else if (flag2 == _TRUE_)
        ppm->A_s = exp(param2)*1.e-10;

Cheers, Thomas

pstoecker commented 3 years ago

Hi,

I wanted to add an alternative solution to what Thomas is proposing. Whereas the solution of Thomas requires you to modify and recompile class you could also just add log10_new_parameter to the for loop starting around line 789 of montepython/data.py

# For all elements in the cosmological parameters from the mcmc list,
# translate any-one that is not directly a CLASS parameter into one.
# The try: except: syntax ensures that the first call
for elem in self.get_mcmc_parameters(['cosmo']):
     # infer h from Omega_Lambda and delete Omega_Lambda
     if elem == 'Omega_Lambda':
         omega_b = self.cosmo_arguments['omega_b']
         omega_cdm = self.cosmo_arguments['omega_cdm']
         Omega_Lambda = self.cosmo_arguments['Omega_Lambda']
         self.cosmo_arguments['h'] = math.sqrt(
             (omega_b+omega_cdm) / (1.-Omega_Lambda))
         del self.cosmo_arguments[elem]
     # infer omega_cdm from Omega_L and delete Omega_L
     elif elem == 'Omega_L':
          ...

This loops through all cosmological parameters and replaces them in case they are defined in one of the elif elem == "my_fancy_parameter" blocks with something reasonable and deleting the old one.

In your case you will need to add the following elif block

for elem in self.get_mcmc_parameters(['cosmo']):
    ...
    elif elem == 'log10_new_parameter':
        self.cosmo_arguments['new_parameter'] = 10.0 ** (self.cosmo_arguments[elem])
        del self.cosmo_arguments[elem]

Note that this will only allow you to use log10_new_parameter in Montepython but not in CLASS whereas the proposed solution of Thomas will also allow you to use the parameter natively in CLASS. In some sense both solution are equivalent; they just differ in where exactly you translate log10_new_parameter into new_parameter.

Best, Patrick

ThomasTram commented 3 years ago

Hi Patrick

Thanks a lot for that alternative suggestion!

Cheers, Thomas

astrogirl1 commented 3 years ago

Hi thank you so much for both the answers, I have tried both and they have worked! 👍