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

Pairs at μ=1.0 #194

Open
mehdirezaie opened this issue Sep 5, 2019 · 18 comments
Open

Pairs at μ=1.0 #194

mehdirezaie opened this issue Sep 5, 2019 · 18 comments
Assignees

Comments

@mehdirezaie
Copy link

Hi! I was wondering how one could have the DDsmu include the pairs that have exactly μ = 1.0. Any help or suggestion are greatly appreciated.

@mehdirezaie mehdirezaie changed the title Pairs at cos(theta)=1.0 Pairs at cos(separation angle)=1.0 Sep 5, 2019
@mehdirezaie mehdirezaie changed the title Pairs at cos(separation angle)=1.0 Pairs at μ=1.0 Sep 5, 2019
@manodeep manodeep self-assigned this Sep 5, 2019
@manodeep
Copy link
Owner

manodeep commented Sep 5, 2019

Hi - thanks for opening the issue and considering Corrfunc.

By design Corrfunc always uses bin comparisons using bin_low <= r < bin_hi and changing that will require you to change a few lines of code within the kernel. For example, in the AVX512F kernel, you will have to change these lines replacing a >= with a > (or an equivalent action)

  • Change to (zpos - last_z1) > max_all_dzhere (Even though it looks like that code is commented out)
  • Change to > instead of >= here
  • Change to min_dz > pimax here
  • Change to *localz1 < target_z here
  • Change to AVX512_MASK_COMPARE_FLOATS(m_mask_left, m_zdiff, m_max_dz,_CMP_GT_OQ); here
  • Change to AVX512_MASK_COMPARE_FLOATS(m_mask_left, m_sqr_zdiff, max_sqr_dz, _CMP_LE_OQ); here
  • You will need an additional check here to make sure that the mu=mu_max bin is added to a valid bin. One solution would be the following:
const AVX512_MASK m_mu_eq_mumax_mask = AVX512_MASK_COMPARE_FLOATS(m_mask, m_sqr_mu, m_sqr_mumax, _CMP_EQ_OQ); //see which mu^2 values are equal to mu_max^2
AVX512_FLOATS m_mubin = AVX512_MASKZ_MULTIPLY_FLOATS(m_mask, m_mu, m_inv_dmu); //perform the normal multiplication as before to calculate the mu bin
const AVX512_FLOATS m_nmu_bins_m1 = AVX512_SET_FLOAT((DOUBLE) (nmu_bins - 1));//the maximum valid mu bin
m_mubin = AVX512_BLEND_FLOATS_WITH_MASK(m_mu_eq_mumax_mask, m_mubin, m_nmu_bins_m1);//set all the mu_max values to fall into the final valid mu-bin

Totally untested code - please check if it works for you. These are the changes for the AVX512F kernel, you will need similar changes if you are using the other kernels.

@mehdirezaie
Copy link
Author

mehdirezaie commented Sep 6, 2019

Thanks a lot, @manodeep. I do not want to completely change bin_low <= mu < bin_hi for all bins, I would like to change it only for the last bin, bin_max_m_1 <= mu <= bin_max. With the current version of the code, we are missing the pairs that lie exactly at mu=1.0.

@manodeep
Copy link
Owner

manodeep commented Sep 6, 2019

Yup. The changes I suggested above are to simply make the last bin include mu_max, and does not affect any other bin. Here are the easier changes in the Fallback kernel

  • Change to if((zpos - last_z1) > max_all_dz) here
  • Change to > instead of >= here
  • Change to min_dz > pimax here
  • Change to *z1 < target_z here
  • Change to > instead of >= here

All the above changes are needed across ALL kernels. You might have noticed that I mentioned these exact changes in the AVX512F kernel previously.

The final change needed in the Fallback kernel is:

  • Add a line mu_bin = mu == mu_max ? nmu_bins-1:mu_bin; after this line

Hope that makes the required changes clearer.

@rainwoodman
Copy link

Are these changes already included in Corrfunc?

@rainwoodman
Copy link

In terms of generality -- what about adding a new bin that counts 1<=mu < inf? This will minimally break backward compatibility, and the caller can decide what to do about the last bin.

