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

Power spectrum responses for SSC #1134

Open
wants to merge 25 commits into
base: master
Choose a base branch
from
Open

Power spectrum responses for SSC #1134

wants to merge 25 commits into from

Conversation

RyoTerasawa
Copy link

In pyccl/pkresponse.py, I implemented the functions to calculate power spectrum responses to the super-survey modes which are responsible for super-sample covariance (SSC), based on arXiv:2310.13330.
We approximate the power spectrum growth response to the super-survey modes by its growth response to the Hubble parameter.

For the galaxy-matter and galaxy-auto power spectra, the halo statistics (halo-matter/halo-auto power spectrum, mass function) are calculated by DarkEmulator.

@coveralls
Copy link

coveralls commented Dec 4, 2023

Pull Request Test Coverage Report for Build 12132306089

Details

  • 304 of 346 (87.86%) changed or added relevant lines in 2 files are covered.
  • No unchanged relevant lines lost coverage.
  • Overall coverage decreased (-0.5%) to 96.946%

Changes Missing Coverage Covered Lines Changed/Added Lines %
pyccl/pkresponse.py 303 345 87.83%
Totals Coverage Status
Change from base Build 11958349869: -0.5%
Covered Lines: 6824
Relevant Lines: 7039

💛 - Coveralls

@damonge
Copy link
Collaborator

damonge commented Feb 21, 2024

@RyoTerasawa can I check what the status of this PR is?

@RyoTerasawa
Copy link
Author

I am waiting for someone to review the code. I haven't asked for a specific person yet, but maybe I should ask a member of the MCPCov group to review it. I am happy to explain the details if needed. I am also working on writing the unit test and benchmark test to be included in this PR.

@carlosggarcia
Copy link
Contributor

@RyoTerasawa @YueNan-c, maybe you want to also implement the DarkEmulator here: https://github.com/LSSTDESC/CCL/blob/master/pyccl/emulators/cosmicemu_pk.py.

Are the unit tests and benchmark ready?

@damonge
Copy link
Collaborator

damonge commented Jun 18, 2024

@RyoTerasawa @carlosggarcia can I check what the status of this is?

@RyoTerasawa
Copy link
Author

@carlosggarcia @damonge
Sorry for the late reply. I just added the unit tests and benchmark, and this pull request is ready to be reviewed.

Copy link
Contributor

@carlosggarcia carlosggarcia left a comment

Choose a reason for hiding this comment

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

Ok. I first round of comments. Some suggestions:

  • Add unit tests per function/method declared to be sure they are doing what they should
  • I would unify the .npy files in a single .npz file. Also, it seems that you are only using the z=0 case, I'd be great if you can compare all redshift, since you have it.
  • Add References to the equations in your paper for later reference (and to make my life easier when reviewing the PR)
  • @damonge, there are integrals and biases computed that I think they might be obtained from the current CCL implementation. Can you have a look? The are in pkresponse.py, where I pinged you, too.

In general, it looks good, though! Good work!

def test_Pmm_resp():
response = Pmm_resp(cosmo, deltah=deltah, lk_arr=lk_arr, a_arr=a_arr)

assert np.all(np.isfinite(response)), "Pmm_resp produced infinity values."
Copy link
Contributor

Choose a reason for hiding this comment

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

What is this "Pmm_resp produced infinity values." doing?

Copy link
Author

Choose a reason for hiding this comment

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

I realized the error message is not needed, so I simply deleted the message.

