Swap out ddfacet degridder for Nifty degridder

[bennahugo]

The time has come for optimization. The goal is to replace the DDFacet degridder with the NIFTY row parallel degridder.

I will incorporate the necessary transforms for non-coplanar faceting without w correction because the nifty gridder w kernels are wrong in the faceting case (l0 and m0 are not 0 in that case).

The approach is to apply the inverse of:
https://github.com/ratt-ru/bullseye/blob/master/bullseye/mo/cpu_gpu_common/baseline_transform_policies.h#L70
to the uvw coordinates per facet (trivially parallel with another numba kernel), and apply
https://github.com/ratt-ru/bullseye/blob/master/bullseye/mo/cpu_gpu_common/phase_transform_policies.h#L141
to the visibilities. This last bit is also trivially parallel, but quite expensive to compute in comparison (a row x chan matrix of complex exponentials). There is a cheat if the channelization is regular which can be checked and applied in the code path.

I will do both with Numba.

[bennahugo]

This defines the faceting scheme implemented by Eric Greisen in AIPS, based on the scheme by Perley and Cornwell 1992 and is the most trivial form of non-coplanar imaging. It will perform very well (better than w-projection approaches), because the number of facets are small because of the specification of regions. It is in essence a variation on the targetted faceting regime.

[mreineck]

I'm not sure whether this issue is still active, but perhaps it is of interest that wgridder (the successor of nifty_gridder) does support shifting the phase center now and also got a few performance tweaks. So if you are still interested in this, please have a look at https://gitlab.mpcdf.mpg.de/mtr/ducc!

[bennahugo]

Hi @mreineck thanks. The approach I have taken has now been adopted as DD calibration system for a number of projects, so this would be useful.

One would need an additional modification to the approach taken in the nifty gridder (or its successor) in that the UVW coordinates should be rotated such that each facet is tangent to the celestial sphere plus a phase shift to the facet centres, therefore no w-stacking or w-projection would be needed.

Since we are doing targetted faceting in small patches this would mean that we only apply anti-aliasing filters and would be more optimal than computing w kernels. See https://github.com/ratt-ru/codex-africanus/blob/master/africanus/gridding/perleypolyhedron/policies/baseline_transform_policies.py and https://github.com/ratt-ru/codex-africanus/blob/master/africanus/gridding/perleypolyhedron/gridder.py#L92 for a working implementation of this

This is one of the core reasons why I haven't yet adopted the nifty gridder so would be nice feature request if you have time.

The successor to this would not be implemented in cubical. @JSKenyon is working on a successor package to cubical based on Dask, although much of the machinery here to mark and split regions would remain the same.

[mreineck]

One would need an additional modification to the approach taken in the nifty gridder (or its successor) in that the UVW coordinates should be rotated such that each facet is tangent to the celestial sphere

Isn't this something that can be done as a preprocessing step? I try to keep the core algorithm as small as possible and try to omly integrate features which cannot be moved outside without a big performance loss.

Since we are doing targetted faceting in small patches this would mean that we only apply anti-aliasing filters and would be more optimal than computing w kernels.

Fine, that would simply mean using the wgridder in narrow-field mode (i.e. do_wgridding=False)

At the moment my main focus is on porting the algorithm to GPUs, and this doesn't leave me a lot of spare time ;-)

[JSKenyon]

Thanks for the replies @mreineck! As @bennahugo mentioned, QuartiCal will hopefully begin to replace CubiCal soon. Degridding is an important feature which we are currently missing but I am very keen to try out the wgridder. I am presently working on some data format stuff but thereafter I may give it a go. @landmanbester already has some thoughts on the matter and (practical concerns notwithstanding) it shouldn't be too tough!

[mreineck]

Great! I'll most likely talk to @landmanbester next week, then we can discuss this in more detail.

[bennahugo]

Hmm.. I don't think one would be able to do this as a pre-processing step -- the phases have to be steered with the original uv coordinates which means one would need to keep NFacet x DATA columns in memory. But thanks for the input

[mreineck]

I wouldn't be surprised if I was completely wrong about the preprocessing ;-)
But next week this shold become clearer in any case!

[landmanbester]

Yeah, let's discuss next week. I would be interested to see how the wgridder compares to non-coplanar facetting in terms of accuracy and speed. Since we are degridding small patches the overhead of including the w-term in the wgridder may be small, especially if low accuracy is required. But let's try to set up a benchmark to compare

