Cobaya: How to use tabulated w(a) from CAMB?

[Dror Berechya]

Hello all,

I am a new user of Cobaya and CAMB, and I face some difficulties:

• The CAMB documentation says one can use the "DarkEnergyEqnOfState" class with a tabulated w(a). However, when I try to give "a" and "w" an array of values (by adding, e.g., "w: [#, #, ..., #]" into the input YAML file to run with "cobaya-run"), I get an error message saying that this parameter can take only numerical values (even when I use "use_tabulated_w: true"). How can I define the tabulated w(a) and feed it into the Cobaya as arrays of values for "a" and "w"?

• I noticed that the "DarkEnergyEqnOfState" class has the "set_w_a_table(a, w)" function. How can I use CAMB's (or Cobaya's) built-in Python functions and methods in Cobaya's input YAML file (to run with Cobaya-run)?

• Is there any advantage of running Cobaya from a Python interpreter compared to running it from Shell with cobaya-run and a YAML file?

• Returning to question number one, can I set up a tabulated w(a) that depends on a sampled parameter? For example, can I set "a: [1, 0.1, 0.001]" and "w: [-x, -0.7x, -0.5x]" where "x" is a sampled parameter?

• Lastly, are there published guidelines on modifying CAMB to define, e.g., new dark energy classes (new dark energy models)? Is there documentation also for the Fortran side of CAMB? (Generically, will I need to modify also the Fortran code if the answer to the previous question requires defining a new class inheriting from "DarkEnergyEqnOfState"?) Are there examples showing how to use CAMB features with Cobaya? (similar to the CAMBdemo.ipynb but showing how to communicate CAMB to Cobaya.)

Thank you in advance for your response!

[Raphael Kou]

Hello,

I am afraid the tabulated w(a) cannot be directly used in Cobaya with the current implementation of CAMB. When the dark energy module is set up, w is expected to be a float, which is the reason why you get the error. If you want to use the tabulated w(a) with Cobaya, you would need to modify the code. One way of doing that is by adding two arrays a and w in argument of the set_params function of dark_energy.py and then to modify that function in order to call the set_w_a_table function.

After having implemented those modifications, you will be able to use the a and w arrays in the yaml file you give to Cobaya. If your w array depends on one x sampled parameter, you should use the vector parameters approach: https://cobaya.readthedocs.io/en/latest/params_prior.html#vector-parameters. This is in any case needed, as Cobaya can only sample scalars.