def Pmm_resp(
cosmo,
deltah=0.02,
extra_parameters={
Copy link
Contributor

Choose a reason for hiding this comment

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

These aren't used anywhere.

Copy link
Author

Choose a reason for hiding this comment

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

I added dependency to extra_parameters.



def set_hmodified_cosmology(cosmo, deltah):
Omega_c = cosmo["Omega_c"]
Copy link
Contributor

Choose a reason for hiding this comment

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

hmm, and what if there are extra parameters given, eg. massive neutrinos? Wouldn't that be relevant? I guess you should make sure the cosmology object is exactly the same, except for the h

Copy link
Author

Choose a reason for hiding this comment

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

In the massive neutrinos' case, m_nu should be the same. I edited the function set_hmodified_cosmology to explicitly make sure the cosmology object is exactly same except for the h, by starting with the copy of the fiducial cosmology object through its keys.

return b2


def darkemu_set_cosmology(emu, cosmo):
Copy link
Contributor

Choose a reason for hiding this comment

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

Add and underscore if these methods are meant to be private (so not used by other libraries); e.g. _darkemu_set_cosmology. Also, add documentation. Same with the other "utility functions".

Copy link
Author

Choose a reason for hiding this comment

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

I made all the utility functions private. Also added documentation.

Copy link
Contributor

Choose a reason for hiding this comment

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

In general, you should add a unit test for each function that you have defined in pkresponse.py to check that they are working properly. For instance, for set_hmodified_cosmology you would like to check that the output cosmologies are exactly the same as the input one, except for the h's.

Copy link
Author

Choose a reason for hiding this comment

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

I added a unit test for all functions defined in pkresponse.py. For set_hmodified_cosmology the unit test checks that the output cosmologies are exactly the same as the input one, except for the h's.


kmin = 1e-2
for ia, aa in enumerate(a_arr):
pk = pk2d.__call__(k_use, aa, cosmo)
Copy link
Contributor

Choose a reason for hiding this comment

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

You don't need to use __call__, you can just call it, pk2d(k_use, aa, cosmo)

Copy link
Author

Choose a reason for hiding this comment

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

Replaced the all __call__ with pk2d() in pkresponse.py.

hbf = halos.HaloBiasTinker10(mass_def=mass_def)

# dark emulator is valid for 0 =< z <= 1.48
if np.any(1.0 / a_arr - 1) > 1.5:
Copy link
Contributor

Choose a reason for hiding this comment

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

if the maximum is 1.48, check for 1.48, not 1.5

Copy link
Author

Choose a reason for hiding this comment

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

Modified the pkresponse.py to check for 1.48.


# dark emulator is valid for 0 =< z <= 1.48
if np.any(1.0 / a_arr - 1) > 1.5:
print("dark emulator is valid for z={:.2f}<1.48")
Copy link
Contributor

Choose a reason for hiding this comment

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

hmm, what does {:.2f} do here? I think you have to remove that.

Copy link
Author

Choose a reason for hiding this comment

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

The {:.2f} was unnecessary. I removed that.

/ ng
)

bgE2 = (
Copy link
Contributor

Choose a reason for hiding this comment

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

@damonge , can this be done with the current methods of the halo model?

Copy link
Contributor

Choose a reason for hiding this comment

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

Same with the other integrals

Copy link
Collaborator

Choose a reason for hiding this comment

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

I think so. @RyoTerasawa can you check:

def _integrate_over_mf(self, array_2):

and
def _integrate_over_mbf(self, array_2):

and
def integrate_over_massfunc(self, func, cosmo, a):

?

Copy link
Author

Choose a reason for hiding this comment

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

I changed the pkresponse.py, and now all the integrals of the mass function are using _integrate_over_mf or _integrate_over_mbf.

# Construct the full path for each data file (z=0)
k_data_path = os.path.join(data_directory_path, "k_h.npy")
k_data_mm_path = os.path.join(data_directory_path, "k_h_mm.npy")
Pmm_resp_data_path = os.path.join(data_directory_path, "Pmm_resp_z0.npy")
Copy link
Contributor

Choose a reason for hiding this comment

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

You have so much more data generated, at different redshifts. Why don't you use it?

Copy link
Author

Choose a reason for hiding this comment

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

I added the tests for higher redshifts, z ~ [0.5, 1.0, 1.5].

@damonge
Copy link
Collaborator

damonge commented Nov 13, 2024

@carlosggarcia @RyoTerasawa can I check what the status of this PR is?

@RyoTerasawa
Copy link
Author

@carlosggarcia @RyoTerasawa can I check what the status of this PR is?

I am currently working on addressing Carlos's comments.

@RyoTerasawa
Copy link
Author

@carlosggarcia Thank you for the comments and suggestions.
I complete addressing these comments except for the integrals.

  1. Added unit tests per function/method declared to be sure they are doing what they should
  2. Unified the .npy files in a single .npz file and now the benchmark test compares all redshift.
  3. Added References to the equations in the paper.
  4. Have not addressed to the integrals in pkresponse.py

@damonge, there are integrals and biases computed that I think they might be obtained from the current CCL implementation. Can you have a look? The are in pkresponse.py, where I pinged you, too.

@damonge
Copy link
Collaborator

damonge commented Dec 4, 2024

@RyoTerasawa I responded to Carlos' comment about the integrals. If possible, it'd be good if CCL does all of these consistently.

@RyoTerasawa
Copy link
Author

@damonge Thank you for the suggestion. I switched to the integral method of pyccl/halos/halo_model.py and it helped to make the code bit simpler.

Merge branch 'master' into SSC
@RyoTerasawa
Copy link
Author

@carlosggarcia, I complete addressing these comments. Although the build failed at the unit test.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants