{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# 1D RM-synthesis weighting\n", "\n", "RM-synthesis is a Fourier transform from $\\lambda^2$ to Faraday depth $\\phi$, so\n", "the RMSF (the response function, analogous to the interferometric dirty beam) is\n", "the Fourier transform of the per-channel *weights* over $\\lambda^2$. The channel\n", "sampling in $\\lambda^2$ is fixed by the band, but the weights are ours to choose,\n", "exactly like uv weighting in aperture-synthesis imaging.\n", "\n", "`rm-lite` provides the interferometric weighting schemes:\n", "\n", "| weight_type | definition | interferometric analogue |\n", "| --- | --- | --- |\n", "| `variance` / `natural` | $1/\\sigma^2$ | natural (max sensitivity) |\n", "| `uniform` | equal per channel | (no weighting) |\n", "| `uniform_lsq` | equal per $\\lambda^2$ cell | uniform (max resolution) |\n", "| `briggs` | robust interpolation | Briggs / robust |\n", "\n", "`briggs` spans the range with a single `robust` knob (large positive $\\to$\n", "`natural`, large negative $\\to$ `uniform_lsq`).\n", "\n", "This notebook has two parts. **Part A** builds the $\\lambda^2$-grid weighting\n", "step by step on ASKAP **RACS** coverage (sub-bands separated by gaps, so the grid\n", "and its occupancy are easy to see). **Part B** compares every scheme on\n", "contiguous **GMIMS-DRAGONS** coverage, where the weighting types differ most.\n", "Every claim below is checked with an `assert`." ] }, { "cell_type": "code", "execution_count": null, "id": "1", "metadata": {}, "outputs": [], "source": [ "from __future__ import annotations\n", "\n", "import logging\n", "\n", "import astropy.units as u\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import polars as pl\n", "from astropy.visualization import quantity_support\n", "from numpy.typing import NDArray\n", "from rm_lite.tools_1d import rmclean, rmsynth\n", "from rm_lite.utils.logging import quiet_logs\n", "from rm_lite.utils.simulate import faraday_slab_spectrum\n", "from rm_lite.utils.synthesis import (\n", " FDFOptions,\n", " briggs_weight,\n", " compute_rmsynth_params,\n", " freq_to_lambda2,\n", " inverse_rmsynth_nufft,\n", " uniform_lsq_weight,\n", ")\n", "\n", "plt.rcParams[\"figure.dpi\"] = 150\n", "_ = quantity_support()\n", "rng = np.random.default_rng(42)\n", "\n", "# rm_lite logs each synth/clean run; `with quiet_logs(logging.ERROR):` around the\n", "# runs below keeps the notebook quiet.\n", "\n", "# Distinct styles so overlapping traces (e.g. variance and uniform) stay visible.\n", "LINE_STYLES = [\"-\", \"--\", \"-.\", \":\", (0, (3, 1, 1, 1)), (0, (5, 2))]\n", "MARKERS = [\".\", \"x\", \"+\", \"^\", \"s\", \"d\"]" ] }, { "cell_type": "markdown", "id": "2", "metadata": {}, "source": [ "A Burn slab gives the Faraday-thick source; a helper measures the RMSF FWHM." ] }, { "cell_type": "code", "execution_count": null, "id": "3", "metadata": {}, "outputs": [], "source": [ "def fwhm_from_curve(phi: NDArray[np.float64], amp: NDArray[np.float64]) -> float:\n", " \"\"\"FWHM of the main lobe from the half-maximum crossings (linear interp).\"\"\"\n", " amp = amp / amp.max()\n", " peak = int(np.argmax(amp))\n", "\n", " def cross(idxs: range) -> float:\n", " prev = peak\n", " for i in idxs:\n", " if amp[i] < 0.5:\n", " frac = (0.5 - amp[prev]) / (amp[i] - amp[prev])\n", " return float(phi[prev] + frac * (phi[i] - phi[prev]))\n", " prev = i\n", " return float(phi[idxs[-1]])\n", "\n", " return cross(range(peak + 1, len(amp))) - cross(range(peak - 1, -1, -1))" ] }, { "cell_type": "markdown", "id": "4", "metadata": {}, "source": [ "## Part A - building the grid weighting (RACS-all)\n", "\n", "ASKAP RACS observes in separated sub-bands (low/mid/high). That gapped coverage\n", "is ideal for showing how `uniform_lsq`/`briggs` are built: lay a virtual\n", "$\\lambda^2$ grid, count the per-cell occupancy, and weight each channel by the\n", "inverse local density. (Channel counts are downsampled here for speed.)" ] }, { "cell_type": "code", "execution_count": null, "id": "5", "metadata": {}, "outputs": [], "source": [ "# ASKAP RACS-all: low/mid/high sub-bands (channel counts downsampled for speed).\n", "bw_low, bw_mid, bw_high = 288, 144, 288\n", "low = np.linspace(943.5 - bw_low / 2, 943.5 + bw_low / 2, 36) * u.MHz\n", "mid = np.linspace(1367.5 - bw_mid / 2, 1367.5 + bw_mid / 2, 9) * u.MHz\n", "high = np.linspace(1655.5 - bw_high / 2, 1655.5 + bw_high / 2, 9) * u.MHz\n", "freq_racs = np.concatenate([low, mid, high]).to(u.Hz).value\n", "lam2_racs = freq_to_lambda2(freq_racs)\n", "\n", "# Sampling in both frequency (top) and the lambda^2 it maps to (bottom): three\n", "# sub-bands separated by gaps, and lambda^2 bunches up at the high-frequency end.\n", "fig, (axf, axl) = plt.subplots(2, 1, figsize=(8, 4), sharex=True)\n", "axf.plot(freq_racs / 1e6, np.ones_like(freq_racs), \"|\", ms=18)\n", "axf.set(yticks=[], ylabel=\"channels\")\n", "axl.plot(freq_racs / 1e6, lam2_racs, \".\")\n", "axl.set(xlabel=\"frequency [MHz]\", ylabel=r\"$\\lambda^2$ [m$^2$]\")\n", "axf.set_title(\"RACS-all coverage: three sub-bands, dense in $\\\\lambda^2$ at high freq\")\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "id": "6", "metadata": {}, "source": [ "### The virtual $\\lambda^2$ grid\n", "\n", "The resolution cell is $\\mathrm{cell} = \\sqrt{3}/\\phi_\\mathrm{max}$ (the RMSF\n", "cell). Each channel is snapped to a cell index; the gaps between sub-bands leave\n", "empty cells." ] }, { "cell_type": "code", "execution_count": null, "id": "7", "metadata": {}, "outputs": [], "source": [ "phi_max = 300.0\n", "cell = np.sqrt(3.0) / phi_max\n", "cell_idx = np.floor((lam2_racs - lam2_racs.min()) / cell).astype(int)\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(lam2_racs, cell_idx, \".\")\n", "ax.set(\n", " xlabel=r\"$\\lambda^2$ [m$^2$]\",\n", " ylabel=\"virtual cell index\",\n", " title=\"each channel snapped to a virtual $\\\\lambda^2$ cell\",\n", ")" ] }, { "cell_type": "markdown", "id": "8", "metadata": {}, "source": [ "### Occupancy $\\to$ weights\n", "\n", "`uniform_lsq` weight $= $ natural weight $/$ local density, where the density is\n", "the total natural weight in a channel's cell over the cell width. So each\n", "*occupied* cell contributes equally, and channels sharing a cell get one weight.\n", "`briggs` blends this with the natural weight via the `robust` knob. With flat\n", "noise the natural weight is 1, so `uniform_lsq` reduces to $1/N_\\mathrm{cell}$." ] }, { "cell_type": "code", "execution_count": null, "id": "9", "metadata": {}, "outputs": [], "source": [ "natural = np.ones_like(lam2_racs) # flat noise -> natural weight = 1\n", "count = np.bincount(cell_idx) # channels per cell\n", "\n", "w_uniform_lsq = uniform_lsq_weight(lam2_racs, natural, cell)\n", "w_briggs0 = briggs_weight(lam2_racs, natural, 0.0, cell)\n", "\n", "fig, (axo, axw) = plt.subplots(1, 2, figsize=(12, 4))\n", "axo.step(np.arange(len(count)), count, where=\"mid\")\n", "axo.set(\n", " xlabel=\"cell index\",\n", " ylabel=\"channels in cell\",\n", " title=\"per-cell occupancy (empty cells = gaps)\",\n", ")\n", "for w, lbl, mk in [\n", " (natural, \"natural (= uniform, flat)\", \".\"),\n", " (w_uniform_lsq, r\"uniform_lsq\", \"o\"),\n", " (w_briggs0, r\"briggs ($r=0$)\", \"x\"),\n", "]:\n", " axw.plot(lam2_racs, w / np.nansum(w), mk, ms=5, ls=\"none\", label=lbl)\n", "axw.set(\n", " xlabel=r\"$\\lambda^2$ [m$^2$]\", ylabel=\"normalised weight\", title=\"resulting weights\"\n", ")\n", "axw.legend(fontsize=\"small\")\n", "fig.tight_layout()\n", "\n", "# Checks: channels in a cell share one weight (no single-channel jump within a\n", "# cell), and each occupied cell contributes equal total weight.\n", "for c in np.unique(cell_idx):\n", " in_cell = w_uniform_lsq[cell_idx == c]\n", " np.testing.assert_allclose(in_cell, in_cell[0], rtol=1e-12)\n", "totals = np.array([w_uniform_lsq[cell_idx == c].sum() for c in np.unique(cell_idx)])\n", "np.testing.assert_allclose(totals, totals[0], rtol=1e-12)\n", "print(\"within-cell uniform + each occupied cell contributes equally: OK\")" ] }, { "cell_type": "markdown", "id": "10", "metadata": {}, "source": [ "## Part B - comparison on GMIMS-DRAGONS\n", "\n", "Now switch to contiguous **GMIMS-DRAGONS** coverage (350-1030 MHz). The wide\n", "low-frequency band gives the weighting real leverage, so `uniform_lsq` narrows\n", "the RMSF substantially. We inject a Faraday-thick Burn slab (thickness\n", "10 rad/m$^2$ at 30 rad/m$^2$), add noise, and compare the schemes." ] }, { "cell_type": "code", "execution_count": null, "id": "11", "metadata": {}, "outputs": [], "source": [ "freq_hz = (np.arange(350, 1030, 1) * u.MHz).to(u.Hz).value\n", "lambda_sq = freq_to_lambda2(freq_hz)\n", "\n", "rm_radm2 = 30.0\n", "delta_rm_radm2 = 10.0\n", "frac_pol = 0.5\n", "psi0_deg = 10.0\n", "model_kw = {\n", " \"frac_pol\": frac_pol,\n", " \"psi0_deg\": psi0_deg,\n", " \"rm_radm2\": rm_radm2,\n", " \"delta_rm_radm2\": delta_rm_radm2,\n", "}\n", "complex_model = faraday_slab_spectrum(lambda_sq, **model_kw)\n", "rms_noise = 0.02\n", "complex_noisy = (\n", " complex_model\n", " + rng.normal(0, rms_noise, freq_hz.size)\n", " + 1j * rng.normal(0, rms_noise, freq_hz.size)\n", ").astype(np.complex128)\n", "complex_err = np.ones_like(complex_noisy) * (rms_noise + rms_noise * 1j)\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(freq_hz, complex_noisy.real, \".\", color=\"tab:red\", alpha=0.5, label=\"Stokes Q\")\n", "ax.plot(freq_hz, complex_noisy.imag, \".\", color=\"tab:blue\", alpha=0.5, label=\"Stokes U\")\n", "ax.plot(freq_hz, complex_model.real, color=\"tab:red\", label=\"model Q\")\n", "ax.plot(freq_hz, complex_model.imag, color=\"tab:blue\", label=\"model U\")\n", "ax.legend(ncol=2)\n", "ax.set(xlabel=rf\"$\\nu$ / {u.Hz:latex_inline}\", ylabel=\"Fractional flux\")" ] }, { "cell_type": "markdown", "id": "12", "metadata": {}, "source": [ "A helper runs RM-synthesis for one weighting and returns the weights, RMSF and\n", "FDF. `compute_rmsynth_params` exposes the actual per-channel weight array." ] }, { "cell_type": "code", "execution_count": null, "id": "13", "metadata": {}, "outputs": [], "source": [ "def run_weighting(**kw: object) -> dict:\n", " \"\"\"Run RM-synthesis for one weighting; collect weights, RMSF, FDF, synth.\"\"\"\n", " opts = FDFOptions(n_samples=50, phi_max_radm2=300.0, **kw) # type: ignore[arg-type]\n", " params = compute_rmsynth_params(freq_hz, complex_noisy, complex_err, opts)\n", " synth = rmsynth.run_rmsynth(\n", " freq_arr_hz=freq_hz,\n", " complex_pol_arr=complex_noisy,\n", " complex_pol_error=complex_err,\n", " n_samples=50,\n", " phi_max_radm2=300.0,\n", " **kw, # type: ignore[arg-type]\n", " )\n", " fdf_arrs, rmsf_arrs = synth[1], synth[2]\n", " phi2 = rmsf_arrs[\"phi2_arr_radm2\"].to_numpy()\n", " rmsf = np.abs(rmsf_arrs[\"rmsf_complex_arr\"].to_numpy().astype(complex))\n", " return {\n", " \"synth\": synth,\n", " \"weight\": params.weight_arr,\n", " \"phi\": fdf_arrs[\"phi_arr_radm2\"].to_numpy(),\n", " \"fdf\": np.abs(fdf_arrs[\"fdf_dirty_complex_arr\"].to_numpy().astype(complex)),\n", " \"phi2\": phi2,\n", " \"rmsf\": rmsf,\n", " \"fwhm\": fwhm_from_curve(phi2, rmsf),\n", " }" ] }, { "cell_type": "markdown", "id": "14", "metadata": {}, "source": [ "## Run every weighting\n", "\n", "We run all the schemes once and reuse them in every plot below: `variance` and\n", "`uniform` (which coincide here because the noise is flat), `uniform_lsq`, and\n", "`briggs` at `robust` $= -1, 0, +1$." ] }, { "cell_type": "code", "execution_count": null, "id": "15", "metadata": {}, "outputs": [], "source": [ "weightings = {\n", " \"variance\": {\"weight_type\": \"variance\"},\n", " \"uniform\": {\"weight_type\": \"uniform\"},\n", " \"uniform_lsq\": {\"weight_type\": \"uniform_lsq\"},\n", " \"briggs r=-1\": {\"weight_type\": \"briggs\", \"robust\": -1.0},\n", " \"briggs r=0\": {\"weight_type\": \"briggs\", \"robust\": 0.0},\n", " \"briggs r=+1\": {\"weight_type\": \"briggs\", \"robust\": 1.0},\n", "}\n", "# TeX display labels (keys above stay plain for indexing/tables).\n", "tex_label = {\n", " \"variance\": \"variance\",\n", " \"uniform\": \"uniform\",\n", " \"uniform_lsq\": \"uniform_lsq\",\n", " \"briggs r=-1\": r\"Briggs ($r=-1$)\",\n", " \"briggs r=0\": r\"Briggs ($r=0$)\",\n", " \"briggs r=+1\": r\"Briggs ($r=+1$)\",\n", "}\n", "with quiet_logs(logging.ERROR):\n", " results = {label: run_weighting(**kw) for label, kw in weightings.items()}" ] }, { "cell_type": "code", "execution_count": null, "id": "16", "metadata": {}, "outputs": [], "source": [ "# Weights vs lambda^2: uniform_lsq / low-robust up-weight the sparse large-lambda^2 end.\n", "fig, ax = plt.subplots()\n", "for (label, res), marker in zip(results.items(), MARKERS, strict=False):\n", " w = res[\"weight\"] / np.nansum(res[\"weight\"])\n", " ax.plot(lambda_sq, w, marker, ms=3, ls=\"none\", label=tex_label[label])\n", "ax.set(\n", " xlabel=rf\"$\\lambda^2$ / {u.m**2:latex_inline}\",\n", " ylabel=\"normalised weight\",\n", " title=r\"Weights vs $\\lambda^2$\",\n", ")\n", "ax.legend(ncol=2, fontsize=\"small\")" ] }, { "cell_type": "markdown", "id": "17", "metadata": {}, "source": [ "### RMSF\n", "\n", "`uniform_lsq` (and low-robust `briggs`) narrow the RMSF main lobe; `variance`,\n", "`uniform` and high-robust `briggs` are the widest, most sensitive response. Over\n", "this wide contiguous band the effect is sizeable: `uniform_lsq` narrows the FWHM\n", "by roughly a third. The `briggs` points run monotonically between the endpoints\n", "(`natural` widest, `uniform_lsq` narrowest)." ] }, { "cell_type": "code", "execution_count": null, "id": "18", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "for (label, res), ls in zip(results.items(), LINE_STYLES, strict=False):\n", " ax.plot(res[\"phi2\"], res[\"rmsf\"] / res[\"rmsf\"].max(), ls=ls, label=tex_label[label])\n", "ax.set(\n", " xlim=(-25, 25),\n", " xlabel=rf\"$\\phi$ / {u.rad / u.m**2:latex_inline}\",\n", " ylabel=\"RMSF (normalised)\",\n", " title=\"RMSF main lobe by weighting\",\n", ")\n", "ax.legend(ncol=2, fontsize=\"small\")\n", "\n", "fwhm = {label: res[\"fwhm\"] for label, res in results.items()}\n", "assert abs(fwhm[\"variance\"] - fwhm[\"uniform\"]) < 0.5, \"flat noise: natural == uniform\"\n", "narrowing = (fwhm[\"uniform\"] - fwhm[\"uniform_lsq\"]) / fwhm[\"uniform\"]\n", "assert narrowing > 0.2, \"uniform_lsq narrows the RMSF by roughly a third\"\n", "# endpoints bracket the sweep: natural widest, uniform_lsq narrowest, briggs\n", "# monotonic between them (tol absorbs sub-cell FWHM measurement jitter).\n", "tol = 0.1\n", "assert fwhm[\"briggs r=+1\"] <= fwhm[\"variance\"] + tol\n", "assert fwhm[\"briggs r=-1\"] <= fwhm[\"briggs r=+1\"] + tol\n", "assert fwhm[\"uniform_lsq\"] <= fwhm[\"briggs r=-1\"] + tol\n", "\n", "pl.DataFrame({\"weighting\": list(fwhm), \"rmsf_fwhm_radm2\": list(fwhm.values())})" ] }, { "cell_type": "markdown", "id": "19", "metadata": {}, "source": [ "### Dirty FDF vs the model slab\n", "\n", "The true Faraday dispersion of a Burn slab is a top-hat over\n", "$[\\mathrm{RM}-\\Delta/2,\\ \\mathrm{RM}+\\Delta/2]$. The dirty FDF is that top-hat\n", "blurred by the RMSF, so the narrower-RMSF weightings give sharper slab edges." ] }, { "cell_type": "code", "execution_count": null, "id": "20", "metadata": {}, "outputs": [], "source": [ "phi_ref = results[\"variance\"][\"phi\"]\n", "model_tophat = (np.abs(phi_ref - rm_radm2) <= delta_rm_radm2 / 2).astype(float)\n", "\n", "fig, ax = plt.subplots()\n", "for (label, res), ls in zip(results.items(), LINE_STYLES, strict=False):\n", " ax.plot(res[\"phi\"], res[\"fdf\"] / res[\"fdf\"].max(), ls=ls, label=tex_label[label])\n", "ax.fill_between(phi_ref, model_tophat, color=\"k\", alpha=0.1, label=\"model slab\")\n", "ax.set(\n", " xlim=(0, 60),\n", " xlabel=rf\"$\\phi$ / {u.rad / u.m**2:latex_inline}\",\n", " ylabel=\"|FDF| (normalised)\",\n", " title=\"Dirty FDF vs the model slab\",\n", ")\n", "ax.legend(ncol=2, fontsize=\"small\")" ] }, { "cell_type": "markdown", "id": "21", "metadata": {}, "source": [ "### RM-CLEAN: recovered flux and the resolution/sensitivity trade\n", "\n", "For a symmetric slab the **flux-weighted mean $\\phi$** locates the source, and the\n", "debiased zeroth moment `mom0_debias` measures the **recovered polarised flux**.\n", "Every weighting recovers the slab centre, but the high-resolution weightings\n", "recover *less* faint flux at a fixed threshold: resolution costs sensitivity." ] }, { "cell_type": "code", "execution_count": null, "id": "22", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "recovered_flux = {}\n", "mean_rm = {}\n", "cleaned = {}\n", "with quiet_logs(logging.ERROR):\n", " for (label, res), ls in zip(results.items(), LINE_STYLES, strict=False):\n", " clean = rmclean.run_rmclean_from_synth(\n", " rm_synth_1d_results=res[\"synth\"],\n", " auto_mask=8,\n", " auto_threshold=1.0,\n", " multiscale=True,\n", " )\n", " cleaned[label] = clean\n", " arrs = clean.fdf_arrs\n", " phi = arrs[\"phi_arr_radm2\"].to_numpy()\n", " clean_fdf = np.abs(arrs[\"fdf_clean_complex_arr\"].to_numpy().astype(complex))\n", " model = np.abs(arrs[\"fdf_model_complex_arr\"].to_numpy().astype(complex))\n", " mean_rm[label] = float(np.sum(phi * model) / np.sum(model))\n", " recovered_flux[label] = float(clean.fdf_parameters[\"mom0_debias\"][0])\n", " assert abs(mean_rm[label] - rm_radm2) < 5.0, \"CLEAN recovers the slab centre\"\n", " ax.plot(phi, clean_fdf / clean_fdf.max(), ls=ls, label=tex_label[label])\n", "ax.fill_between(phi, model_tophat, color=\"k\", alpha=0.1, label=\"model slab\")\n", "ax.set(\n", " xlim=(0, 60),\n", " xlabel=rf\"$\\phi$ / {u.rad / u.m**2:latex_inline}\",\n", " ylabel=\"|clean FDF| (normalised)\",\n", " title=\"Cleaned FDF vs the model slab\",\n", ")\n", "ax.legend(ncol=2, fontsize=\"small\")\n", "\n", "assert recovered_flux[\"variance\"] >= recovered_flux[\"uniform_lsq\"]\n", "\n", "pl.DataFrame(\n", " {\n", " \"weighting\": list(recovered_flux),\n", " \"mean_rm_radm2\": list(mean_rm.values()),\n", " \"mom0_recovered\": list(recovered_flux.values()),\n", " }\n", ")" ] }, { "cell_type": "markdown", "id": "23", "metadata": {}, "source": [ "The recovered `mom0` compared with the true input polarised fraction\n", "($\\mathrm{mom0}_\\mathrm{true} = f_\\mathrm{pol} = 0.5$, the integrated\n", "polarisation of the slab). Every weighting under-recovers the diffuse thick\n", "source at this threshold, and the higher-resolution weightings recover the\n", "least." ] }, { "cell_type": "code", "execution_count": null, "id": "24", "metadata": {}, "outputs": [], "source": [ "flux_table = pl.DataFrame(\n", " {\n", " \"weighting\": list(recovered_flux),\n", " \"mom0_recovered\": list(recovered_flux.values()),\n", " \"mom0_true\": frac_pol,\n", " }\n", ").with_columns(\n", " (pl.col(\"mom0_recovered\") / pl.col(\"mom0_true\")).alias(\"fraction_recovered\")\n", ")\n", "flux_table" ] }, { "cell_type": "markdown", "id": "25", "metadata": {}, "source": [ "### Depolarisation vs slab width and the Faraday max-scale\n", "\n", "RM-synthesis is blind to Faraday structure broader than the **max-scale**,\n", "$\\phi_\\mathrm{max\\text{-}scale} \\approx \\pi / \\lambda^2_\\mathrm{min}$ (set by\n", "the shortest wavelength, i.e. the top of the band). A Burn slab thicker than this\n", "is resolved out and depolarises. Each slab is normalised to the **same total\n", "polarised flux** (the integral over Faraday depth), so its height scales as\n", "$f_\\mathrm{pol}/\\Delta$; that total is $|P(\\lambda^2{=}0)| = f_\\mathrm{pol}$\n", "in this `sinc` model, independent of width. We sweep `robust` from $-2$ to $+2$\n", "and plot the recovered debiased `mom0` for each width. Wider slabs sit lower\n", "(more depolarised); the dashed line is the input flux." ] }, { "cell_type": "code", "execution_count": null, "id": "26", "metadata": {}, "outputs": [], "source": [ "faraday_max_scale = float(np.pi / lambda_sq.min())\n", "print(f\"Faraday max-scale ~ {faraday_max_scale:.0f} rad/m^2\")\n", "\n", "\n", "def recover_mom0(delta_rm_radm2: float, robust: float) -> float:\n", " \"\"\"Recovered debiased mom0 for a slab of given width and briggs robust.\"\"\"\n", " model = faraday_slab_spectrum(\n", " lambda_sq, frac_pol, psi0_deg, rm_radm2, delta_rm_radm2\n", " )\n", " noisy = (\n", " model\n", " + rng.normal(0, rms_noise, freq_hz.size)\n", " + 1j * rng.normal(0, rms_noise, freq_hz.size)\n", " ).astype(np.complex128)\n", " err = np.ones_like(noisy) * (rms_noise + rms_noise * 1j)\n", " synth = rmsynth.run_rmsynth(\n", " freq_arr_hz=freq_hz,\n", " complex_pol_arr=noisy,\n", " complex_pol_error=err,\n", " n_samples=10,\n", " phi_max_radm2=400.0,\n", " weight_type=\"briggs\",\n", " robust=robust,\n", " )\n", " clean = rmclean.run_rmclean_from_synth(\n", " rm_synth_1d_results=synth,\n", " auto_mask=8,\n", " auto_threshold=1.0,\n", " )\n", " return float(clean.fdf_parameters[\"mom0_debias\"][0])\n", "\n", "\n", "robust_grid = np.arange(-2.0, 2.0 + 1e-6, 0.25)\n", "slab_widths = [\n", " 1,\n", " 5,\n", " 10.0,\n", " 30.0,\n", " 60.0,\n", " 120.0,\n", "]\n", "mom0_grid = {}\n", "with quiet_logs(logging.ERROR): # keep the many synth/clean runs quiet\n", " for width in slab_widths:\n", " mom0_grid[width] = [recover_mom0(width, r) for r in robust_grid]\n", "\n", "fig, ax = plt.subplots()\n", "for (width, mom0), ls in zip(mom0_grid.items(), LINE_STYLES, strict=False):\n", " ax.plot(\n", " robust_grid,\n", " mom0,\n", " ls=ls,\n", " marker=\".\",\n", " ms=4,\n", " label=rf\"$\\Delta$RM={width:.0f} ({width / faraday_max_scale:.2f} max-scale)\",\n", " )\n", "ax.axhline(frac_pol, color=\"k\", ls=\":\", label=\"input flux\")\n", "ax.set(\n", " xlabel=\"Briggs robust\",\n", " ylabel=\"recovered mom0 (debiased)\",\n", " title=\"Depolarisation vs slab width across the robust sweep\",\n", ")\n", "ax.legend(fontsize=\"small\", loc=\"upper left\")\n", "\n", "# Thin slabs are recovered; slabs wider than the max-scale depolarise away.\n", "assert max(mom0_grid[slab_widths[0]]) > 0.6 * frac_pol, \"thin slab recovered\"\n", "assert max(mom0_grid[slab_widths[-1]]) < 0.25 * frac_pol, \"wide slab depolarised\"\n", "print(\"Depolarisation checks passed.\")" ] }, { "cell_type": "code", "execution_count": null, "id": "27", "metadata": {}, "outputs": [], "source": [ "# RMSF resolution across the same robust sweep: FWHM moves smoothly from\n", "# uniform_lsq (large -robust) to natural (large +robust), flattening by |r|~2.\n", "with quiet_logs(logging.ERROR):\n", " robust_fwhm = [\n", " run_weighting(weight_type=\"briggs\", robust=r)[\"fwhm\"] for r in robust_grid\n", " ]\n", "\n", "fwhm_natural = results[\"variance\"][\"fwhm\"]\n", "fwhm_lsq = results[\"uniform_lsq\"][\"fwhm\"]\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(robust_grid, robust_fwhm, \".-\", label=\"Briggs\")\n", "ax.axhline(fwhm_natural, color=\"k\", ls=\":\", label=\"natural\")\n", "ax.axhline(fwhm_lsq, color=\"grey\", ls=\"--\", label=r\"uniform_lsq\")\n", "ax.set(\n", " xlabel=\"Briggs robust\",\n", " ylabel=rf\"RMSF FWHM / {u.rad / u.m**2:latex_inline}\",\n", " title=\"RMSF resolution across the robust sweep\",\n", ")\n", "ax.legend(fontsize=\"small\")\n", "\n", "# Higher robust -> wider RMSF (less resolution); +/-2 reach the endpoints.\n", "assert robust_fwhm[-1] >= robust_fwhm[0], \"higher robust is wider\"\n", "assert abs(robust_fwhm[0] - fwhm_lsq) < 0.5, \"r=-2 approximates uniform_lsq\"\n", "assert abs(robust_fwhm[-1] - fwhm_natural) < 0.5, \"r=+2 approximates natural\"\n", "print(\"RMSF resolution checks passed.\")" ] }, { "cell_type": "markdown", "id": "28", "metadata": {}, "source": [ "### Extrapolating each clean model back to $\\lambda^2 = 0$\n", "\n", "Inverse RM-synthesis turns each clean model FDF back into $Q/U$ over $\\lambda^2$.\n", "On a contiguous grid from $0$ (the intrinsic, unrotated polarisation) to the band\n", "maximum, comparing with the true model shows: **in band** the reconstruction\n", "tracks the model for every weighting (to within the noise); **at $\\lambda^2 =\n", "0$**, well outside the band, a Faraday-thick model is unconstrained and diverges\n", "from the truth." ] }, { "cell_type": "code", "execution_count": null, "id": "29", "metadata": {}, "outputs": [], "source": [ "# Contiguous lambda^2 grid from 0 to the band maximum (the observed samples\n", "# are irregular; a dense grid gives smooth reconstruction curves).\n", "lambda_sq_ext = np.linspace(0.0, lambda_sq.max(), 500)\n", "model_ext = faraday_slab_spectrum(lambda_sq_ext, **model_kw)\n", "in_band = (lambda_sq_ext >= lambda_sq.min()) & (lambda_sq_ext <= lambda_sq.max())\n", "\n", "fig, (ax0, ax1, ax2) = plt.subplots(3, 1, sharex=True, figsize=(6, 9))\n", "ax0.plot(lambda_sq_ext, model_ext.real, \"k\", lw=2, label=\"model Q\")\n", "ax1.plot(lambda_sq_ext, model_ext.imag, \"k\", lw=2, label=\"model U\")\n", "ax2.plot(lambda_sq_ext, np.abs(model_ext), \"k\", lw=2, label=\"model pI\")\n", "band_rms = {}\n", "with quiet_logs(logging.ERROR):\n", " for (label, clean), ls in zip(cleaned.items(), LINE_STYLES, strict=False):\n", " arrs = clean.fdf_arrs\n", " phi = arrs[\"phi_arr_radm2\"].to_numpy()\n", " fdf_model = arrs[\"fdf_model_complex_arr\"].to_numpy().astype(complex)\n", " lam0 = float(clean.fdf_parameters[\"lam_sq_0_m2\"][0])\n", " recon = inverse_rmsynth_nufft(fdf_model, lambda_sq_ext, phi, lam0)\n", " band_rms[label] = float(\n", " np.sqrt(np.mean(np.abs(recon[in_band] - model_ext[in_band]) ** 2))\n", " )\n", " ax0.plot(lambda_sq_ext, recon.real, ls=ls, label=tex_label[label])\n", " ax1.plot(lambda_sq_ext, recon.imag, ls=ls, label=tex_label[label])\n", " ax2.plot(lambda_sq_ext, np.abs(recon), ls=ls, label=tex_label[label])\n", " print(f\"{label:12s} in-band reconstruction RMS = {band_rms[label]:.4f}\")\n", "ax0.axvline(0.0, color=\"grey\", ls=\":\")\n", "ax1.axvline(0.0, color=\"grey\", ls=\":\")\n", "ax2.axvline(0.0, color=\"grey\", ls=\":\")\n", "\n", "for i, ax in enumerate((ax0, ax1, ax2)):\n", " ax.axvline(\n", " lambda_sq.min(),\n", " color=\"grey\",\n", " ls=\"--\",\n", " label=\"measured band\" if i == 0 else None,\n", " )\n", " ax.axvline(lambda_sq.max(), color=\"grey\", ls=\"--\")\n", "\n", "ax0.set(ylabel=\"Q\")\n", "ax1.set(ylabel=\"U\")\n", "ax2.set(xlabel=rf\"$\\lambda^2$ / {u.m**2:latex_inline}\", ylabel=\"pI\")\n", "ax0.legend(ncol=2, fontsize=\"small\")\n", "\n", "# Every weighting reconstructs the in-band spectrum to within the noise.\n", "for label, rms in band_rms.items():\n", " assert rms < 3 * rms_noise, f\"{label} in-band reconstruction diverges\"\n", "print(\"Extrapolation checks passed.\")" ] }, { "cell_type": "markdown", "id": "30", "metadata": {}, "source": [ "## Summary\n", "\n", "- `briggs` tunes continuously between `natural` (max sensitivity) and\n", " `uniform_lsq` (max resolution); narrower-RMSF weightings recover less faint flux\n", " and reconstruct the spectrum slightly worse.\n", "- `uniform_lsq` is interferometric **uniform** weighting on the $\\lambda^2$ grid:\n", " each channel is weighted by the inverse local density (per-cell natural-weight\n", " occupancy over the cell width), so each occupied cell contributes equally and\n", " channels sharing a cell get one weight. The weight steps where the true\n", " sampling density changes (gaps, channelisation) - correct inverse-density\n", " weighting, not aliasing. The resolution gain is coverage-dependent (large on\n", " the wide GMIMS band, above).\n", "- Extrapolating a clean model to $\\lambda^2 = 0$ is reliable only within the\n", " observed band; a Faraday-thick source is unconstrained beyond it.\n", "- **Weighting does not correct an uncorrected Stokes $I$ spectrum.** It reshapes\n", " the RMSF but cannot remove the multiplicative $I(\\lambda^2)$ that biases the\n", " effective RMSF; divide out Stokes $I$ (fractional polarisation) for that." ] } ], "metadata": { "kernelspec": { "display_name": "rm-lite (3.11.13)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.13" } }, "nbformat": 4, "nbformat_minor": 5 }