Test against DM implementation of CTInterpolator #84

arunkannawadi · 2024-07-30T21:22:49Z

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.

descwl_coadd/tests/test_coadd_and_psf_coadd_agree.py

esheldon · 2024-07-31T14:30:58Z

descwl_coadd/tests/test_coadd_correct.py

@@ -78,7 +78,7 @@ def test_coadd_image_correct(crazy_wcs, crazy_obj):
    world_origin = galsim.CelestialCoord(0 * galsim.degrees, 0 * galsim.degrees)

    aff = galsim.PixelScale(scale).affine()
-    aff = aff.withOrigin(
+    aff = aff.shiftOrigin(


Why is shiftOrigin more appropriate than withOrigin here?

This has been deprecated since GalSim v2.3. I just noticed it when looking at the GitHub Actions logs when I had a failure in the beginning.

https://github.com/GalSim-developers/GalSim/blob/f6499baeadf433a34f98f5c918b2580e6d28bf7b/docs/older.rst#L92

withOrigin is basically a call to shiftOrigin with a deprecation warning:

https://github.com/GalSim-developers/GalSim/blob/f6499baeadf433a34f98f5c918b2580e6d28bf7b/galsim/wcs.py#L699-L702

beckermr · 2024-07-31T14:45:41Z

descwl_coadd/tests/test_coadd.py

+        assert (coadd_data['coadd_exp'].image.array ==
+                dm_coadd_data['coadd_exp'].image.array).all()


This level of agreement is surprising in that you'd need the set of good pixels around the bad ones to match exactly.

Let's add an assert that some data was actually interpolated in the test. Otherwise things would appear to agree but it'd be because we didn't do any interpolation.

The logic of selecting the bad pixels is exactly the same and involves only integer tuples. The floating point operations are all done by the scipy function, so I am not surprised that the agreement is exact. But I'll add a test to see that some pixels are actually interpolated.

You mentioned in the pr description that the dm version was more efficient since it used some sort of spanset data structure / algorithm. I assumed this might not give the same results in all cases. Is that not true?

Yes, I said that we have a different implementation of this function in the DM codebase: https://github.com/LSSTDESC/descwl_coadd/blob/master/descwl_coadd/interp.py#L13-L14

Once we get the location of good pixels and bad pixels, we call the scipy.CloughTocher2DInterpolator the same way as done here.

FWIW, here's the C++ implementation of the function that behaves identical to the function I posted above: https://github.com/lsst/meas_algorithms/blob/main/src/CloughTocher2DInterpolatorUtils.cc

Good catch with the pixels not being interpolated in the first place, Matt. I'm going to do some more local testing and raise it again later. Marking it as a draft for now.

descwl_coadd/tests/test_coadd.py

arunkannawadi · 2024-09-06T20:51:29Z

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.

Since this uses LSST DM code, which has GPL v3 LICENSE, this repo will have to have the same license.

arunkannawadi · 2024-09-16T16:32:10Z

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

esheldon · 2024-09-16T18:39:16Z

descwl_coadd/tests/test_coadd.py

+                    actual=dm_coadd_data['coadd_exp'].image.array,
+                    desired=coadd_data['coadd_exp'].image.array,
+                    atol=0.0,
+                    rtol=1653.16,


Why do we allow such a large rtol?

All of these values are empirically set based on observed differences and are intended to act as benchmarking rather than tests. For small pixel values close to zero, this shouldn't be surprising but this is just how sensitive CTInterpolator is to the choice of coordinate frames.

Reducing rtol and running the tests gives me this.

this test is better specified in atol anyways - the values cross zero

descwl_coadd/tests/test_coadd.py

beckermr · 2024-09-16T21:18:57Z

descwl_coadd/tests/test_coadd.py

+                    actual=dm_coadd_data['coadd_exp'].image.array,
+                    desired=coadd_data['coadd_exp'].image.array,
+                    atol=0.0,
+                    rtol=1653.16,


this test is better specified in atol anyways - the values cross zero

arunkannawadi · 2024-09-16T21:28:53Z

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 · 2024-09-16T21:46:42Z

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 · 2024-09-16T22:00:45Z

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 · 2024-09-16T22:03:04Z

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

arunkannawadi · 2024-09-25T13:43:10Z

@esheldon @beckermr am I good to merge this?

esheldon · 2024-09-25T14:00:04Z

LGTM

arunkannawadi · 2024-09-30T18:58:02Z

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 · 2024-09-30T19:00:25Z

Thanks and sorry for missing the formal lgtm!

arunkannawadi · 2024-09-30T19:10:11Z

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 · 2024-09-30T19:18:30Z

Yep - will get fixed by conda-forge/admin-requests#1094

Just merged - will be maybe 30 mins

arunkannawadi force-pushed the DM_CT_interpolator branch 4 times, most recently from 59bb716 to 93d047d Compare July 30, 2024 23:50

arunkannawadi marked this pull request as ready for review July 30, 2024 23:54

arunkannawadi requested a review from esheldon July 30, 2024 23:54

esheldon reviewed Jul 31, 2024

View reviewed changes

beckermr requested changes Jul 31, 2024

View reviewed changes

arunkannawadi force-pushed the DM_CT_interpolator branch 3 times, most recently from 10b2217 to b6a6c0f Compare July 31, 2024 16:49

beckermr reviewed Jul 31, 2024

View reviewed changes

descwl_coadd/tests/test_coadd.py Outdated Show resolved Hide resolved

arunkannawadi marked this pull request as draft July 31, 2024 20:14

arunkannawadi force-pushed the DM_CT_interpolator branch from cfa08c7 to a513517 Compare September 6, 2024 20:51

arunkannawadi force-pushed the DM_CT_interpolator branch 3 times, most recently from 568ee55 to c029a46 Compare September 16, 2024 15:09

arunkannawadi added 5 commits September 16, 2024 12:31

Avoid deprecated warning from GalSim

dea4c7a

Add bad columns to interpolator smoke test

296599d

Run interpolator on just the maskedImage

5510037

Test against DM implementation of CTInterpolator

344079f

Create LICENSE

65f1e41

Since this uses LSST DM code, which has GPL v3 LICENSE, this repo will have to have the same license.

arunkannawadi force-pushed the DM_CT_interpolator branch from a28ea7b to 65f1e41 Compare September 16, 2024 16:31

arunkannawadi marked this pull request as ready for review September 16, 2024 16:31

arunkannawadi requested review from esheldon and beckermr September 16, 2024 16:32

esheldon reviewed Sep 16, 2024

View reviewed changes

beckermr requested changes Sep 16, 2024

View reviewed changes

Check mfrac for both values of bad_columns

99259e7

arunkannawadi requested a review from beckermr September 16, 2024 21:31

Remove comparisons involving only rtol

eef0f00

arunkannawadi merged commit 2e771d7 into master Sep 30, 2024
1 check passed

arunkannawadi deleted the DM_CT_interpolator branch September 30, 2024 18:58

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Test against DM implementation of CTInterpolator #84

Test against DM implementation of CTInterpolator #84

arunkannawadi commented Jul 30, 2024 •

edited

Loading

esheldon Jul 31, 2024

arunkannawadi Jul 31, 2024

arunkannawadi Jul 31, 2024

beckermr Jul 31, 2024 •

edited

Loading

arunkannawadi Jul 31, 2024

beckermr Jul 31, 2024

arunkannawadi Jul 31, 2024

arunkannawadi Jul 31, 2024

arunkannawadi Jul 31, 2024

arunkannawadi commented Sep 6, 2024 •

edited

Loading

arunkannawadi commented Sep 16, 2024

esheldon Sep 16, 2024

arunkannawadi Sep 16, 2024

arunkannawadi Sep 16, 2024

beckermr Sep 16, 2024

beckermr Sep 16, 2024

arunkannawadi commented Sep 16, 2024

beckermr commented Sep 16, 2024

arunkannawadi commented Sep 16, 2024

beckermr commented Sep 16, 2024

arunkannawadi commented Sep 25, 2024

esheldon commented Sep 25, 2024

arunkannawadi commented Sep 30, 2024

beckermr commented Sep 30, 2024

arunkannawadi commented Sep 30, 2024

beckermr commented Sep 30, 2024

		assert (coadd_data['coadd_exp'].image.array ==
		dm_coadd_data['coadd_exp'].image.array).all()

Test against DM implementation of CTInterpolator #84

Test against DM implementation of CTInterpolator #84

Conversation

arunkannawadi commented Jul 30, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

beckermr Jul 31, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

arunkannawadi commented Sep 6, 2024 • edited Loading

arunkannawadi commented Sep 16, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

arunkannawadi commented Sep 16, 2024

beckermr commented Sep 16, 2024

arunkannawadi commented Sep 16, 2024

beckermr commented Sep 16, 2024

arunkannawadi commented Sep 25, 2024

esheldon commented Sep 25, 2024

arunkannawadi commented Sep 30, 2024

beckermr commented Sep 30, 2024

arunkannawadi commented Sep 30, 2024

beckermr commented Sep 30, 2024

arunkannawadi commented Jul 30, 2024 •

edited

Loading

beckermr Jul 31, 2024 •

edited

Loading

arunkannawadi commented Sep 6, 2024 •

edited

Loading