[mreineck]

I just noticed something interesting: the necessary phase adjustments for the visibilities seem to be practically identical for both the shifting (https://github.com/mreineck/ducc/blob/1548bad07fac2f4331b551d7e239401d572a88e0/src/ducc0/wgridder/wgridder.h#L1028) and the rotating (https://github.com/ratt-ru/codex-africanus/blob/29c463c7e79eb8d2a2703d3381448d96c6840eb3/africanus/gridding/perleypolyhedron/policies/phase_transform_policies.py#L49) approach. So with a bit of luck we might be able to pull this of with only preprocessing after all ...
Just one question: are the uvw in phase_rotate the original ones or already the ones rotated by uvw_rotate?

[bennahugo]

Hi @mreineck. I had a quick look through the codebase wrt. to the shifting :-:
It looks like you have the correct idea here, but the following looks incorrect to me: https://github.com/mreineck/ducc/blob/1548bad07fac2f4331b551d7e239401d572a88e0/src/ducc0/wgridder/wgridder.h#L1028
This assumes a narrow angle approximation on the shift which will break for RA near the celestial poles (here the pixel sizes are especially non-linear in RA). The correct way to do this is to use cosines. A linear stepping, as currently implemented, will show as offsets from the requested phase centres as a function of declination.

As for the question:
Yes they are the original -- https://github.com/ratt-ru/codex-africanus/blob/29c463c7e79eb8d2a2703d3381448d96c6840eb3/africanus/gridding/perleypolyhedron/gridder.py#L83-L93
These are rephased from the original w / n direction.

Btw:
I was under the mistaken impression that you don't rephase visibilities, from my discussion with @landmanbester. Here it looks like you do though? Is it possible to combine this with w-stacking to limit the number of wplanes needed for a particular facet -- like is done in DDF with wprojection kernels computed around the rephased centre of each of the images. If so one can perhaps forgo the non-coplanar approach that I indicated for a coplanar approach with limited w planes.
As I indicated to @landmanbester one cannot simply turn w-corrections off without doing something about the geometry -- this will lead to requiring many many small facets far away from the original phase centre of the array to limit the error at the edge of the FoV of each facet.
Edit here is the figure that shows the geometry problem arising from turning w corrections off while keeping the image planes co-planar (Perley 1999):

[mreineck]

Hi @bennahugo, the code you point out is an implementation of the tiling approach discussed in section 2.3 of https://arxiv.org/pdf/1407.1943. This involves rephasing the visibilities, and indeed I had forgotten that I do this already. That's why I had (unjustified) reservations about adding this kind of functionality further up in this issue.
However it seems that the rephasing required for your approach is different from what I'm doing at the moment. As of now, the w-gridder does not care at all about RA and Dec of the observation, and therefore does not even have this information. If we need absolute positions on the celestial sphere for rotating the facets, we need further extensions to the interface.

Is it possible to combine this with w-stacking to limit the number of wplanes needed for a particular facet -- like is done in DDF with wprojection kernels computed around the rephased centre of each of the images.

I think that this was exactly the point of this method. At least that's how I interpret the section of the wsclean paper. If that addresses your concern and we can stick to this method, I'd of course be extremely happy :-)
Running a few real-world tests is probably the best way forward.

[mreineck]

(I should perhaps add that section 2.3 of the wsclean paper was not originally meant for tiling a large image into smaller parts, but it is now being used for that purpose in wsclean.)

[bennahugo]

Thanks @mreineck . The reservation I have is how you interpret delta l and delta m in the wsclean paper -- it is not just a linear scale stepping with the pixel size -- it is a difference of cosines. This makes little difference when radec is on or near the equator on the equatorial frame, however it gives large offsets at the poles as I explained (a requested offset in RA of a degree or so can be off by many arcminutes at say declination of 80 with the former approach).

Eqn 8 does indeed make the w correction dependent on the rephased centre, which is the correct thing to do to limit the number of planes.

It looks to me like the only change required is then to fix the interpretation of the rephasing operation?

[bennahugo]

Re tests, as suggested to @landmanbester would be to:

Target offsets (say up to 5 degrees from the original phase centre) in the sky for high declination positions on the sphere. The error at the edge of the FoV can then be given by evaluating the backward step and looking at the smearing- induced drop in amplitude. One would place a (unity) point source at the edge of the facet, itself placed far away and evaluate the smearing level.

