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

Bug in two_point_clustering.tpcf_multipole #651

Closed
johannesulf opened this Issue Jan 28, 2017 · 3 comments

Comments

Projects
None yet
3 participants
@johannesulf
Copy link
Contributor

johannesulf commented Jan 28, 2017

@duncandc As mentioned in #649, based on some science tests I believe that the multipoles calculated by halotools are a factor of 2 too low. Mathematically, the multipoles are the coefficients for expanding xi(s, mu) into a series of Legendre polynomials P_l (mu), where l is the order. Take the following code.

import numpy as np
from halotools.mock_observables import tpcf_multipole

xi_mu_s_fake = np.ones((49, 49))
mu_bins = np.linspace(0, 1, 50)
print tpcf_multipole(xi_mu_s_fake, mu_bins, order=0)

We calculate the monopole (order = 0) of a 2PCF that is 1 at all radii and LOS angles. Clearly, the monopole should be 1 at all radii because the 0th order Legendre polynomial is also 1. However, I get an array of 0.5. A factor of 2 too low. Let's look at the source code of tpcf_multipole.

# process inputs
s_mu_tcpf_result = np.atleast_1d(s_mu_tcpf_result)
mu_bins = np.atleast_1d(mu_bins)
order = int(order)

# calculate the center of each mu bin
mu_bin_centers = (mu_bins[:-1]+mu_bins[1:])/(2.0)

# get the Legendre polynomial of the desired order.
Ln = legendre(order)

 # numerically integrate over mu
 result = (2.0*order + 1.0)/2.0 *\
        np.sum(s_mu_tcpf_result * np.diff(mu_bins) * Ln(mu_bin_centers), axis=1)

This code calculates the coefficients/multipoles by integrating over mu. The problem is that mu is only defined in [0, 1] in halotools. The source code of tpcf_multipole explicitly states that and s_mu_tpcf also only accepts [0, 1]. However, the integration should be going from -1 to 1 in mu, not 0 to 1. That means all multipoles of even order (monopole, quadropole, hexadecapole...) where the Legendre polynomial is even are currently underestimated by a factor of 2. On the other hand, all multipoles of odd order now give large non-zero values. However, they should be zero because their Legendre polynomials are odd and the 2PCF even as a function of mu. Does that make sense?

An easy fix would be to change the last line to

 # numerically integrate over mu
 result = (2.0*order + 1.0)/2.0 *\
        np.sum(s_mu_tcpf_result * np.diff(mu_bins) * (Ln(mu_bin_centers) + Ln(-mu_bin_centers)), axis=1)

@aphearin aphearin added this to the v0.5 milestone Feb 1, 2017

@aphearin aphearin self-assigned this Feb 1, 2017

@aphearin aphearin assigned duncandc and unassigned aphearin Feb 26, 2017

@duncandc

This comment has been minimized.

Copy link
Contributor

duncandc commented Mar 9, 2017

I think @johannesulf's fix here is exactly the right thing to do. I can implement this, or Johannes can since its his good eye...

@aphearin

@duncandc

This comment has been minimized.

Copy link
Contributor

duncandc commented Mar 9, 2017

This is also a good test. The odd multipoles should be zero to within numerical precision.

johannesulf added a commit to johannesulf/halotools that referenced this issue Mar 9, 2017

aphearin added a commit that referenced this issue Mar 10, 2017

Merge pull request #698 from johannesulf/tpcf_multipole_dev
Fixed bug in tpcf_multipole.py, addressing #651
@aphearin

This comment has been minimized.

Copy link
Contributor

aphearin commented Mar 22, 2017

After some discussion with @duncandc, we decided to include a brute force comparison of npairs_s_mu to a pure python implementation in order to consider this resolved. See #711.

At some point, it would be useful to find someone who has an independently written multipole calculation that uses the distant observer approximation, but since we do not know of such a code, we'll stick with improving npairs_s_mu testing for the time being.

@aphearin aphearin closed this Mar 22, 2017

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment