Make software compatible with Visibilities class

changelog/67.feature.rst

@@ -0,0 +1 @@
+Make the software compatible with :class:`~xrayvision.visibility.Visibilities`.

examples/stix.py

@@ -26,7 +26,7 @@
 stix_data = pickle.load(urllib.request.urlopen("https://pub099.cs.technik.fhnw.ch/demo/stix_vis.pkl"))
 
 time_range, energy_range, offset, stix_vis = stix_data
-stix_vis.phase_centre = [0, 0] * apu.arcsec
+stix_vis.phase_center = [0, 0] * apu.arcsec
 stix_vis.offset = offset
 
 ###############################################################################

xrayvision/clean.py

@@ -21,7 +21,7 @@
 from sunpy.map.map_factory import Map
 
 from xrayvision.imaging import vis_psf_image, vis_to_map
-from xrayvision.visibility import Visibility
+from xrayvision.visibility import Visibilities
 
 __all__ = ["clean", "vis_clean", "ms_clean", "vis_ms_clean"]
 
@@ -98,7 +98,7 @@ def clean(
     #     raise ValueError('')
     pad = [0 if x % 2 == 0 else 1 for x in dirty_map.shape]
 
-    # Assume beam, map phase_centre is in middle
+    # Assume beam, map phase_center is in middle
     beam_center = (dirty_beam.shape[0] - 1) / 2.0, (dirty_beam.shape[1] - 1) / 2.0
     map_center = (dirty_map.shape[0] - 1) / 2.0, (dirty_map.shape[1] - 1) / 2.0
 
@@ -173,7 +173,7 @@ def clean(
 
 @u.quantity_input
 def vis_clean(
-    vis: Visibility,
+    vis: Visibilities,
     shape: Quantity[u.pix],
     pixel_size: Quantity[u.arcsec / u.pix],
     clean_beam_width: Optional[Quantity[u.arcsec]] = 4.0,
@@ -403,7 +403,7 @@ def ms_clean(
 
 
 def vis_ms_clean(
-    vis: Visibility,
+    vis: Visibilities,
     shape: Quantity[u.pix],
     pixel_size: Quantity[u.arcsec / u.pix],
     scales: Optional[Iterable],

xrayvision/imaging.py

@@ -6,7 +6,7 @@
 from sunpy.map import GenericMap, Map
 
 from xrayvision.transform import dft_map, idft_map
-from xrayvision.visibility import Visibility
+from xrayvision.visibility import Visibilities
 
 __all__ = [
     "get_weights",
@@ -24,7 +24,7 @@
 WEIGHT_SCHEMES = ("natural", "uniform")
 
 
-def get_weights(vis: Visibility, scheme: str = "natural", norm: bool = True) -> np.ndarray:
+def get_weights(vis: Visibilities, scheme: str = "natural", norm: bool = True) -> np.ndarray:
     r"""
     Return spatial frequency weight factors for each visibility.
 
@@ -47,7 +47,7 @@ def get_weights(vis: Visibility, scheme: str = "natural", norm: bool = True) ->
         raise ValueError(f"Invalid weighting scheme {scheme}, must be one of: {WEIGHT_SCHEMES}")
     weights = np.sqrt(vis.u**2 + vis.v**2).value
     if scheme == "natural":
-        weights = np.ones_like(vis.vis, dtype=float)
+        weights = np.ones_like(vis.visibilities, dtype=float)
 
     if norm:
         weights /= weights.sum()
@@ -97,11 +97,11 @@ def image_to_vis(
     *,
     u: Quantity[apu.arcsec**-1],
     v: Quantity[apu.arcsec**-1],
-    phase_centre: Optional[Quantity[apu.arcsec]] = (0.0, 0.0) * apu.arcsec,
+    phase_center: Optional[Quantity[apu.arcsec]] = (0.0, 0.0) * apu.arcsec,
     pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1.0 * apu.arcsec / apu.pix,
-) -> Visibility:
+) -> Visibilities:
     r"""
-    Return a Visibility created from the image and u, v sampling.
+    Return a Visibilities object created from the image and u, v sampling.
 
     Parameters
     ----------
@@ -111,27 +111,27 @@ def image_to_vis(
         Array of u coordinates where the visibilities will be evaluated
     v :
         Array of v coordinates where the visibilities will be evaluated
-    phase_centre :
-        The coordinates the phase_centre.
+    phase_center :
+        The coordinates the phase_center.
     pixel_size :
         Size of pixels, if only one value is passed, assume square pixels (repeating the value).
 
     Returns
     -------
     :
-        The new visibility object
+        The new Visibilities object
 
     """
     pixel_size = validate_and_expand_kwarg(pixel_size, "pixel_size")
     if not (apu.get_physical_type((1 / u).unit) == ANGLE and apu.get_physical_type((1 / v).unit) == ANGLE):
         raise ValueError("u and v must be inverse angle (e.g. 1/deg or 1/arcsec")
-    vis = dft_map(image, u=u, v=v, phase_centre=phase_centre, pixel_size=pixel_size)
-    return Visibility(vis, u=u, v=v, offset=phase_centre)
+    vis = dft_map(image, u=u, v=v, phase_center=phase_center, pixel_size=pixel_size)
+    return Visibilities(vis, u=u, v=v, phase_center=phase_center)
 
 
 @apu.quantity_input()
 def vis_to_image(
-    vis: Visibility,
+    vis: Visibilities,
     shape: Quantity[apu.pix] = (65, 65) * apu.pixel,
     pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1 * apu.arcsec / apu.pix,
     scheme: str = "natural",
@@ -161,15 +161,21 @@ def vis_to_image(
     shape = shape.to(apu.pixel)
     weights = get_weights(vis, scheme=scheme)
     bp_arr = idft_map(
-        vis.vis, u=vis.u, v=vis.v, shape=shape, weights=weights, pixel_size=pixel_size, phase_centre=vis.phase_centre
+        vis.visibilities,
+        u=vis.u,
+        v=vis.v,
+        shape=shape,
+        weights=weights,
+        pixel_size=pixel_size,
+        phase_center=vis.phase_center,
     )
 
     return bp_arr
 
 
 @apu.quantity_input
 def vis_psf_map(
-    vis: Visibility,
+    vis: Visibilities,
     *,
     shape: Quantity[apu.pix] = (65, 65) * apu.pixel,
     pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1 * apu.arcsec / apu.pix,
@@ -203,7 +209,7 @@ def vis_psf_map(
 
 @apu.quantity_input()
 def vis_psf_image(
-    vis: Visibility,
+    vis: Visibilities,
     *,
     shape: Quantity[apu.pix] = (65, 65) * apu.pixel,
     pixel_size: Quantity[apu.arcsec / apu.pix] = 1 * apu.arcsec / apu.pix,
@@ -237,14 +243,19 @@ def vis_psf_image(
     # Make sure psf is always odd so power is in exactly one pixel
     shape = [s // 2 * 2 + 1 for s in shape.to_value(apu.pix)] * shape.unit
     psf_arr = idft_map(
-        np.ones(vis.vis.shape) * vis.vis.unit, u=vis.u, v=vis.v, shape=shape, weights=weights, pixel_size=pixel_size
+        np.ones(vis.visibilities.shape) * vis.visibilities.unit,
+        u=vis.u,
+        v=vis.v,
+        shape=shape,
+        weights=weights,
+        pixel_size=pixel_size,
     )
     return psf_arr
 
 
 @apu.quantity_input()
 def vis_to_map(
-    vis: Visibility,
+    vis: Visibilities,
     shape: Quantity[apu.pix] = (65, 65) * apu.pixel,
     pixel_size: Optional[Quantity[apu.arcsec / apu.pix]] = 1 * apu.arcsec / apu.pixel,
     scheme: Optional[str] = "natural",
@@ -278,7 +289,7 @@ def vis_to_map(
 
 
 @apu.quantity_input()
-def generate_header(vis: Visibility, *, shape: Quantity[apu.pix], pixel_size: Quantity[apu.arcsec / apu.pix]) -> dict:
+def generate_header(vis: Visibilities, *, shape: Quantity[apu.pix], pixel_size: Quantity[apu.arcsec / apu.pix]) -> dict:
     r"""
     Generate a map head given the visibilities, pixel size and shape
 
@@ -296,8 +307,8 @@ def generate_header(vis: Visibility, *, shape: Quantity[apu.pix], pixel_size: Qu
     :
     """
     header = {
-        "crval1": (vis.offset[1]).to_value(apu.arcsec),
-        "crval2": (vis.offset[0]).to_value(apu.arcsec),
+        "crval1": (vis.phase_center[1]).to_value(apu.arcsec),
+        "crval2": (vis.phase_center[0]).to_value(apu.arcsec),
         "cdelt1": (pixel_size[1] * apu.pix).to_value(apu.arcsec),
         "cdelt2": (pixel_size[0] * apu.pix).to_value(apu.arcsec),
         "ctype1": "HPLN-TAN",
@@ -312,9 +323,9 @@ def generate_header(vis: Visibility, *, shape: Quantity[apu.pix], pixel_size: Qu
 
 
 @apu.quantity_input()
-def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 / apu.arcsec]) -> Visibility:
+def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 / apu.arcsec]) -> Visibilities:
     r"""
-    Return a Visibility object created from the map, sampling it at give `u`, `v` coordinates.
+    Return a Visibilities object created from the map, sampling it at give `u`, `v` coordinates.
 
     Parameters
     ----------
@@ -328,7 +339,7 @@ def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 /
     Returns
     -------
     :
-        The new visibility object
+        The new Visibilities object
 
     """
     if not apu.get_physical_type(1 / u) == ANGLE and apu.get_physical_type(1 / v) == ANGLE:
@@ -346,6 +357,8 @@ def map_to_vis(amap: GenericMap, *, u: Quantity[1 / apu.arcsec], v: Quantity[1 /
     if "cdelt2" in meta:
         new_psize[0] = float(meta["cdelt2"])
 
-    vis = image_to_vis(amap.data, u=u, v=v, pixel_size=new_psize * apu.arcsec / apu.pix)
-    vis.offset = new_pos * apu.arcsec
+    vis = image_to_vis(
+        amap.quantity, u=u, v=v, pixel_size=new_psize * apu.arcsec / apu.pix, phase_center=new_pos * apu.arcsec
+    )
+
     return vis

xrayvision/mem.py

@@ -27,7 +27,7 @@
     "mem",
 ]
 
-from xrayvision.visibility import Visibility
+from xrayvision.visibility import Visibilities
 
 logger = get_logger(__name__, "DEBUG")
 
@@ -81,8 +81,8 @@ def _get_fourier_matrix(vis, shape=(64, 64) * apu.pix, pixel_size=(4.0312500, 4.
     The complex Fourier matrix
     """
     m, n = shape.to("pix")
-    y = generate_xy(m, phase_centre=0 * apu.arcsec, pixel_size=pixel_size[1])
-    x = generate_xy(n, phase_centre=0 * apu.arcsec, pixel_size=pixel_size[0])
+    y = generate_xy(m, phase_center=0 * apu.arcsec, pixel_size=pixel_size[1])
+    x = generate_xy(n, phase_center=0 * apu.arcsec, pixel_size=pixel_size[0])
     x, y = np.meshgrid(x, y, indexing="ij")
     uv = np.vstack([vis.u, vis.v])
     # Check apu are correct for exp need to be dimensionless and then remove apu for speed
@@ -192,8 +192,8 @@ def _prepare_for_optimise(pixel, shape, vis):
     # 'Hv' is composed by the union of its real and imaginary part
     Hv = np.concatenate([ReHv, ImHv], axis=-1)
     # Division of real and imaginary part of the visibilities
-    ReV = np.real(vis.vis)
-    ImV = np.imag(vis.vis)
+    ReV = np.real(vis.visibilities)
+    ImV = np.imag(vis.visibilities)
     # 'Visib' is composed by the real and imaginary part of the visibilities
     Visib = np.concatenate([ReV, ImV], axis=-1)
     # Standard deviation of the real and imaginary part of the visibilities
@@ -256,7 +256,7 @@ def _get_mean_visibilities(vis, shape, pixel):
     v = np.zeros(n_vis) * (1 / apu.arcsec)
     den = np.zeros(n_vis)
     weights = np.zeros_like(vis.amplitude_error**2)
-    visib = np.zeros_like(vis.vis)
+    visib = np.zeros_like(vis.visibilities)
     for ip in range(n_vis):
         # what about 0.5 pix offset
         i = ru[ip].to_value("pix").astype(int)
@@ -269,15 +269,15 @@ def _get_mean_visibilities(vis, shape, pixel):
             # we save in 'u' and 'v' the u and v coordinates of the first frequency that corresponds
             # to the position (i, j) of the discretization of the (u,v)-plane 'iuarr'
 
-            visib[count] = vis.vis[ip]
+            visib[count] = vis.visibilities[ip]
             weights[count] = vis.amplitude_error[ip] ** 2.0
             den[count] = 1.0
             iuarr[i, j] = count
 
             count += 1
         else:
             # Save the sum of the visibilities that correspond to the same position (i, j)
-            visib[iuarr[i, j].astype(int)] += vis.vis[ip]
+            visib[iuarr[i, j].astype(int)] += vis.visibilities[ip]
             # Save the number of the visibilities that correspond to the same position (i, j)
             den[iuarr[i, j].astype(int)] += 1.0
             # Save the sum of the variances of the amplitudes of the visibilities that
@@ -296,7 +296,7 @@ def _get_mean_visibilities(vis, shape, pixel):
     # correspond to the same position in the discretization of the (u,v)-plane
     weights = np.sqrt(weights[:count]) / den
 
-    return SimpleNamespace(u=u, v=v, vis=visib, amplitude_error=weights)
+    return SimpleNamespace(u=u, v=v, visibilities=visib, amplitude_error=weights)
 
 
 def _proximal_entropy(y, m, lamba, Lip, tol=10**-10):
@@ -601,7 +601,7 @@ def _get_percent_lambda(vis):
     if ibig.size < 2:
         snr_value = -1
     else:
-        snr_value, _ = resistant_mean((np.abs(vis.vis[ibig]) / vis.amplitude_error[ibig]).flatten(), 3)
+        snr_value, _ = resistant_mean((np.abs(vis.visibilities[ibig]) / vis.amplitude_error[ibig]).flatten(), 3)
 
     # TODO magic numbers
     percent_lambda = 2 / (snr_value**2 + 90)
@@ -611,7 +611,7 @@ def _get_percent_lambda(vis):
 
 @apu.quantity_input
 def mem(
-    vis: Visibility,
+    vis: Visibilities,
     shape: Quantity[apu.pix],
     pixel_size: Quantity[apu.arcsec / apu.pix],
     *,