manodeep / Corrfunc

⚡️⚡️⚡️Blazing fast correlation functions on the CPU.
https://corrfunc.readthedocs.io
MIT License
166 stars 51 forks source link

w(theta) brute-force test failure #107

Closed manodeep closed 7 years ago

manodeep commented 7 years ago

Most likely the failures are on linux - have not seen any OSX failures yet. The raw pair counts for the DR brute-force tests are off by a few counts.

Status OS Compilers
Fail Ubuntu 16.04.4 gcc 4.8, 4.9, 5.4
Fail CentOS 6.5 icc 15
Pass CentOS 6.8 gcc 4.6.2, 4.9.3, 5.1.0
Pass CentOS 6.8 gcc 5.2.0
Pass OSX 10.11.6 (El Capitan) gcc 4.8.5, 5.4.0, 6.3.0, clang 3.8.1, 3.9.1, 4.0.0

@lgarrison

lgarrison commented 7 years ago

Appears to be fixed in 5eda1511147369f0fb33947408fea0961045618f. The "developer mode" tests introduced there should detect any regressions. We may still want to investigate what exactly caused the failures in the first place.

manodeep commented 7 years ago

Since the only relevant change in https://github.com/manodeep/Corrfunc/commit/5eda1511147369f0fb33947408fea0961045618f was to fix the weights typo and that change was also present in the master branch you were testing, I don't understand how https://github.com/manodeep/Corrfunc/commit/5eda1511147369f0fb33947408fea0961045618f could fix this issue!

That being said, I am not sure how to proceed either.

lgarrison commented 7 years ago

It appears the change of

const DOUBLE dz = *z1 - zpos;

to

const DOUBLE dz = (*z1) * zpos; 

in mocks/wtheta/countpairs_theta_mocks_kernels.c.src is the reason the new commit passes. Was this a bugfix?

manodeep commented 7 years ago

Of course! I forgot that bruteforce calls the fallback kernel.

The previous construct was incorrect but the current one isn't quite correct either. I am now pretty sure that that dz equation is why the PR #108 is failing.

Here's the math: z = sin(dec)

We need to find the first z1 such that |dec1 - decpos| < thetamax; where the following conditions hold true:

1.  -1 <= x1, y1, z1, xpos, ypos, zpos <= 1
2. x1^2 + y1^2 + z1^2 = 1; xpos^2 + ypos^2 + zpos^2 = 1
3. zpos and z1 are sorted
4. x1^2 + y1^2 = cos^2(dec1); xpos^2 + ypos^2 = cos^2(decpos)

The original implementation was:

dz = sin(dec1) - sin(decpos)
=> dz = 2 sin((dec1 - decpos)/2.0) * cos ((dec1 + decpos)/2.0)
=> dz <= 2 sin((dec1 - decpos)/2.0) since the angle argument in the second term satisfies the inequality 0 <= 0.5 * (dec1 + decpos) <= 90.0, the second term itself must be >=0 and <= 1.

Which wasn't really giving a constraint on |dec1 - decpos|.

The current one is directly computing one term in the cos(theta) expression for the angular separation and iff all terms have the same sign (i.e., sin(ra1) * sin(rapos) >= 0), then the inequality works out.

Let me see if I can work out the necessary math for the general case. Or we can simply take out that entire section and remove the sorting on z within gridlink_mocks_theta_ra_dec_DOUBLE and gridlink_mocks_theta_dec_DOUBLE functions.

manodeep commented 7 years ago

I am happy to close this issue for now and focus on fixing PR #108

@lgarrison - what do you think?

lgarrison commented 7 years ago

Agreed!