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

test: reproduce Pober+2014 results #82

Merged
merged 22 commits into from
Aug 1, 2023
Merged

Conversation

steven-murray
Copy link
Contributor

Adds a new tutorial notebook that tries to reproduce (some of) the results of https://arxiv.org/pdf/1310.7031.pdf. It also adds a few cells that run explicit tests so that our tests will fail if something is not right.

While the plots look roughly similar to Figure 5, it seems the new vesion has a deficit of SNR overall. I'm not sure why this is, but it can be up to a ~35% difference, and gets larger the more small-k modes are being used (i.e. it's larger for optimistic than moderate than pessimistic). There are a number of "small" changes in the current vesion to do with accuracy etc, so I think this version probably does a better job than the old, rather than it being a bug. But I can't be quite sure.

I would go further to investigate but it's super duper hard to run the old code now since it's Python 2.7.

@review-notebook-app
Copy link

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

@codecov
Copy link

codecov bot commented Oct 4, 2022

Codecov Report

Merging #82 (a1d232c) into fix-issue-80 (31c2ed2) will decrease coverage by 0.08%.
The diff coverage is 100.00%.

@@               Coverage Diff                @@
##           fix-issue-80      #82      +/-   ##
================================================
- Coverage         98.42%   98.34%   -0.08%     
================================================
  Files                13       14       +1     
  Lines               951      967      +16     
================================================
+ Hits                936      951      +15     
- Misses               15       16       +1     
Files Changed Coverage Δ
py21cmsense/config.py 100.00% <ø> (ø)
py21cmsense/__init__.py 100.00% <100.00%> (ø)
py21cmsense/_utils.py 100.00% <100.00%> (ø)
py21cmsense/antpos.py 100.00% <100.00%> (ø)
py21cmsense/beam.py 98.07% <100.00%> (+0.11%) ⬆️
py21cmsense/conversions.py 100.00% <100.00%> (ø)
py21cmsense/data/__init__.py 100.00% <100.00%> (ø)
py21cmsense/observation.py 95.62% <100.00%> (-1.16%) ⬇️
py21cmsense/observatory.py 98.24% <100.00%> (+0.01%) ⬆️
py21cmsense/sensitivity.py 97.90% <100.00%> (+0.37%) ⬆️
... and 2 more

@steven-murray steven-murray self-assigned this Oct 4, 2022
@@ -0,0 +1,721 @@
{
Copy link
Contributor

Choose a reason for hiding this comment

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

It does seem good. For a check, the beam was defined as a Gaussian with a sigma of sigma = 0.0643 at 150 MHz


Reply via ReviewNB

@@ -0,0 +1,721 @@
{
Copy link
Contributor

Choose a reason for hiding this comment

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

This isn't quite right. The code at the time used aipy cal files, which contained a longitude and latitude. For this paper, the array was actually in Green Bank!

'loc': ('38:25:59.24', '-79:51:02.1'), # Green Bank, WV


Reply via ReviewNB

@steven-murray
Copy link
Contributor Author

@jpober I finally got round to updating this. However, it seemed to push our results further away from each other (still around 30% difference in SNR). Still, the shape of the SNR and the noise power is the same in the new code as your original.

It would be nice to be able to run your code to check things out, but I'm afraid I lost that ability a while ago.

@jpober
Copy link
Contributor

jpober commented Nov 4, 2022

Oh, darn, that's unfortunate. There's a chance I may have an old laptop where my code is just ready to go. I'm traveling right now, but I'll look into it and see when I'm back.

@jpober
Copy link
Contributor

jpober commented Nov 17, 2022

Unfortunately it doesn't look like I can run 21cmSense off my old machine either. I can try to set up an Anaconda environment that can run my old repo, but it's going to be a little while longer yet...

@jpober
Copy link
Contributor

jpober commented Jan 17, 2023

I spent some time the other day trying to get a working environment to run my old code, but installing aipy is just too painful. Should we maybe have a call and figure out how we want to close this loop?

@steven-murray
Copy link
Contributor Author

@jpober I think this is ready to have one more look over. I updated the "reproduce_pober2014" notebook after your fix, and I'm now finding that the new code is always slightly more sensitive than the plots in the paper, but always less than 15%, and roughly equally for opt/pess/moderate.

I'm happy with this, but I would like to verbally explain in the notebook what the remaining differences are. I know there is the cosmological calculations, but is there anything else that you had to tweak manually? You mentioned at some point that the obs duration had to be hand-set to 35 minutes, whereas in the notebook I calculate it to be 38 minutes. This is a ~10% effect if that's true.

@jpober
Copy link
Contributor

jpober commented Jul 14, 2023

Here is a complete list of changes I had to make to the test-pober-2014 branch get the results shown in #70 (in addition to the actual bug fix).

-In the observation yaml:

  1. Change obs_duration to 2077.38 (from 10.7)
  2. Change n_channels to 82 (from 60)

-In the observatory yaml:

  1. Change the hex_num to 11 (from 3)
  2. Change the beam frequency to 135 (from 150)
  3. Change the latitude to Green Bank, i.e. 0.6707845 (from 0.536189)

-In the sensitivity yaml:

  1. Change the theory model to Legacy21cmFAST (from EOS2021)
    -In antpos.py:
  2. Revert back to the version in main and hardcore the antenna separation to 12.12 m (instead of sin(60 degrees)*separation)

-In conversions.py:

  1. Replace the correct cosmological distance calculations in dL_dth and dL_df with the approximate relations from Furlanetto et al. 2006.

-In observation.py:

  1. Change the system temperature calculation by (i) changing the spectral_index to 2.6 (from 2.55); (ii) changing tsky_amplitude to 60000 (from 260000); and (iii) changing tsky_ref_freq to 300 (from 150)

That's a pretty long list, but I think more than accounts for any remaining discrepancies.

@steven-murray
Copy link
Contributor Author

@jpober thank you, that is very helpful! Those changes were to reproduce your memo which used a slightly different version of the code (and some different parameters) than the paper, right? I would just switch the demo notebook to reproduce the memo, but it seems more useful to reproduce the paper, since it has more plots -- the memo doesn't have much data to compare to (on the surface). We could, on the other hand, reproduce the memo results if we just included some test data that you produce with your old version of the code (i.e. the data that went into the plots in #70), and have a statement in the demo that what we're comparing to is data from a memo, but that corresponds closely to the paper. What do you think?

@jpober
Copy link
Contributor

jpober commented Jul 14, 2023

I think comparing with the memo is probably better. 21cmSense didn't actually exist when I wrote the 2014 paper, just the code that eventually became 21cmSense. I can generate some output files with the memo-matching code. What in particular are you imagining? Anything beyond the noise power spectrum?

@steven-murray
Copy link
Contributor Author

I think the noise power (or noise + sample variance, if that's easy enough) is good. If I have that data, I can more easily make plots that showcase how close (or otherwise) the results are. Thanks!

@steven-murray steven-murray merged commit ff23b13 into fix-issue-80 Aug 1, 2023
18 of 19 checks passed
@steven-murray steven-murray deleted the test-pober-2014 branch August 1, 2023 22:00
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.

2 participants