IP-2b: small correctness fixes (F12, F23, D3, 585#18)

changelog.d/+ip2b-small-fixes.default.md

@@ -0,0 +1 @@
+`MultipoleField` (with `r_min_mas=0`): the taper that makes the field single-valued at the origin is now `tanh(r/length_r)^|m|` instead of `tanh(r/spacing)^|m|`, so it no longer changes with `n_r`. With the default `length_r` (five node spacings) the region where orders m >= 1 are suppressed grows from about one to about five spacings, so prior draws and every existing `MultipoleField` fit change. Also: `MultipoleDisk.suggest_n_r` emits a `UserWarning` when it steps down below the nodes it wants (return value unchanged); `ToeplitzVisData(max_offset_mas=0.0)` gives a zero bound instead of the nuisance default; already-projected `OIData` observables with NaN, infinite or non-positive errors now raise `ValueError`.

docs/api/models/notes.md

@@ -394,10 +394,15 @@ evidence proper.
 single-valued at $r=0$. A smooth 2D field has Fourier coefficients
 $a_m(r)\sim r^{|m|}$ (the Zernike radial polynomials, Born & Wolf, *Principles
 of Optics*, section 9.2), so order $m$ of $\eta$ is multiplied by
-$\tanh(r/\Delta r)^{|m|}$, with $\Delta r$ the node spacing. The factor is
-$1 - O(2|m|\,e^{-2r/\Delta r})$, so the prior is stationary away from the
-centre (the variance of order $m$ is within $4|m|e^{-10}$ of $\sigma^2$ beyond
-five spacings); the latents stay whitened and only $m=0$ is free at $r=0$. It is
+$\tanh(r/\ell_r)^{|m|}$, with $\ell_r$ the radial correlation length, the
+scale on which the field varies (a physical length, so the prior does not
+depend on the grid). The variance of order $m$ is $\tanh(r/\ell_r)^{2|m|}$ of
+$\sigma^2$: it is within $4|m|e^{-10}$ of $\sigma^2$ beyond five correlation
+lengths, but high orders recover slowly (for $n_\theta = 64$, order 32 keeps
+about 10 per cent of its variance at $r = 2\ell_r$). The default
+$\ell_r$ is five node spacings, so the region where the prior is not
+stationary is about five spacings, not one, and a fitted $\ell_r$ also sets
+how much $m \ge 1$ structure the inner disk can carry; the latents stay whitened and only $m=0$ is free at $r=0$. It is
 applied only when `r_min_mas = 0`; an annulus keeps the plain prior.
 
 **Geometry.** By default the geometry parameters belong to the field: they are

docs/multipole_disk_rings.md

@@ -426,7 +426,7 @@ rms error of ln I0 over 2.0-29.0 mas (truth within e^-6 of its peak): GP 0.33, p
 
 ## The multipole field prior
 
-A `MultipoleField` is the same Gaussian-process idea on a polar grid in the disk's frame, here inclined 50° at position angle 30°. It is a residual that follows the disk's geometry. The top row shows three prior draws. The bottom shows the field's azimuthal Fourier coefficients: the cosine amplitude $a_m(r)$ of each harmonic $m = 0, \dots, 3$ in $\eta = \sum_m a_m \cos m\theta + b_m \sin m\theta$, against disk-plane radius, each on its own offset track (the $m = 0$ coefficient is the azimuthal mean). Every harmonic is a smooth function of radius with correlation length $\ell_r$. At the centre the $m \geq 1$ harmonics are tapered to zero, as $\tanh^m(r/\Delta r)$ with $\Delta r$ the radial node spacing, so that the field has a single value at the origin; this is why those tracks pinch to zero at $r = 0$. Higher harmonics have smaller amplitudes, with a scale set by $\ell_\theta$.
+A `MultipoleField` is the same Gaussian-process idea on a polar grid in the disk's frame, here inclined 50° at position angle 30°. It is a residual that follows the disk's geometry. The top row shows three prior draws. The bottom shows the field's azimuthal Fourier coefficients: the cosine amplitude $a_m(r)$ of each harmonic $m = 0, \dots, 3$ in $\eta = \sum_m a_m \cos m\theta + b_m \sin m\theta$, against disk-plane radius, each on its own offset track (the $m = 0$ coefficient is the azimuthal mean). Every harmonic is a smooth function of radius with correlation length $\ell_r$. At the centre the $m \geq 1$ harmonics are tapered to zero, as $\tanh^m(r/\ell_r)$, so that the field has a single value at the origin; this is why those tracks pinch to zero at $r = 0$. Higher harmonics have smaller amplitudes, with a scale set by $\ell_\theta$.
 
 ```python
 import jax

notebooks/multipole_disk_rings.ipynb

@@ -600,7 +600,7 @@
    "source": [
     "## The multipole field prior\n",
     "\n",
-    "A `MultipoleField` is the same Gaussian-process idea on a polar grid in the disk's frame, here inclined 50° at position angle 30°. It is a residual that follows the disk's geometry. The top row shows three prior draws. The bottom shows the field's azimuthal Fourier coefficients: the cosine amplitude $a_m(r)$ of each harmonic $m = 0, \\dots, 3$ in $\\eta = \\sum_m a_m \\cos m\\theta + b_m \\sin m\\theta$, against disk-plane radius, each on its own offset track (the $m = 0$ coefficient is the azimuthal mean). Every harmonic is a smooth function of radius with correlation length $\\ell_r$. At the centre the $m \\geq 1$ harmonics are tapered to zero, as $\\tanh^m(r/\\Delta r)$ with $\\Delta r$ the radial node spacing, so that the field has a single value at the origin; this is why those tracks pinch to zero at $r = 0$. Higher harmonics have smaller amplitudes, with a scale set by $\\ell_\\theta$."
+    "A `MultipoleField` is the same Gaussian-process idea on a polar grid in the disk's frame, here inclined 50° at position angle 30°. It is a residual that follows the disk's geometry. The top row shows three prior draws. The bottom shows the field's azimuthal Fourier coefficients: the cosine amplitude $a_m(r)$ of each harmonic $m = 0, \\dots, 3$ in $\\eta = \\sum_m a_m \\cos m\\theta + b_m \\sin m\\theta$, against disk-plane radius, each on its own offset track (the $m = 0$ coefficient is the azimuthal mean). Every harmonic is a smooth function of radius with correlation length $\\ell_r$. At the centre the $m \\geq 1$ harmonics are tapered to zero, as $\\tanh^m(r/\\ell_r)$, so that the field has a single value at the origin; this is why those tracks pinch to zero at $r = 0$. Higher harmonics have smaller amplitudes, with a scale set by $\\ell_\\theta$."
    ]
   },
   {

src/virgil/fields.py

@@ -573,13 +573,16 @@ def polar_values(self):
         )
         if self.r_min_mas == 0:
             # A smooth single-valued field has a_m(r) ~ r^|m| at the origin
-            # (Born & Wolf 1999, sec. 9.2): scale order m by tanh(r/spacing)^|m|,
-            # which is 1 - O(e^(-2r/spacing)), so the prior stays stationary
-            # away from the centre and only m = 0 survives at r = 0.
+            # (Born & Wolf 1999, sec. 9.2). The power law holds out to the
+            # scale on which the field varies, its radial correlation length
+            # ell (a physical length: the node spacing is only the grid's, so
+            # the taper would change with n_r). Scale order m by
+            # tanh(r/ell)^|m|, which is 1 - O(e^(-2r/ell)), so the prior stays
+            # stationary away from the centre and only m = 0 survives at r = 0.
             n_theta = field.shape[1]
             j = np.arange(n_theta)
             m = np.minimum(j, n_theta - j)
-            taper = np.tanh(nodes / self.spacing)[:, None] ** m[None, :]
+            taper = np.tanh(nodes / self.length_r)[:, None] ** m[None, :]
             field = np.real(
                 np.fft.ifft(taper * np.fft.fft(field, axis=1), axis=1)
             )

src/virgil/models.py

@@ -24,6 +24,7 @@
 import dataclasses
 import math
 import textwrap
+import warnings
 from typing import Any
 
 import equinox as eqx
@@ -53,7 +54,7 @@
     undo_elliptical_transf_spat_freq,
 )
 from . import _elr
-from ._deprecate import renamed
+from ._deprecate import external_stacklevel, renamed
 from .harmonics import cone_coefficients, rim_coefficients
 from ._utils import (
     FWHM_PER_SIGMA,
@@ -4744,8 +4745,8 @@ def suggest_n_r(self, data):
         phase of a ``height``'s lift (at the current parameter values, which
         must be concrete), rounded up to a
         multiple of 8 (but at most the 144 nodes buildable for a Gauss-Laguerre
-        tail). A starting point only: check convergence by
-        doubling ``n_r``.
+        tail), with a warning when it has to step down. A starting point
+        only: check convergence by doubling ``n_r``.
 
         Parameters
         ----------
@@ -4768,12 +4769,20 @@ def suggest_n_r(self, data):
         n = int(8 * onp.ceil(max(n, 16) / 8))
         # A Gauss-Laguerre tail (ExpTail, alone or inside a RadialSum/Gap)
         # has at most 150 nodes; step down in 8s until the rule builds.
+        wanted = n
         while n > 16:
             try:
                 self.profile.with_rule(n)
                 break
             except ValueError:
                 n -= 8
+        if n < wanted:
+            warnings.warn(
+                f"suggest_n_r wants {wanted} radial nodes but the profile "
+                f"cannot build that many; using {n}: the Hankel quadrature "
+                "may be under-resolved on the longest baselines.",
+                stacklevel=external_stacklevel(),
+            )
         return n
 
     def _alpha(self, r):

src/virgil/oidata.py

@@ -2027,14 +2027,22 @@ def _check_bandwidth(bandwidth):
 def _good_samples(values, errors, flag, n_samples, name, operator=False):
     """Mask of unflagged, usable samples, or ``None`` if all are good.
 
-    Already-projected data (size not ``n_samples``) are not checked.
+    Already-projected data (size not ``n_samples``) cannot be flagged per
+    sample, so they must be finite with positive finite errors.
     """
     if values.size != n_samples:
         if flag is not None:
             raise ValueError(
                 f"{name}_flag was given, but {name} has {values.size} values "
                 f"for {n_samples} samples (it looks already projected)."
             )
+        if not onp.all(_valid_samples(values, errors)):
+            raise ValueError(
+                f"{name} has {values.size} values for {n_samples} samples (it "
+                "looks already projected), and some are not finite or have a "
+                "non-positive or non-finite error: projected data cannot be "
+                "flagged, so remove them first."
+            )
         return None
     if operator:
         # A raw error of zero is usable: only the projected errors enter, and

src/virgil/radio_toeplitz.py

@@ -263,7 +263,11 @@ def __init__(
         bound = 0.0
         if free:
             # The largest image-frame shift (pixels) of the offsets' prior box.
-            shift_mas = max_offset_mas or nuisance._max_offset_mas
+            shift_mas = (
+                nuisance._max_offset_mas
+                if max_offset_mas is None
+                else max_offset_mas
+            )
             if shift_mas is None:
                 free_obs = onp.arange(nuisance.n_obs) != nuisance.ref_obs
                 unit = nuisance.offset_unit_mas[free_obs]
@@ -308,9 +312,6 @@ def __init__(
             last_l += [False] * (n_chunks - 1) + [True] * min(n_chunks, 1)
         idx = onp.concatenate(idx_l).reshape(-1, chunk)
         valid = onp.concatenate(valid_l).reshape(-1, chunk)
-        chunk_group = onp.asarray(group_l, dtype=onp.int32)
-        last = onp.asarray(last_l, dtype=bool)
-        nf_b = grid.grid_size(2 * npix, kp)
         nf_d = grid.grid_size(npix + 2 * math.ceil(bound), kp)
         margin = grid.margin(nf_d, kp) if free else None
         if free and tau > margin:
@@ -327,30 +328,26 @@ def chunked(x):
         self.plan = dict(
             idx=np.asarray(used[idx]),
             valid=np.asarray(valid),
-            group=np.asarray(chunk_group),
-            last=np.asarray(last),
+            group=np.asarray(group_l, dtype=onp.int32),
+            last=np.asarray(last_l, dtype=bool),
             tx=chunked(tx),
             ty=chunked(ty),
             h=chunked(h),
             rot=chunked(rot),
         )
-        self.data = data
-        self.mask = mask
-        self.npix = npix
+        self.data, self.mask, self.npix = data, mask, npix
         self.pixel_scale_mas = float(pix)
         self.rotation_deg = float(image.rotation_deg)
         self.blur_fwhm_mas = image.blur_fwhm_mas
         self.image_offset_mas = (dra, ddec)
-        self.free_offsets = free
-        self.max_shift_pix = float(bound)
-        self.nyquist_fraction = tau
-        self.margin = margin
+        self.free_offsets, self.max_shift_pix = free, float(bound)
+        self.nyquist_fraction, self.margin = tau, margin
         self.signature = _signature(nuisance)
         self.kp = kp
-        self.sizes = (n_groups, nf_b, nf_d)
-        count = onp.bincount(group, minlength=n_groups).astype(onp.float64)
+        self.sizes = (n_groups, grid.grid_size(2 * npix, kp), nf_d)
         self.arrays = dict(
-            self._compute(data.vis, data.weight), count=np.asarray(count)
+            self._compute(data.vis, data.weight),
+            count=np.asarray(counts, dtype=onp.float64),
         )
 
     def _compute(self, vis, weight):

tests/test_durations.json

@@ -79,6 +79,7 @@
  "tests/test_north_angle.py": 11.8,
  "tests/test_nufft.py": 46.7,
  "tests/test_observables.py": 28.2,
+ "tests/test_oidata_projected.py": 0.5,
  "tests/test_oifits.py": 14.9,
  "tests/test_orbit_cells.py": 280.0,
  "tests/test_orbit_coordinates.py": 12.2,

tests/test_multipole_disk.py

@@ -577,3 +577,13 @@ def test_a_profile_has_no_weight_and_a_sum_term_carries_it():
     assert onp.allclose(_cvis(heavy, u, v), _cvis(only_b, u, v), atol=1e-4)
     with pytest.raises(TypeError, match="not a profile"):
         RadialSum(a=(1.0, GaussianRing(3.0, 0.5)))
+
+
+def test_suggest_n_r_warns_when_it_steps_down():
+    class _Data:
+        u = onp.array([57.0])
+        v = onp.array([0.0])
+        wavel = onp.array([1e-6])
+
+    with pytest.warns(UserWarning, match="suggest_n_r wants"):
+        MultipoleDisk(ExpTail(0.0, 10.0)).suggest_n_r(_Data)

tests/test_multipole_field.py

@@ -82,13 +82,13 @@ def test_covariance_is_the_radial_matern_times_the_multipolekernel_spectrum(
     k_r = onp.exp(-0.5 * d * d) if onp.isinf(nu) else _matern(d, 1.0, nu)
     expected = sigma**2 * onp.kron(k_r, _azimuthal_covariance(N_T, ell_t, nu))
     # The prior is the plain separable one with order m of node k scaled by
-    # tanh(r_k / spacing)^|m|: T_k is the circulant with that symbol.
+    # tanh(r_k / ell_r)^|m|: T_k is the circulant with that symbol.
     mm = onp.minimum(onp.arange(N_T), N_T - onp.arange(N_T))
     taper = onp.stack(
         [
             onp.real(
                 onp.fft.ifft(
-                    onp.tanh(rk / field.spacing) ** mm[:, None]
+                    onp.tanh(rk / ell_r) ** mm[:, None]
                     * onp.fft.fft(onp.eye(N_T), axis=0),
                     axis=0,
                 )
@@ -130,9 +130,10 @@ def test_the_azimuthal_spectrum_is_multipolekernels(nu_theta):
 
 
 def test_latents_are_whitened_and_zero_latent_is_the_template():
-    # An empirical check that standard-normal latents give variance sigma^2.
+    # An empirical check that standard-normal latents give variance sigma^2
+    # (an annulus: the origin taper of a disk lowers the variance near r = 0).
     sigma = 2.0
-    field = _field(sigma=sigma)
+    field = _field(sigma=sigma, r_min_mas=2.0)
     draws = onp.random.default_rng(1).normal(size=(400, N_R, N_T))
     values = jax.vmap(lambda z: field.set("latent", z).polar_values())(draws)
     assert float(np.std(values)) == pytest.approx(sigma, rel=0.1)
@@ -789,14 +790,20 @@ def test_the_multipole_kernel_spectrum_is_finite_for_a_long_angle():
 def test_the_field_is_single_valued_at_the_centre_and_stationary_away_from_it():
     n_r, n_theta = 12, 16
     z = onp.random.default_rng(0).normal(size=(600, n_r, n_theta))
-    kw = dict(image_shape=(8, 8), r_max_mas=11.0, r_min_mas=0.0, sigma=2.0)
+    kw = dict(
+        image_shape=(8, 8),
+        r_max_mas=11.0,
+        r_min_mas=0.0,
+        sigma=2.0,
+        length_r=1.0,
+    )
     vals = onp.stack(
         [onp.asarray(MultipoleField(zi, **kw).polar_values()) for zi in z]
     )
     # r = 0: no azimuthal variation.
     assert onp.abs(vals[:, 0, :] - vals[:, 0, :1]).max() < 1e-5
-    # Per-order variance: gone at node 0 for m >= 1, tanh(r/spacing)^(2m)
-    # at node 1, and unchanged eight spacings out (node 11 is the outermost
+    # Per-order variance: gone at node 0 for m >= 1, tanh(r/ell)^(2m)
+    # at node 1 (ell = one spacing), and unchanged eight spacings out (node 11 is the outermost
     # node, node 8 is eight spacings out).
     var = (onp.abs(onp.fft.fft(vals, axis=2)) ** 2).mean(axis=0) / n_theta
     assert var[0, 3] < 1e-8
@@ -819,3 +826,24 @@ def test_an_annulus_keeps_the_plain_prior():
         )
     )
     assert onp.allclose(ring.polar_values(), expected, atol=1e-5)
+
+
+def test_the_origin_taper_is_set_by_the_correlation_length_not_the_grid():
+    # The same physical field (same ell_r, same r_max) sampled on two radial
+    # grids: the order-m variance at a given radius follows tanh(r/ell)^2m,
+    # whatever n_r is. Nodes at r = 1 exist on both grids (spacings 1, 0.5); there main's
+    # spacing taper gives 0.19 (n_r = 9) and 0.80 (n_r = 17) against 0.0098.
+    ell, m, n_theta = 2.0, 3, 16
+    var = {}
+    for n_r in (9, 17):
+        z = onp.random.default_rng(0).normal(size=(800, n_r, n_theta))
+        kw = dict(image_shape=(8, 8), r_max_mas=8.0, sigma=1.0, length_r=ell)
+        vals = onp.stack(
+            [onp.asarray(MultipoleField(zi, **kw).polar_values()) for zi in z]
+        )
+        v = (onp.abs(onp.fft.fft(vals, axis=2)) ** 2).mean(axis=0) / n_theta
+        var[n_r] = (v[(n_r - 1) // 8, m], v[-1, m])  # r = 1 and r = 8
+    ratios = [a / b for a, b in var.values()]
+    expected = (onp.tanh(1.0 / ell) / onp.tanh(8.0 / ell)) ** (2 * m)
+    for r in ratios:
+        assert r == pytest.approx(expected, rel=0.2)

tests/test_oidata_projected.py

@@ -0,0 +1,35 @@
+import numpy as onp
+import pytest
+
+from virgil.oidata import OIData
+
+
+def _record(vis, d_vis):
+    n, m = 8, 5
+    rng = onp.random.default_rng(0)
+    return {
+        "u": rng.uniform(1.0, 5.0, n),
+        "v": rng.uniform(1.0, 5.0, n),
+        "wavel": onp.full(n, 2.2e-6),
+        "vis": vis,
+        "d_vis": d_vis,
+        "vis_mat": rng.normal(size=(m, n)),
+    }
+
+
+def test_projected_data_build_when_finite():
+    OIData(_record(onp.full(5, 0.9), onp.full(5, 0.01)))
+
+
+@pytest.mark.parametrize(
+    "vis, d_vis",
+    [
+        (onp.array([0.9, onp.nan, 0.9, 0.9, 0.9]), onp.full(5, 0.01)),
+        (onp.full(5, 0.9), onp.array([0.01, onp.inf, 0.01, 0.01, 0.01])),
+        (onp.full(5, 0.9), onp.array([0.01, -0.01, 0.01, 0.01, 0.01])),
+        (onp.full(5, 0.9), onp.array([0.01, 0.0, 0.01, 0.01, 0.01])),
+    ],
+)
+def test_projected_data_must_be_finite_with_positive_errors(vis, d_vis):
+    with pytest.raises(ValueError, match="already projected"):
+        OIData(_record(vis, d_vis))

tests/test_radio_toeplitz.py

@@ -1074,3 +1074,13 @@ def both(sigma):
     assert t1 == pytest.approx(d1, rel=1e-12)
     assert t2 == pytest.approx(d2, rel=1e-12)
     assert t1 - t2 == pytest.approx(d1 - d2, rel=1e-10)
+
+
+def test_a_zero_max_offset_is_honoured_not_replaced_by_the_nuisance_bound(
+    small,
+):
+    data, nuis, scene, template = small
+    zero = data.toeplitz(template, nuis, offsets="spread", max_offset_mas=0.0)
+    default = data.toeplitz(template, nuis, offsets="spread")
+    assert zero.max_shift_pix == 0.0
+    assert default.max_shift_pix > 0.0