LSSTDESC / descwl_coadd

Tools for coadding using the LSST DM software stack
GNU General Public License v3.0
2 stars 2 forks source link

Test against DM implementation of CTInterpolator #84

Closed arunkannawadi closed 3 weeks ago

arunkannawadi commented 2 months ago

Since the interpolator objects can now be specified, I have added a test in this PR to compare that the CT interpolation implemented in DM behaves identical to the implementation here. The implementation differs in that i) the computationally expensive search function under the hood that finds good pixels around bad pixels is implemented in C++ instead of numba for policy reasons and ii) it uses the SpanSet data structure to find the bad pixels more efficiently than the multipass brute-force search implemented here. The test here ensures that the end results are the same and will highlight any regressions in the future.

arunkannawadi commented 1 month ago

It turns out that the Delaunay triangulation happening under the hood is not invariant under x-y flip (not too surprising). The DM implementation uses (x, y) notation which is flipped relative to the [row, col] notation used with the numpy arrays here. This results in slightly different values in the output of scipy.CloughTocher2DInterpolator. I've included a PR in meas_algorithms to support either modes. Both modes are self-consistent and I don't think we can say that one is better than the other. The tests should run fine once a weekly release with that PR is available on stackvana.

arunkannawadi commented 1 month ago

w37 appeared on stackvana over the weekend, which means this PR is ready to be reviewed and merged.

arunkannawadi commented 1 month ago

Why does values crossing zero matter? rtol is the ratio between the absolute value of the difference to the absolute value of the desired value.

beckermr commented 1 month ago

Here is the formula from the numpy docs:

absolute(a - b) <= (atol + rtol * absolute(b))

If b is close to zero and atol=0, then rtol may need to be quite large to acomodate a modest change in the absolute value. For example, if a=1e-10 and b=1e-16, then you'd need rtol > |a-b|/b ~ 9.99e5 for the test to pass. On the other hand, an absolute tolerance of ~1e-7 might be more sensible if you think both numbers effectively mean zero.

In this case, the image values are of order -1 to 1 and sometimes quite small. Thus we hit the math above near zero quite a bit, causing the need for a large rtol. Using an absolute tolerance of a few tenths is more sensible and expresses better what is happening.

arunkannawadi commented 1 month ago

That's all consistent with what I told Erin. There's a test with atol as well in addition. If this rtol is too high, then this is as good as not having a test with rtol anyway. Is the concern that a future version of numpy or scipy might make this test fail?

beckermr commented 1 month ago

Yes, I suspect the test based on atol will be more robust.

arunkannawadi commented 1 month ago

@esheldon @beckermr am I good to merge this?

esheldon commented 1 month ago

LGTM

arunkannawadi commented 3 weeks ago

I'll merge this despite the stale 'Changes requested' from @beckermr since i) the comments have been addressed and ii) this is largely an addition to unit tests and not a change to the behavior of the code

beckermr commented 3 weeks ago

Thanks and sorry for missing the formal lgtm!

arunkannawadi commented 3 weeks ago

Well, we do need your help now. GHA is failing with

InvalidSpec: The package "conda-forge/noarch::mpi==1.0=mpich" is not available for the specified platform

that was not the case about a week ago.

beckermr commented 3 weeks ago

Yep - will get fixed by https://github.com/conda-forge/admin-requests/pull/1094

Just merged - will be maybe 30 mins