The only advantage I can think of, running Cobaya from a Python interpreter is when you want to use the model object (https://cobaya.readthedocs.io/en/latest/models.html) to interact with the sampler, but there could have others.

Ultimately, you can find guidelines on modifying CAMB here: https://camb.readthedocs.io/en/latest/modifying_code.html. Those guidelines should help you modifying both the python and fortran codes.

I hope this helps and answers most of your questions.

[Dror Berechya]

Thank you so much for your answer; it is of great help!

For clarification, by saying,

If you want to use the tabulated w(a) with Cobaya, you would need to modify the code.

Do you mean I need to modify CAMB or Cobaya? (I think CAMB, as I didn't find "dark_energy.py" in Cobaya's source code.)

Thank you very much again.

[Raphael Kou]

Yes indeed, I mean modifying CAMB.

[Dror Berechya]

Thank you very much for your valuable aid!

I modified the set_parms() function of CAMB's dark_energy.py like so:

Code: Select all

def set_params(self, a=[0.2,0.4,0.6,0.8,1], b=[-1,-0.9,-0.8,-0.9,-1], w=-1, wa=0, cs2=1.0):
        """
         Set the parameters so that P(a)/rho(a) = w(a) = w + (1-a)*wa

        :param w: w(0)
        :param wa: -dw/da(0)
        :param cs2: fluid rest-frame sound speed squared
        """
        self.w = w
        self.a = a
        self.b = b
        self.wa = wa
        self.cs2 = cs2
        self.validate_params()
        self.set_w_a_table(a,b)

I named the scale factor and the EoS parameter as "a" and "b" because I didn't want to mess up with any of CAMB's original parameters. This modification seems to take effect only when I provide "b: [#,#,...,#]" in the "extra arguments" section of Cobaya's input YAML file. However, it doesn't work if I try to put the "b" list in the "params" section of the YAML file, which is needed if I want to solve question number 4 (in the original post) with the "Vector parameters" approach. (Also, I don't understand why it doesn't read the default for "b" I specified in the set_params() definition but only the default for "a".)

Do you know what the problem is and what I should do?

Conceptual question: When Cobaya runs, how does it call CAMB's parameters and functions? Does it call something like "DarkEnergy.set_params()"?

Thank you so much for the help,
Dror

[Raphael Kou]

I am surprised it doesn't work. What is the error you get? Following your fourth question in the original message, did you try implementing the vector b as:

Code: Select all

"x": {"prior": {"min": 0, "max": 1}, "proposal": 0.1, "drop": True},
"b": {"value": lambda x: [-x, -0.7*x, -0.5*x], "derived": False}

What do you mean by "I don't understand why it doesn't read the default for "b""? If you add print(b) in the set_params function, is b not set?

About your last question, there is a camb.py file in cobaya/theories/camb where CAMB's parameters and functions are called, including the set_params functions.

[Dror Berechya]

Thank you for your response!

To clarify, I get results that differ from the LCDM results when I specify "b" in the "extra_args" section of the input YAML file like so:

Code: Select all

theory:
  camb:
    path: .../cobaya/cosmo_packages/code/CAMB_modified_2
    extra_args:
      halofit_version: mead
      bbn_predictor: PArthENoPE_880.2_standard.dat
      lens_potential_accuracy: 1
      num_massive_neutrinos: 1
      nnu: 3.044
      theta_H0_range:
      - 20
      - 100
      b: [-1,-0.9,-0.8,-0.9,-1]

However, if I don't specify "b" in the YAML file, the results I get are identical to the LCDM results, even though I have "def set_params(self, a=[0.2,0.4,0.6,0.8,1], b=[-1,-0.9,-0.8,-0.9,-1], w=-1, wa=0, cs2=1.0)" in the definition of set_params(). In any case, I don't need to specify "a" in the YAML file to get results that differ from LCDM.

• Why doesn't Cobaya read "b" from the default in the definition of set_params() such that I must manually specify "b" in the YAML file to get the expected results?

• And then, why does Cobaya read "a" from the default in the definition of set_params() such that I don't need to specify it in the YAML file?

Next, when I try to specify "b" in the "params" section of the YAML file, I get the following error:

[model] *ERROR* Could not find anything to use input parameter(s) {'b'}.

• Why does it happen?

Now, regarding the vector parameters, to run:

Code: Select all

theory:
  camb:
    path: .../cobaya/cosmo_packages/code/CAMB_modified_4
    extra_args:
      halofit_version: mead
      bbn_predictor: PArthENoPE_880.2_standard.dat
      lens_potential_accuracy: 1
      num_massive_neutrinos: 1
      nnu: 3.044
      theta_H0_range:
      - 20
      - 100
      x:
        prior:
          min: 0
          max: 1
        proposal: 0.1
        drop: true
      b:
        value: 'lambda x: [-1*x,-1*x,-1*x,-1*x,-1*x]'
        derived: false

I had to modify the code to accept a dictionary for b:

Code: Select all

def set_params(self, a=[0.2,0.4,0.6,0.8,1], b=[-1,-0.9,-0.8,-0.9,-1], w=-1, wa=0, cs2=1.0):
        """
         Set the parameters so that P(a)/rho(a) = w(a) = w + (1-a)*wa

        :param w: w(0)
        :param wa: -dw/da(0)
        :param cs2: fluid rest-frame sound speed squared
        """
        if isinstance(b, dict):
            b = b['value']
        else:
            b = b
        
        self.w = w
        self.a = a
        self.b = b
        self.wa = wa
        self.cs2 = cs2
        self.validate_params()
        self.set_w_a_table(a,b)

When I run it with "--test", I get:

[mcmc] Getting initial point... (this may take a few seconds)
[model] *ERROR* Could not find random point giving finite posterior after 1080 tries

Upon checking with "--debug", I think the problem is with the "lambda" function of "b":

2024-06-12 18:40:01,358 [camb.transfers] Ignored error at evaluation and assigned 0 likelihood (set 'stop_at_error: True' as an option for this component to stop here and print a traceback). Error message: ValueError("could not convert string to float: 'lambda x: [-1*x,-1*x,-1*x,-1*x,-1*x]'")
2024-06-12 18:40:01,358 [model] Calculation failed, skipping rest of calculations
2024-06-12 18:40:01,358 [model] *ERROR* Could not find random point giving finite posterior after 1080 tries

• Do you know what the problem is, and what should I do to fix it?

Lastly, I do have a camb.py file, but I didn't touch it:

Code: Select all

#!/usr/bin/env python

# Python command line CAMB reading parameters from a .ini file
# an alternative to fortran binary compiled into fortran/camb using "make camb".
# To use .ini files from a python script use camb.run_ini or camb.read_ini
# If you have installed the camb package, you can just use "camb params.ini" without using this script.

from camb._command_line import run_command_line

run_command_line()

• Should I modify this as well? What is this file, and what is its function?

Thanks a lot; your help is valued :)

[Antony Lewis]

Try "b=None" rather than setting to a default list? (can check inside function if b is None and assign default array if you want to).

Raphael meant camb.py in cobaya.

[Dror Berechya]

Thanks for the clarification. Yes, I have found the camb.py file in cobaya/theories/camb; does one need to modify it to make Cobaya read parameters for CAMB from the "params:" section in the YAML input file? If so, which part of this code does one need to modify?

My "set_params()" function in dark_energy.py (in cosmo_packages/code/CAMB/camb/) is defined as:

Code: Select all

