Add on the fly rephasing and beam interpolation/reprojection to hci worker (#139)

* don't depend on tart2ms

* switch hci worker to ray

* == instead of > typo in hci

* switch hci worker to ray. fix initial rephasing issues. add basic beam interpolation and reprojection code

* do not change the order/number of antennas during UVW recomputation

* add ant to index mapping in uvw computation

* allow implicit_stokes in degrid

* initial beam interp+reprojection functional

* add multi-ms support to hci worker + beam zarr outputs

* fix output names and pass new ra dec through to fits headers

* add fid to init outputs

* add separate beam interpolation function to perform both spatial and pol rotation

* move target compuattion to after rephasing + beam reproject

* check npix is even if passed in from cli

* fix mask typo

* rephase fiddle factors

* rephase like chgcentre

* fix beam orientation

* remove debugging statements

* add implicit_stokes=True when calling convert in degrid

* add frequency axis to zarr outputs of hci worker

* more generic ant to index mapping in synthesize_uvw

* depend on africanus master

* add stack option to hci to produce a single zarr cube instead of slices

* modify mean and flag computation to not rewrite all data vars

* add log messages for scaffold creation

* add log message for reduction to compute mean

* add stacked-ds as an hci output

* do not flag planes in hci (breifast also makes diagnostic plots)

* change wsum -> weight + remove flag and initialise nonzero in hci

* sync remote

* remove unintialised flag from outputs of hci

* single hci output

* bester -> landmanbester in dependabot config

* update beam-model info string

Author: landmanbester

.github/dependabot.yml

@@ -9,9 +9,9 @@ updates:
       time: "09:00"
     open-pull-requests-limit: 10
     reviewers:
-      - "@bester"
+      - "@landmanbester"
     assignees:
-      - "@bester"
+      - "@landmanbester"
     commit-message:
       prefix: "ci"
       include: "scope"

pfb/operators/gridder.py

@@ -335,7 +335,8 @@ def image_data_products(dsl,
                         do_residual=True,
                         do_weight=True,
                         do_noise=False,
-                        do_beam=False):
+                        do_beam=False,
+                        min_padding=1.7):
     '''
     Function to compute image space data products in one go
 
@@ -474,11 +475,20 @@ def image_data_products(dsl,
     if robustness is not None:
         numba_threads = np.maximum(nthreads, 1)
         numba.set_num_threads(numba_threads)
+        # we need to compute the weights on the padded grid
+        # but we don't have control over the optimal gridding
+        # parameters so assume a minimum
+        nx_pad = int(np.ceil(min_padding*nx))
+        if nx_pad%2:
+            nx_pad += 1
+        ny_pad = int(np.ceil(min_padding*ny))
+        if ny_pad%2:
+            ny_pad += 1
         counts = _compute_counts(uvw,
                                  freq,
                                  mask,
                                  wgt,
-                                 nx, ny,
+                                 nx_pad, ny_pad,
                                  cellx, celly,
                                  uvw.dtype,
                                  k=0,
@@ -492,7 +502,7 @@ def image_data_products(dsl,
             freq,
             wgt,
             mask,
-            nx, ny,
+            nx_pad, ny_pad,
             cellx, celly,
             robustness,
             usign=1.0 if flip_u else -1.0,
@@ -711,26 +721,23 @@ def compute_residual(dsl,
     # do not apply weights in this direction
     convim = np.zeros_like(dirty)
     for c in range(ncorr):
-        try:
-            model_vis = dirty2vis(
-                uvw=uvw,
-                freq=freq,
-                dirty=beam[c]*model[c],
-                pixsize_x=cellx,
-                pixsize_y=celly,
-                center_x=x0,
-                center_y=y0,
-                flip_u=flip_u,
-                flip_v=flip_v,
-                flip_w=flip_w,
-                epsilon=epsilon,
-                do_wgridding=do_wgridding,
-                nthreads=nthreads,
-                divide_by_n=False,  # incorporate in smooth beam
-                sigma_min=1.1, sigma_max=3.0,
-                verbosity=0)
-        except:
-            import ipdb; ipdb.set_trace()
+        model_vis = dirty2vis(
+            uvw=uvw,
+            freq=freq,
+            dirty=beam[c]*model[c],
+            pixsize_x=cellx,
+            pixsize_y=celly,
+            center_x=x0,
+            center_y=y0,
+            flip_u=flip_u,
+            flip_v=flip_v,
+            flip_w=flip_w,
+            epsilon=epsilon,
+            do_wgridding=do_wgridding,
+            nthreads=nthreads,
+            divide_by_n=False,  # incorporate in smooth beam
+            sigma_min=1.1, sigma_max=3.0,
+            verbosity=0)
 
         vis2dirty(
             uvw=uvw,

pfb/parser/hci.yaml

@@ -95,6 +95,10 @@ inputs:
       - double
     info:
       Gridding precision
+  beam-model:
+    dtype: str
+    info:
+      Path to beam model as an xarray dataset backed by zarr
   field-of-view:
     dtype: float
     abbreviation: fov
@@ -205,6 +209,11 @@ inputs:
     info:
       Rephase visibilities to this phase center.
       Should be a string in the format 'HH:MM:SS,DD:MM:SS' (note the , delimiter)
+  stack:
+    dtype: bool
+    default: false
+    info:
+      Stack everything into a single xarray dataset
 
   _include:
     - (.)gridding.yml
@@ -218,9 +227,6 @@ outputs:
     dtype: Directory
     must_exist: false
     mkdir: false
-  # optional-fits:
-  #   implicit: '{current.output-filename}_{current.product}.fits'
-  #   required: false
 
 policies:
   pass_missing_as_none: true