Fix WL functions in pyccl.halos.profiles

[hsinfan1996]

In pyccl.halos.profiles, convergence, shear, reduced_shear and magnification can now handle array_like a_source.

[coveralls]

Pull Request Test Coverage Report for Build 2510576231

• 6 of 6 (100.0%) changed or added relevant lines in 1 file are covered.
• No unchanged relevant lines lost coverage.
• Overall coverage increased (+0.004%) to 97.104%

---

Totals
Change from base Build 2463247276:	0.004%
Covered Lines:	4493
Relevant Lines:	4627

---

💛 - Coveralls

[damonge]

@hsinfan1996 thanks for this. Just so I'm clear, what fails exactly when you make a_source an array?

[hsinfan1996]

@hsinfan1996 thanks for this. Just so I'm clear, what fails exactly when you make a_source an array?

In convergence and shear, the first step self.projected(cosmo, r, M, a_lens, mass_def) requires a_lens to be a float. The next step sigma_critical(cosmo, a_lens, a_source) requires a_lens and a_source to have the same dimension, so both of them have to be floats given that a_lens is a float. If a_source is an array, sigma_critical throws an error.

[damonge]

OK, I see! For practical applications, how inconvenient is it that a_lens can only be a float? Just trying to figure out if it's worth trying to enable that in this PR.

[hsinfan1996]

OK, I see! For practical applications, how inconvenient is it that a_lens can only be a float? Just trying to figure out if it's worth trying to enable that in this PR.

From what I know, that a_lens can only be a float does not affect any relevant calculations in CLMM.

[damonge]

Humm, but if I understand correctly, from now on, a_source will always need to have the same shape as r, right? I think this is very non-desirable. Could we force both a_source and a_lens to be scalars in that case?

[hsinfan1996]

Humm, but if I understand correctly, from now on, a_source will always need to have the same shape as r, right? I think this is very non-desirable. Could we force both a_source and a_lens to be scalars in that case?

For array_like r, a_source can still be a float and the output will still have the same shape as r. The result is the same as passing array_like a_source with the dimension of r.

[damonge]

OK, but the only two possibilities are:
a) a_source is a float
b) a_source is an array of the same shape as r.

I think option b) above is not useful at all (and is not documented as such), so I can see two ways forward:

1. We force a_source and a_lens to be always floats.
2. We allow a_source to be an array, but then the output of these functions should be a 2D array of shape [len(a_source), len(r)], so that r and a_source don't need to be linked up.

Option 1 above seems to be the simplest one, so I'd advocate for that unless you feel like giving the second one a go.

[hsinfan1996]

OK, but the only two possibilities are: a) a_source is a float b) a_source is an array of the same shape as r.

I think option b) above is not useful at all (and is not documented as such), so I can see two ways forward:

1. We force a_source and a_lens to be always floats.
2. We allow a_source to be an array, but then the output of these functions should be a 2D array of shape [len(a_source), len(r)], so that r and a_source don't need to be linked up.

Option 1 above seems to be the simplest one, so I'd advocate for that unless you feel like giving the second one a go.

Actually, b) is what I wanted to achieve in this PR and is already very useful to us. In addition, it only requires a few lines of code added to pyccl to allow array_like a_source, which is said to be supported in the documentation. I agree that we can make the functions more versatile by allowing different shapes of a_source and r, but I think it will need c-level changes (not 100% sure).

[damonge]

Humm, I don't really understand how b) is useful at all. Why is it useful to have the function return e.g. convergence as a function of (r,a_source) where each r corresponds to a different a_source?

Wouldn't it be more useful to have it as a function of both r and a_source independently of each other?

[damonge]

(I don't think the latter would require c-level modifications, BTW, that's why I'm asking)

[eduardojsbarroso]

For what I understand, in CLMM, we want to make computations considering several sources with different projected radius but in the same redshift, therefore same a_source, or to make computations for several sources in different redshifts. That is why we wanted the case b). Either each projected radius must have a pair, thus len(r) = len(a_source), or we are giving several projected radius but we want to analyze all sources in the same redshift, which then requires a_source = [a_source]*len_r.