def set_params(self, a=[0.2,0.4,0.6,0.8,1], x=0.5, w=-1, wa=0, cs2=1.0):
        """
         Set the parameters so that P(a)/rho(a) = w(a) = w + (1-a)*wa

        :param w: w(0)
        :param wa: -dw/da(0)
        :param cs2: fluid rest-frame sound speed squared
        """
        self.x = x
        
        b = [-0.2*x,-0.4*x,-0.6*x,-0.8*x,-1*x]
        
        self.w = w
        self.a = a
        self.wa = wa
        self.cs2 = cs2
        self.validate_params()
        self.set_w_a_table(a,b)

Now, I can call the "x" parameter with a prior within the "params:" section of the YAML file and the "a" parameter as a list only in the "extra_args:" section. It seems I can call list parameters only in the "extra_args:" section (as, e.g., "a: [#,...,#]") and numerical parameters in either the "extra_args:" or "params:" sections. For now, it is enough for my needs, but I would like to know what I should change to be able to include lists in the "params:" section.

Important: Is my modification of the "set_params()" function enough to define a new dark energy model, or do I need to modify other parts?

I have also found a camb.py file in cosmo_packages/code/CAMB/camb/; it has those parts in it:

Code: Select all

def do_set(setter):
        kwargs = {kk: params[kk] for kk in getfullargspec(setter).args[1:] if kk in params}
        used_params.update(kwargs)
        if kwargs:
            if verbose:
                logging.warning('Calling %s(**%s)' % (setter.__name__, kwargs))
            setter(**kwargs)

    # Note order is important: must call DarkEnergy.set_params before set_cosmology if setting theta rather than H0
    # set_classes allows redefinition of the classes used, so must be called before setting class parameters
    do_set(cp.set_accuracy)
    do_set(cp.set_classes)
    do_set(cp.DarkEnergy.set_params)
    do_set(cp.set_cosmology)
    do_set(cp.set_matter_power)
    do_set(cp.set_for_lmax)
    do_set(cp.InitPower.set_params)
    do_set(cp.NonLinearModel.set_params)

And

Code: Select all

def set_params_cosmomc(p, num_massive_neutrinos=1, neutrino_hierarchy='degenerate', halofit_version='mead',
                       dark_energy_model='ppf', lmax=2500, lens_potential_accuracy=1, inpars=None):
    """
    get CAMBParams for dictionary of cosmomc-named parameters assuming Planck 2018 defaults

    :param p: dictionary of cosmomc parameters (e.g. from getdist.types.BestFit's getParamDict() function)
    :param num_massive_neutrinos: usually 1 if fixed mnu=0.06 eV, three if mnu varying
    :param neutrino_hierarchy: hierarchy
    :param halofit_version: name of the soecific Halofit model to use for non-linear modelling
    :param dark_energy_model: ppf or fluid dark energy model
    :param lmax: lmax for accuracy settings
    :param lens_potential_accuracy: lensing accuracy parameter
    :param inpars: optional input CAMBParams to set
    :return:
    """
    pars = inpars or model.CAMBparams()
    if p.get('alpha1', 0) or p.get('Aphiphi', 1) != 1:
        raise ValueError('Parameter not currrently supported by set_params_cosmomc')
    pars.set_cosmology(H0=p['H0'], ombh2=p['omegabh2'], omch2=p['omegach2'], mnu=p.get('mnu', 0.06),
                       omk=p.get('omegak', 0), tau=p['tau'], deltazrei=p.get('deltazrei', None),
                       nnu=p.get('nnu', constants.default_nnu), Alens=p.get('Alens', 1.0),
                       YHe=p.get('yheused', None), meffsterile=p.get('meffsterile', 0),
                       num_massive_neutrinos=num_massive_neutrinos, neutrino_hierarchy=neutrino_hierarchy)
    pars.InitPower.set_params(ns=p['ns'], r=p.get('r', 0), As=p['A'] * 1e-9, nrun=p.get('nrun', 0),
                              nrunrun=p.get('nrunrun', 0))
    pars.set_dark_energy(w=p.get('w', -1), wa=p.get('wa', 0), dark_energy_model=dark_energy_model)
    pars.set_for_lmax(lmax, lens_potential_accuracy=lens_potential_accuracy)
    pars.NonLinearModel.set_params(halofit_version=halofit_version)
    pars.WantTensors = pars.InitPower.has_tensors()
    return pars

Do I need to change "do_set(cp.DarkEnergy.set_params)", "pars.set_dark_energy(w=p.get('w', -1), wa=p.get('wa', 0), dark_energy_model=dark_energy_model)", or anything else?

[Antony Lewis]

Did you try putting "b=None" as an argument to set_params for the list (with vector mapping set up in cobaya)?

[Dror Berechya]

Thank you, Antony!
My modification works fine for my needs (at least for the moment).
Now, is it enough to modify the "set_params()" function in cosmo_packages/code/CAMB/camb/ to define a new dark energy model with a tabulate Equation of State, or do I need to modify other parts? Is there a way to verify or ensure that Cobaya runs only with my tabulated dark energy model for all its calculations?

Thanks again,
Dror