@manodeep
Copy link
Owner

@rainwoodman Your second point made me realise that the final bin within the pairs already count these "overflow" separations. Just that those bin-counts are never returned from the API. For every s bin, the corresponding mu_bin == nmu_bins contains the relevant information -- see here

@mehdirezaie Would that work for what you need? In that case, you will have to alter the for loop in this line to include mu == nmu_bins. The malloc here should already be allocating the necessary (extra) space.

@mehdirezaie
Copy link
Author

Sorry for the slow response @manodeep ... have been applying for postdocs. Yes, that would be much better, especially it would not break the backward compatibility as Yu mentioned. I will work on this over this weekend, and send the code for review.

@manodeep
Copy link
Owner

No worries @mehdirezaie.

Regarding the change, the behaviour you need is different from the remainder of the pair-counters - so I don't think it will be prudent to merge into the main code-base. @lgarrison What do you think?

@lgarrison
Copy link
Collaborator

I think it would be good to provide an option to have the last bin boundary closed. Actually, I think it makes sense to have it be the default behavior, but in the interest of backwards compatibility, I would suggest adding a flag called last_bin_closed (or something like that) with a default value of False, allowing users to opt-in to the new behavior if they want. I realize providing this flexibility would come with the cost of increased kernel complexity, though.

I think adding an extra bin with the "overflow" separations is potentially a more breaking change. Some users will have assumed that the results array is of length nbins, not nbins+1 (reasonably so), and this will break their scripts.

@manodeep
Copy link
Owner

@lgarrison That is a great suggestion -- should be reasonably easy to provide the numpy-like behaviour to the user.

And agree on the overflow bins.

@mehdirezaie
Copy link
Author

Would it not be better to have an option to have the extra bin values to be added to the last bin? I think this would demand less modification to the original code since the pairs at mu==1 are already computed, and we could just add the values of the extra bin to the last bin if needed. This will maintain the dimension of the results array.

@manodeep
Copy link
Owner

@mehdirezaie While that might be okay for DDsmu since there is a natural maxima for the "pi" value (i.e., when mu_max == 1.0), but for all other cases the overflow bin spans the range [pi_max/mu_max, infinity] and adding the pair-counts of the overflow bin to the last bin would not be meaningful.

@lgarrison
Copy link
Collaborator

lgarrison commented Dec 16, 2019 via email

@manodeep
Copy link
Owner

@lgarrison Yup - that's what I understood. Your suggestion is that the last bin also contains pairs with rsep == rmax, or mu==mu_max etc. That's just a one-line check, and an addition to pair-counts, weight-counts etc (I do have reservations about floating point equality tests against non-zero numbers but we can work that out later).

My opposition to adding the overflow bin pair-counts is that the last bin then contains pairs with rsep >= rmax, or pi >= pi_max etc, which can add arbitrary number of pairs, and the results would likely depend on the bin-refine-factors.

@rainwoodman
Copy link

rainwoodman commented Mar 1, 2020 via email

@manodeep
Copy link
Owner

@rainwoodman I am not sure what you mean with numpy histogram putting out-of-range samples in some bins; my understanding is that all out-of-range data are simply not counted. In any case, I would argue that if the user wanted to include the pimax, then they could pass pimax + FLT_EPSILON, and then the pimax values would be included.

However, since there is an upper limit to mu_max, the case is different for DDsmu. Corrfunc already carries overflow bins along both s, and mu. Since the only possible overflow in mu is mu==1.0, simply adding in those bin counts to the final bin would be fine (i.e., what I suggested in this comment).

I am all for @lgarrison suggestion to have a new keyword that makes the final bin a closed interval, but we may have to decide if we only roll that out for the two DDsmu, DDsmu_mocks routines.

@lgarrison What do you suggest?

@lgarrison
Copy link
Collaborator

I think that rolling out a last_mu_bin_closed flag for just the two DDsmu CFs is fine. I think that perhaps the flag could even be on by default, as it's a small change. But I'd like to know what others think.

@rainwoodman
Copy link

rainwoodman commented Apr 13, 2020 via email

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

No branches or pull requests

4 participants