CAMB: passing a new CAMBparams parameter to recfast.f90

Use of Cobaya. camb, CLASS, cosmomc, compilers, etc.
Post Reply
Joshua Benabou
Posts: 5
Joined: September 14 2022
Affiliation: University of California Berkeley

CAMB: passing a new CAMBparams parameter to recfast.f90

Post by Joshua Benabou » September 14 2022

Hello,

I have added a new parameter in model.py and model.f90 at the bottom of the CAMBparams in model.f90, via

Code: Select all

real(dl)  :: newparam = 0._dl. 
I would like to be able to have the user modify this parameter , via e.g modifying the setter set_cosmology to use

Code: Select all

pars = camb.CAMBparams()
pars.set_cosmology(newparam= 1)
This is easy enough. But I would then like the Recombination module recfast.f90 to take newparam as an input. I see that recfast.f90 has the subroutine TRecfast_ReadParams(this, Ini) which reads a few input parameters from some input file.

However I don't want to read from any file, I simply want to pass newparam from python directly, once the variable has been set, presumably at the level of recombination.py. How can I do this easily?

Thanks alot!

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

Re: CAMB: passing a new CAMBparams parameter to recfast.f90

Post by Antony Lewis » September 14 2022

You can add it to the TRecfast type instead (Recfast in python).

Joshua Benabou
Posts: 5
Joined: September 14 2022
Affiliation: University of California Berkeley

Re: CAMB: passing a new CAMBparams parameter to recfast.f90

Post by Joshua Benabou » September 14 2022

Dear Anthony,

Thanks for your reply.

Can you please elaborate? I have not yet understood how CAMB interfaces between python and fortran. Where is the TRecfast' type defined? Are you suggesting to add newparam as one of the fields in

Code: Select all

class Recfast(RecombinationModel)
?
If so, Im a bit confused it seems those fields are initialized in the fortran file, not in recombination.py. I did try to add

Code: Select all

 ("newparam", c_double)
to the end of _fields_ in class Recfast(RecombinationModel), and similarly add

Code: Select all

    real(dl) :: fa = 1D1 
at the end of the type, extends(TRecombinationModel) :: TRecfast block in recfast.f90. I can set this new parameter by modifying camb.CAMBparams() just fine, but calling

Code: Select all

 camb.get_results(pars)
causes my Jupyter notebook to crash, even though newparam is not being used anywhere yet.

Perhaps you could point me to a similar Python class which calls fortran elsewhere in CAMB which already does what im trying to do.

thanks much!

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

Re: CAMB: passing a new CAMBparams parameter to recfast.f90

Post by Antony Lewis » September 14 2022

Should be OK - make sure you add in fortran after the last member currently accessed in python, e.g. after wGauss2 in TRecfast. Parameters must in the same order in fortran and python and sequential (but not all fortran type parameters have to be accessible from python if they are defined at the end after the last one defined in python).

Joshua Benabou
Posts: 5
Joined: September 14 2022
Affiliation: University of California Berkeley

Re: CAMB: passing a new CAMBparams parameter to recfast.f90

Post by Joshua Benabou » September 14 2022

Dear Anthony,

Right, that is the setup I have currently. The issue is, running camb.get_results(pars) works fine if I don't do anything else and dont modify pars from what camb.CAMBparams spits out, and I see that the recombination code uses the value of newparam initialized in recfast.f90 (real(dl) :: newparam = 1D1 )

However accessing pars.Recomb shows that newparam has initial value zero. So Im not sure how Python initializes this.

When I then try to modify in my Jupyter notebook pars.Recomb.newparam = 1 (or any other value) and rerun camb.get_results(pars), my notebook crashes, but no error is given so Im not sure what is going on.

I also see that first changing pars.Recomb.newparam = 1, then running camb.get_results(pars) ( so only running this once) also crashes.

Further, running camb.get_results(pars) twice without changing pars, also crashes.

thanks again.

Joshua Benabou
Posts: 5
Joined: September 14 2022
Affiliation: University of California Berkeley

Re: CAMB: passing a new CAMBparams parameter to recfast.f90

Post by Joshua Benabou » September 15 2022

Dear Anthony,

Here is a minimal example demonstrating the problem:

Without any source code changes, running the following code block twice in a Jupyter notebook works fine:

Code: Select all

pars = camb.CAMBparams()
pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122, mnu=0.06, omk=0, tau=0.06)
pars.InitPower.set_params(As=2e-9, ns=0.965, r=0)
pars.set_for_lmax(2500, lens_potential_accuracy=0);
results = camb.get_results(pars)
print("Done!")
If I just add the following line right after line 267 in recfast.f90 (i.e right after where 'wGauss2' is defined), then running the above code block once is fine, but twice causes my Jupyter to crash.

Code: Select all

real(dl) :: newparam = 1D-22
Running this in a python scipt, I get Segmentation fault.

Thanks.

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

Re: CAMB: passing a new CAMBparams parameter to recfast.f90

Post by Antony Lewis » September 15 2022

Not sure, crashes are usually having things in the wrong order, wrong type, or not rebuilding the binary consistently/not re-starting notebook kernels.

Joshua Benabou
Posts: 5
Joined: September 14 2022
Affiliation: University of California Berkeley

Re: CAMB: passing a new CAMBparams parameter to recfast.f90

Post by Joshua Benabou » September 15 2022

Dear Anthony,

I have followed your instructions verbatim concerning the ordering (I added the new parameter definition immediately after the last variable which is also defined in the Python wrapper). The type is the same as 'wGauss2' so this should not be a problem.

I am running

python setup.py make
pip install -e .

everytime after making changes, and then restarting the Jupyter kernel.

Could this segfault issue be fortran - version dependent? Any tips on how to debug this?

thanks again.

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

Re: CAMB: passing a new CAMBparams parameter to recfast.f90

Post by Antony Lewis » September 16 2022

Not very easy debugging binary-python interface if it works OK from pure fortran. Debug fortran build may help. Can also check camb is loading from where you think it is, and maybe trying python setup.py clean before rebuilding. Seems unlikely to be a version issue.

Post Reply