I feel like the case where, for example, r = [1,2,3] and a_source = [0.5, 0.4] is not very usefull because it is not clear what this stands for. When we give an array r and a float a_source, we can consider that all the sources are in the same redshift, but in the case where we give an array a_source that does not match the array r, how can we know in which redshift are all the sources represented by r?

We did implement these changes in CLMM after importing the CCL functions but we've come up with two problems: 1) the vectorization on our end let the code very slow. 2)Even though the calculations were feasible, we got a failled test because a_source was an array. And that is why we opened this issue.

[damonge]

Thanks @eduardojsbarroso
OK, just wanna make sure I understand this (and I'm afraid I still don't), so when you say:

I feel like the case where, for example, r = [1,2,3] and a_source = [0.5, 0.4] is not very usefull because it is not clear what this stands for

This would give you the profile at r=[1, 2, 3] for two sources, one at a=0.5 and one at a=0.4. Isn't this more useful than e.g. having the profile for e.g. r=1 at a=0.4, r=2 at a=0.5, and r=3 at a=0.6?

I had assumed you'll want to do some line-of-sight integrals that involve integrating over source redshifts for the same lens, and that's why I thought that having r and a_source be independent would be useful, but I'm probably wrong. Can I check what calculation is it that benefits from having r and a_source coupled?

[damonge]

(I should say, I don't want to delay this getting merged, just want to make sure I'm on the same page!)

[hsinfan1996]

Thanks @eduardojsbarroso OK, just wanna make sure I understand this (and I'm afraid I still don't), so when you say:

I feel like the case where, for example, r = [1,2,3] and a_source = [0.5, 0.4] is not very usefull because it is not clear what this stands for

This would give you the profile at r=[1, 2, 3] for two sources, one at a=0.5 and one at a=0.4. Isn't this more useful than e.g. having the profile for e.g. r=1 at a=0.4, r=2 at a=0.5, and r=3 at a=0.6?

I had assumed you'll want to do some line-of-sight integrals that involve integrating over source redshifts for the same lens, and that's why I thought that having r and a_source be independent would be useful, but I'm probably wrong. Can I check what calculation is it that benefits from having r and a_source coupled?

https://github.com/LSSTDESC/CLMM/blob/fc813c348fceab0b009f39bd4472af244fba9693/clmm/theory/parent_class.py#L490
Our calculation is almost identical to what is in CCL. The difference is that passing array_like z_src (it will be converted to a_source for CCL) will work. We are trying to replace this function (_eval_convergence) with pyccl.halos.convergence for our code when using CCL as the backend. In our use case, one source has one position (r) at most, so the dimension of z_src (a_source) is coupled with that of r.

I still agree that allowing a_source and r to have different shapes makes the function more versatile, but this is beyond what we need for now. I was thinking that this PR would be an easy fix to make the function work as stated in the documentation. Most CCL functions (if not all) in the background module require a_lens and a_source to have the same dimension. In addition, vectorized functions (e.g. angular_diameter_distance_vec) only work with 1D arrays. Thus, if we want to pass a 2D array of a_source, this probably involves flattening the input and reshaping the output in python level (so no c level change).

[damonge]

Thanks @hsinfan1996 @eduardojsbarroso , if this is actually what you guys need, let's move on with it. Thanks for indulging me!

Just one final query to fix the docs while we're at it and we can merge

[hsinfan1996]

Thanks @hsinfan1996 @eduardojsbarroso , if this is actually what you guys need, let's move on with it. Thanks for indulging me!

Just one final query to fix the docs while we're at it and we can merge

@damonge Thank you for reviewing and the conversation. @nikfilippas Thank you.

[damonge]

LGTM!