Why is diffuse only tested to tolerance of 1e-2?

[dannyjacobs]

I am making a summary of all our unit tests and notice that the analytic diffuse is only toleranced to 1e-2. This seems rather large. Does it have to do with the number of terms used in the analytic expansion?

Here is the line https://github.com/RadioAstronomySoftwareGroup/pyuvsim/blob/9179bad3573856cefc03b7adfc9e5e23e8b54520/pyuvsim/tests/test_run.py#L129

Apologies to @aelanman, @bhazelton and @mkolopanis if this was already discussed in the diffuse test PR (#309). I read that thread and didn't see any mention but could have missed it.

[aelanman]

I don't think this was discussed. Playing around with it now, the largest error is on the non-flat monopole model, and this error drops if I increase the expansion order, so I think it is the expansion order.

[dannyjacobs]

The rest of the point source sims (with the exception of zenith for some reason) are tested to 1e-8. Can you get it that low with enough terms and still run in a reasonable time?

[aelanman]

If I increase the order to 50, it reaches the same error level as the other models, which is around 1e-4. I see a similar error when I increase the nside to 256 the error remains the same. This might take some time to track down.

[dannyjacobs]

To which other models do you refer? All the point source tests test to the default allclose atol of 1e-8 except for the original "pt_src_zenith" which is at 5e-3 for some reason. What tolerance do you expect given the current status of the diffuse analysis?

[aelanman]

The other diffuse models (the test is parametrized). The nonzero-w test with the monopole is the only solution there that has an expansion order. The others have maximum errors of around 1e-4. Increasing the expansion order on the monopole solution brings its error down to about there.

[dannyjacobs]

Oh, I think I see. The test iterates through the various models https://github.com/RadioAstronomySoftwareGroup/pyuvsim/blob/9179bad3573856cefc03b7adfc9e5e23e8b54520/pyuvsim/tests/test_run.py#L68 If you set https://github.com/RadioAstronomySoftwareGroup/pyuvsim/blob/9179bad3573856cefc03b7adfc9e5e23e8b54520/pyuvsim/tests/test_run.py#L85 to 50, all tests pass at 1e-4. The other functional forms, cosza, don't pass below 1e-4 either. Is that what you expect?

[aelanman]

It's unclear. Based on the results I'm getting with healvis for the paper, the biggest errors for the cosza model are around 1e-6 for nside 256, so it's not unreasonable to see errors around 1e-4 for a lower res simulation.

(This plot shows errors vs. baseline length at three resolutions. These are errors relative to the value at |u| = 0)

The best theoretical error bounds I have are shown in the following table. They're the error bounds on a Monte Carlo integration of the RIME with the same number of pixels (maximized over all baselines). This is, of course, an overestimate of the error, since HEALPix pixels are not random. Steven and I have done a lot of reading on other approaches to getting error bounds, but so far none have been very successful.

Sorry for the formatting.

model  nside  bound
cosza 128 0.0018414265108688966
gauss_a0.01 128 7.557956298980507e-06
gauss_a0.5 128 0.00045461018164161916
gauss_a1.57 128 0.0015982193802132737
monopole 128 0.0031894396159931728
polydome 128 0.0012054956881261793
projgauss_a0.05 128 7.971105723083991e-05
projgauss_a0.5 128 0.0007699381020527679
projgauss_a2.5 128 0.001730271451150872
xysincs_a64.0 128 6.190377869690528e-05

[steven-murray]

@aelanman how did you derive the numbers in the final table?

@dannyjacobs yeah this is a tricky problem. As @aelanman said, we have tried a few approaches to generating upper limits on the errors given a set of input nodes of the integration, but none seem to match what we're getting. We're thinking now to defer that to a future paper, and be content with presenting the models themselves and just showing their actual errors for heuristics. Potentially then these tests could be updated to be more specific/accurate.

One thing to mention is that it is not clear that choosing healpix nodes for integration gives the minimum error. With a theoretical upper limit model we might be able to make that claim (or point to a better pixelisation).

[aelanman]

@steven-murray Those numbers come from the error bound on MC integration : error < Var{f} / sqrt(N). N is, of course, just the number of sample points. Var{f} is estimated as max_{u,v} ( var[ f(l, m) exp(- 2 π i (ul + vm) ] ), where var is taking the variance of its argument over the sample points and it's maximized over all baselines that were used in simulations.

It's clear that there has been a lot of effort to find better pixellizations of spheres for integration, but for this paper I really think we just want to know that the remaining error is reasonable.

[dannyjacobs]

According to the updated analysis, the residuals no longer diverge under any tested resolution.

[bhazelton]

I'm not sure where we are on this, @dannyjacobs last comment suggests improvement but I'm not sure how much. I played around with this some, I set the non-flat monopole order to 50 and I could get the tests to pass with the following tolerances: ('monopole', 6e-4), ('cosza', 2e-4), ('quaddome', 1e-4), ('monopole-nonflat', 3e-4) If I increase the side to 256, the tests run a little slower but I can get them to pass with: ('monopole', 3e-4), ('cosza', 2e-4), ('quaddome', 1e-4), ('monopole-nonflat', 3e-4)

Is this in line with expectations? They seem a little higher than the plots in Danny's last comment.

[steven-murray]

So, in the paper we didn't really come to any solution for how to express the expected residual error as a function of sky-model and nside. Turns out that is wickedly tricky. I think that evaluating the analytic solutions to the highest precision that is feasible is a good idea here, and then we just play around empirically to test the accuracy.