search

manodeep / Corrfunc

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

Pairs at μ=1.0 #194

Open mehdirezaie opened 4 years ago

mehdirezaie commented 4 years ago

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.

manodeep commented 4 years ago

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)

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 commented 4 years ago

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 commented 4 years ago

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

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:

Hope that makes the required changes clearer.

rainwoodman commented 4 years ago

Are these changes already included in Corrfunc?

rainwoodman commented 4 years ago

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 commented 4 years ago

@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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

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 commented 4 years ago

@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 commented 4 years ago

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 commented 4 years ago

@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 commented 4 years ago

I think I'd agree with including separations in r that match the end of the last bin in general (someday!). And to be clear, that's what I was proposing the new option do with pi values that match the end of the last pi bin.

On Mon, Dec 16, 2019 at 5:25 PM Manodeep Sinha notifications@github.com wrote:

@mehdirezaie https://github.com/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.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/manodeep/Corrfunc/issues/194?email_source=notifications&email_token=ABLA7S5AA75NBXRKW3VSO3TQY757NA5CNFSM4IUCY7N2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHAKQVI#issuecomment-566274133, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABLA7S6BJMYGPFFUKMDBY2DQY757NANCNFSM4IUCY7NQ .

-- Lehman Garrison lgarrison@flatironinstitute.org Flatiron Research Fellow, Cosmology X Data Science Group Center for Computational Astrophysics, Flatiron Institute lgarrison.github.io

manodeep commented 4 years ago

@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 commented 4 years ago

On Mon, Dec 16, 2019 at 3:04 PM Manodeep Sinha notifications@github.com wrote:

@lgarrison https://github.com/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.

But the caller can decide to drop it with [:-1] --not exactly a huge pain. numpy's histogram functions also put the out of range samples to these extra bins, so an acute user shouldn't be too surprised.

—

You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/manodeep/Corrfunc/issues/194?email_source=notifications&email_token=AABBWTEQTWU6QGAAYYMR2J3QZACQFA5CNFSM4IUCY7N2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHAN4LQ#issuecomment-566287918, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABBWTHM54YDYT6OUJL3423QZACQFANCNFSM4IUCY7NQ .

manodeep commented 4 years ago

@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 commented 4 years ago

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 commented 4 years ago

Yes. You are right. I confused histogram with digitize.

On Mon, Apr 13, 2020 at 5:54 AM Manodeep Sinha notifications@github.com wrote:

@rainwoodman https://github.com/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 https://github.com/manodeep/Corrfunc/issues/194#issuecomment-555739507).

I am all for @lgarrison https://github.com/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 https://github.com/lgarrison What do you suggest?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/manodeep/Corrfunc/issues/194#issuecomment-612885985, or unsubscribe https://github.com/notifications/unsubscribe-auth/AABBWTEDPR4BBVOZVQMTHT3RMMDO5ANCNFSM4IUCY7NQ .