I don't think real data is actually a good test to quantify error though -- I would stick to simulation here

[bennahugo]

I can outline the algorithm for you here:
If you only have an offset in number of pixels use the WCS library to get the corresponding ra dec in radians of that pixel w.r.t the original image centre, see https://docs.astropy.org/en/stable/wcs/index.html.

Once you have the offset position in radians then use the identities to derive the cosines off the original phase direction. See https://github.com/ratt-ru/fundamentals_of_interferometry/blob/master/3_Positional_Astronomy/3_4_direction_cosine_coordinates.ipynb for a full discussion.

[mreineck]

Maybe we are talking at cross purposes here ...
I'm pretty sure that the code in its current form is doing exactly what it is meant to do, even far away from the (original) phase centre.
Given a measurement set, one could in principle ask the gridder to generate an image of the full hemisphere, and one would specify a desired angular pixel size. But this pixel size will of course only correspond to that angle near the center of the image; further away the resolution will become worse.
Computing the image of the full hemisphere is of course very expensive. What the wgridder allows you to do is to specify any rectangular part of the full image (by specifying the number of pixels in and y direction, and a shift in l and m), and it will compute this part as efficiently as it can. But the provided pixel size still is the pixel size at the original phase center!
This approach looks a bit weird at first, but it allows you to compute a full sky image in arbitrary individual tiles and just piece them together at the end without needing any transformation.

[bennahugo]

@mreineck I'm referring to getting the correct projected point on the linearly stepped image plane, not changing the stepping on the plane from linear to non-linear. The requested l,m coordinates are not radians (only in a small angle approximation is a radian approximately the same as a cosine). As implemented now the targetted offset will give you the wrong projected position on the linearly-stepped image plane.

[mreineck]

I agree completely that the l,m coordinates are not radians, they are Cartesian coordinates, and l**2+m**2=1 at the horizon.

Perhaps it is easier to describe what I mean using code: https://github.com/mreineck/ducc/blob/1548bad07fac2f4331b551d7e239401d572a88e0/python/test/test_wgridder.py#L196 shows the unit test that makes sure our faceting gridder is doing what we expect it to do, even for very wide fields (i.e. no matter whether we grid a certain field in one go or in 2x2 or 2x4 tiles, we get the same result within the specified error ranges).

If your interpretation of delta_l and delta_m is different, would it be possible to translate it into ours before calling the gridder? I'd rather not introduce a dependency on WCS for just this purpose (having no external dependencies feels really good :-)

[bennahugo]

Ah wait ok I misunderstood the code - I see you pass the requested positions in: https://github.com/mreineck/ducc/blob/1548bad07fac2f4331b551d7e239401d572a88e0/src/ducc0/wgridder/wgridder.h#L1540 Yes so as long as they are correctly passed in it should be fine - the WCSlib dependencies can live on the exterior. The way to test that targetted faceting works is to take a high declination dataset take a off-axis source, predict it with a DFT using, e.g. Meqtrees, and target the source position. Cheers,

…

On Mon, Mar 28, 2022 at 11:40 AM mreineck ***@***.***> wrote: I agree completely that the l,m coordinates are not radians, they are Cartesian coordinates, and l*2+m*2=1 at the horizon. Perhaps it is easier to describe what I mean using code: https://github.com/mreineck/ducc/blob/1548bad07fac2f4331b551d7e239401d572a88e0/python/test/test_wgridder.py#L196 shows the unit test that makes sure our faceting gridder is doing what we expect it to do, even for very wide fields (i.e. no matter whether we grid a certain field in one go or in 2x2 or 2x4 tiles, we get the same result within the specified error ranges). If your interpretation of delta_l and delta_m is different, would it be possible to translate it into ours before calling the gridder? I'd rather not introduce a dependency on WCS for just this purpose (having no external dependencies feels really good :-) — Reply to this email directly, view it on GitHub <#353 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB4RE6R7AMFST2JQLDXRZ3DVCF5BFANCNFSM4LEIUZMQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- -- Benjamin Hugo PhD. student, Centre for Radio Astronomy Techniques and Technologies Department of Physics and Electronics Rhodes University Junior software developer Radio Astronomy Research Group South African Radio Astronomy Observatory Black River Business Park Observatory Cape Town

[mreineck]

Perfect! I don't have enough knowledge to set up the proposed test, but I'll be happy to hear about the results and help out with any problems that might turn up.