Hermitian pairs of visibilities for imaging and fitting

[iancze]

It would be worthwhile reviewing/auditing how we treat raw visibilities and their Hermitian pairs through the fitting and imaging processes.

We know that Hermitian pairs must be included for any routine that uses an iFFT to synthesize an image, otherwise the image will not turn out to be real.

However, it's less clear whether they should be included when we are using an image to forward model the visibilities and evaluating a likelihood on those, because there is the potential for "double counting" measurements. I don't think it makes a technical difference (since all measurements would be double counted) and would just be a scale factor adjustment relative to other regularizer settings. However, it would be good to come to a consensus on this and document it thoroughly.

[iancze]

The chi_square and log_likelihood routines should not accept Hermitian pairs, this would be equivalent to double-counting the data. This doesn't change the minimum of the loss function surface, but does change the parameter uncertainties that would be inferred.

It doesn't matter whether this happens in the "normalized log likelihood" functions, since the normalization over data points makes the calculation the same as whether it has Hermitian pairs or not. These normalized loss functions should not be used for inference, anyway.

[jeffjennings]

Note to potentially come back to: in #132 , the RandomCellSplitGridded class takes in a griddedDataset, which contains Hermitian pairs, but the selection/splitting of data in RandomCellSplitGridded doesn't maintain that Hermitian property.

[iancze]

I've come to the conclusion that it would be best that any routine that is evaluating a likelihood function (or a loss function) should not be involved with Hermitian pairs, whether they are "loose" or whether they are gridded.

If Hermitian pairs were to be included in the likelihood function calculation, then, in order for the calculation to be accurate, the simple \chi^2 sum of data - model would need to be replaced with a multi-variate Gaussian likelihood that included covariances between data points. There would be a perfect covariance between each point and its Hermitian pair.

The one case where Hermitian pairs do need to be included is in any routine that is creating an image using the inverse FFT. This applies to all of the routines involved in making dirty images. It does not apply to any forward-modeling routine involved in RML, like FourierCube or the NuFFT.

[iancze]

Merging #156 has gone a long way towards resolving this issue, but I think a little bit of work still remains regarding documentation. I think the loose-visibilities.md tutorial should be updated to be more specific about the difference in gridding operations performed by the DirtyImager vs. DataAverager, and how only the DataAverager is the correct one for inference. This tutorial should also have a discussion about the perfect covariance matrix between points and their Hermitian pairs, as described above.

[iancze]

#248 makes significant progress on this issue by consolidating docs in various places.

I discovered/remembered that mpol.datasets.DartBoard and its tests have numerous references to treating Hermitian pairs. This may have made sense when datasets.GriddedDataset contained Hermitian pairs, but after #156, we made the decision that gridding.DataAverager and datasets.GriddedDataset will not contain Hermitian pairs (they should only be present internally to gridding.DirtyImager).

The tests still seem to pass but there's a good change that Dartboard.build_grid_mask_from_cells contains incorrect logic.

mask = np.zeros_like(self.cartesian_qs, dtype="bool")

# uses about a Gb..., and this only 256x256
for cell_index in cell_index_list:
    qi, pi = cell_index
    q_min, q_max = self.q_edges[qi : qi + 2]
    p0_min, p0_max = self.phi_edges[pi : pi + 2]
    # also include Hermitian values
    p1_min, p1_max = self.phi_edges[pi : pi + 2] - np.pi

    # whether or not the q and phi values of the coordinate array
    # fit in the q cell and *either of* the regular or Hermitian phi cell
    ind = (
        (self.cartesian_qs >= q_min)
        & (self.cartesian_qs < q_max)
        & (
            ((self.cartesian_phis > p0_min) & (self.cartesian_phis <= p0_max))
            | ((self.cartesian_phis > p1_min) & (self.cartesian_phis <= p1_max))
        )
    )

    mask[ind] = True

return mask

And we should update the documentation and internal comments to be clear about Hermitian pairs in these routines.