Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Halo model for power spectrum calculation #339

Merged
merged 152 commits into from Aug 21, 2018
Merged

Halo model for power spectrum calculation #339

merged 152 commits into from Aug 21, 2018

Conversation

villarrealas
Copy link
Collaborator

This branch works to add a halo model calculation and refactor some of the code in ccl_massfunc as a result of this. This will allow a new method of power spectrum calculation.

I think one key thing prior to pushing to master will be to develop benchmarks of this for validation purposes.

Copy link
Collaborator

@philbull philbull left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Minor changes needed to tidy this up (e.g. removing comments, fixing other comments), and maybe a few things to rename to avoid namespace clashes. But after that, this is good to go.


double gz, g0, nu, delta_c, a_form;
double Mpiv, A, B, C;

switch(label){

// something something crazy Bullock
// Bullock et al. (2001; xxxx.xxxx; Delta = Virial)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This needs to be tidied up.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

@@ -181,7 +184,7 @@ static double window_function(ccl_cosmology *cosmo, double m, double k, double a
rho_matter = ccl_rho_x(cosmo, 1., ccl_species_m_label, 1, status);

// The halo concentration for this mass and scale factor
c = halo_concentration(cosmo, m, a, odelta, Duffy2008_virial, status);
c = ccl_halo_concentration(cosmo, m, a, odelta, Duffy2008_virial, status);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be fixed to the Duffy version?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, I did a fairly substantial re-write of things to change this. The halo_concentration is now elevated to the config for the cosmology object, putting on the same status as mass_function. The user can now choose combinations of mass function and halo concentration that they would like to use for the halo-model power-spectrum calculation. The only snag is that the halo-model calculation is fixed to assume that haloes are defined by the virial definition, rather than say M200 or M500. The calculation enforces consistency though, so you cannot use a mass function for haloes with M500 in tandem with a concentration relation for haloes defined by M200. For the power spectrum calculation you have to use relations that are defined for virial haloes. I think this is now done!

}

/*----- ROUTINE: nu_mass -----
//TODO: Remove this function if no longer useful
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be removed?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed; done

*/
double nu_mass(ccl_cosmology *cosmo, double halomass, double a, int *status);
*/
// double nu_mass(ccl_cosmology *cosmo, double halomass, double a, int *status);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove this if no longer needed.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

.gitignore Outdated
@@ -61,6 +70,16 @@ pyccl.egg*
# Ignore VCS files
.idea/*



# Ignore document pdf and bibliography
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will probably break the ccl_paper branch...

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, I removed it. Done.

// Equation (12) in arXiv: 9901122
// Derived using the peak-background split applied to the mass function in the same paper
// Note that Sheth & Tormen (1999) use nu=(dc/sigma)^2 whereas we use nu=dc/sigma
// TODO: Mead: Somehow enforce that this should only work with a virial-density Delta_v
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

TODO: Does this need to be addressed?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It had already been addressed. Only the comment remained. Done


/*----- ROUTINE: ccl_massfunc -----
INPUT: ccl_cosmology * cosmo, double halo mass in units of Msun, double scale factor
TASK: returns halo mass function as dn / dlog10 m
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need to mention units, e.g. is it in comoving or physical Mpc^-3, is it in h^-1 units etc.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

src/ccl_power.c Outdated
}

// TODO: log10 could be taken already in the macros.
// TODO: 1E-5 should be a macro
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Address these TODO items.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These comments were inherited from the ccl_sigmaR function, from which the ccl_sigmaV function was copied. I don't know what they refer to: I don't know what log10 could be taken, what relevance 1E-5 is, and it looks to me like gsl integration success is already checked for. Maybe @villarrealas knows something more... ?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not entirely sure what these TODO items are for in sigmaR (I forget who wrote ccl_power off the top of my head), but I believe both of these look like they are already resolved.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What about // TODO: we should check for integration success ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Apparently this is sorted, so I've now done

# halo concentration
#assert_( all_finite(ccl.halo_concentration(cosmo, mass_scalar, a)) )
#assert_( all_finite(ccl.halo_concentration(cosmo, mass_list, a)) )
#assert_( all_finite(ccl.halo_concentration(cosmo, mass_array, a)) )
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove commented code.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

@@ -732,8 +784,8 @@ def test_valid_transfer_combos():

assert_raises(ValueError, ccl.Cosmology, transfer_function='emulator',
matter_power_spectrum='linear', **params)
assert_raises(ValueError, ccl.Cosmology, transfer_function='boltzmann',
matter_power_spectrum='halomodel', **params)
#assert_raises(ValueError, ccl.Cosmology, transfer_function='boltzmann',
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why commented out?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, a bit of a long story here:

First; previously there were halo_model and halomodel cases for ccl_nonlin_matter_power, neither of which were wired up to anything. If you called ccl_nonlin_matter_power with each method I think you just got the linear power back. I removed one of them so we are now left with halo_model only. I couldn't think of any good reason to have both.

Second; I looked at wiring the halo_model option up to the new calculation from this new halomod branch, but this is impossible I think because ccl_halomod.c uses ccl_power.h, so I can't use the correct routines without entering some sort of infinite calling loop. So at the moment the halo_model option remains useless.

This sort of relates to one of the issues I raised previously (#337) and also sort of to #438

Generally, I think it would be better to have a source file dedicated to getting the linear power spectrum (ccl_linear_power.c ?) and then a separate file(s) dedicated to whatever non-linear method you wanted to use... and then one final source file that had a way of choosing between whatever power spectrum calculation the user needed, be that linear, emu, halomodel, halofit... whatever.

@philbull
Copy link
Collaborator

philbull commented Aug 8, 2018

Also, I see that we have a test failing due to a precision issue:

TEST 58/65 halomod:model_2 [FAIL]
  ERR: /home/travis/build/LSSTDESC/CCL/tests/ccl_test_halomod.c:156  expected 1.388688e-01, got 1.390115e-01 (diff -1.427826e-04, tol 1.388688e-04)

@philbull
Copy link
Collaborator

philbull commented Aug 8, 2018

@elisachisari suggests that this might be caused by the recent changes to the default set of physical constants.

@elisachisari
Copy link
Collaborator

elisachisari commented Aug 8, 2018

Here is a concrete set of tasks that need to be addressed concerning this PR.

  • Address review by @philbull
  • Fix failure in CTEST, most likely due to the new physical constants adopted by CCL @damonge
  • Upload benchmark code that produces the benchmark to the relevant branch @damonge
  • Make a plot for the paper comparing to the benchmarks - add this in the paper and in the jupypter notebook Jupyter notebook with benchmark plots for paper #369. @tilmantroester
  • Address comments in the paper itself: making sure notation is homogeneous with respect to HMF section, add implementation and validation details.

@tilmantroester
Copy link
Contributor

Benchmark plot is now in the notebook and in the paper.

Copy link
Collaborator

@philbull philbull left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is ready to go!

@elisachisari elisachisari merged commit 73f6413 into master Aug 21, 2018
@beckermr beckermr deleted the halo_model branch December 12, 2018 13:03
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

9 participants