search

rmjarvis / TreeCorr

Code for efficiently computing 2-point and 3-point correlation functions. For documentation, go to
http://rmjarvis.github.io/TreeCorr/
Other
97 stars 37 forks source link

multiply-occurring objects bias results low #134

Closed SandraUnruh closed 2 years ago

SandraUnruh commented 2 years ago

I got unexpected shear-shear correlation results with a catalog including galaxies that can appear more than once. To test it I just compared results from a "normal" catalog of galaxies with one where I put in every galaxy twice and I got the following example result:

check_treecorr

The signal of the doubled catalog is consistently suppressed by ~30% for low theta while for large theta they roughly agree. Unfortunately, I could not really pinpoint where the calculation goes wrong. It is, e.g., not just the number counts per bin that are wrong.

My current work-around is to shift the position of every single duplicated galaxy by <0.5 arc seconds in a random way. Then the original and duplicated catalog agree perfectly.

I thought, however, you might find an easy fix for this or at least mention it in the documentation (maybe I just missed it though ;)

Cheers, Sandra Unruh

rmjarvis commented 2 years ago

Hi Sandra. I couldn't reproduce this with a straightforward test that duplicated all the galaxies. I'm getting answers that are the same to machine precision, so there must be something in your script that I'm not including in my test.

Can you please send me a self-contained script that shows the effect? Either post here, or send me by email or Slack.

SandraUnruh commented 2 years ago

Dear Mike,

I put together a minimal script that reproduces my error. It can be found here

https://www.dropbox.com/sh/8ko0en0xz7edp89/AABJKY_IpdhcXsTzyRwVF4hQa?dl=0 https://www.dropbox.com/sh/8ko0en0xz7edp89/AABJKY_IpdhcXsTzyRwVF4hQa?dl=0

with a 26MB catalog of input data (which shouldn’t matter but still). Please change the data directory (data_dir) to yours, the it should run through.

I use these versions of Python and TreeCorr: TreeCorr.version = ‚4.2.3' Python 3.8.10 | packaged by conda-forge

I really hope it is no stupid mistake from my side! Thank you for your time, I enjoy working with TreeCorr a lot. Cheers, Sandra

Am 23.09.2021 um 18:31 schrieb Mike Jarvis @.***>:

Hi Sandra. I couldn't reproduce this with a straightforward test that duplicated all the galaxies. I'm getting answers that are the same to machine precision, so there must be something in your script that I'm not including in my test.

Can you please send me a self-contained script that shows the effect? Either post here, or send me by email or Slack.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/rmjarvis/TreeCorr/issues/134#issuecomment-925972766, or unsubscribe https://github.com/notifications/unsubscribe-auth/AKUJROO3TJPSDHRKNGI7O43UDNI7TANCNFSM5EUAQKFA. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

rmjarvis commented 2 years ago

Thanks Sandra. That was perfect. This was indeed a bug in TreeCorr, not in your code.

The root problem was in the step of parallel transporting shears to the center of a cell to take the average. When all the shears were at the same place, some calculations could sometimes be numerically very tiny but not zero, when perfect math would have them be 0.0. This then shows up in a denominator. The numerator is also essentially zero, so you get 0/0, which meant the phase factor exp(i phi) was some random value, not equal to 1, which is what you should have when you don't actually need to parallel transport at all (since all the points are at the center of the cell).

I had a check for exactly 0.0 to skip the phase factor, which worked when there was only one object in a cell. But in your case with the duplicated objects, rounding errors prevented it from being exact. Essentially sin(phi) was calculated as something of order 3.e-19/1.e-17 (it varies for each cell -- probably some were worse than this), which was not exactly 0, and so the shears were rotated a little bit. The rotation was random, which had the systematic effect of diluting the signal, thus explaining the decrease you saw for the duplicated catalog.

Anyway, this is probably more information than you care about, but suffice to say, I very much appreciate the bug report! I'll start the process now of merging this and cutting a new release with the fix. Should be out sometime tomorrow.

rmjarvis commented 2 years ago

v4.2.6 is released with this fix. Thanks again Sandra.