manodeep
Pairs at μ=1.0
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.
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 > pimaxhere - Change to
*localz1 < target_zhere - 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_maxbin 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.
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.
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 > pimaxhere - Change to
*z1 < target_zhere - 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.
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.
@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.
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.
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?
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.
@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.
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.
@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.
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 [email protected] 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 [email protected] Flatiron Research Fellow, Cosmology X Data Science Group Center for Computational Astrophysics, Flatiron Institute lgarrison.github.io
@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.
On Mon, Dec 16, 2019 at 3:04 PM Manodeep Sinha [email protected] 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 .
@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?
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.
Yes. You are right. I confused histogram with digitize.
On Mon, Apr 13, 2020 at 5:54 AM Manodeep Sinha [email protected] 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 .