Pairs at μ=1.0

[mehdirezaie]

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]

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]

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]

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]

Are these changes already included in Corrfunc?

[rainwoodman]

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]

@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]

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]

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]

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]

@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]

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]

@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]

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 ***@***.***> 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 <#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

[manodeep]

@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]

On Mon, Dec 16, 2019 at 3:04 PM Manodeep Sinha ***@***.***> 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 <#194?email_source=notifications&email_token=AABBWTEQTWU6QGAAYYMR2J3QZACQFA5CNFSM4IUCY7N2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEHAN4LQ#issuecomment-566287918>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AABBWTHM54YDYT6OUJL3423QZACQFANCNFSM4IUCY7NQ> .

[manodeep]

@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]

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]

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

…

On Mon, Apr 13, 2020 at 5:54 AM Manodeep Sinha ***@***.***> 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 <#194 (comment)>). 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 <#194 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AABBWTEDPR4BBVOZVQMTHT3RMMDO5ANCNFSM4IUCY7NQ> .