{ "cells": [ { "cell_type": "markdown", "id": "81693eb7", "metadata": {}, "source": [ "# 12 — Minkowski Tensors\n", "\n", "**Purpose:** Quantify morphological *anisotropy* in CIB and tSZ patches using\n", "rank-2 Minkowski tensors — the tensorial generalisation of the scalar Minkowski\n", "functionals computed in `docs/05_plots.ipynb`.\n", "\n", "For each intensity threshold ν, the map is binarised to the excursion set\n", "K = {x : T(x) > ν} and a rank-2 tensor W is computed from the geometry of K.\n", "Eigendecomposing W yields:\n", "- **β = λ_min / λ_max ∈ [0, 1]**: anisotropy index (1 = isotropic, 0 = maximally elongated)\n", "- **θ ∈ (-π/2, π/2]**: orientation of the major axis\n", "\n", "Three tensor types are computed to expose their different sensitivities:\n", "- **W012** (W^{0,2}_1): boundary normal tensor — sums n⊗n over boundary pixels.\n", " Most sensitive to the *shape of cluster boundaries*. Recommended default.\n", "- **W200** (W^{2,0}_0): area inertia tensor — sums r⊗r over interior pixels.\n", " Measures *elongation of filled regions*.\n", "- **W201** (W^{2,0}_1): boundary position tensor — sums r⊗r over boundary pixels.\n", " Hybrid between W012 and W200.\n", "\n", "Comparing β(ν) and θ(ν) between Agora, DDPM, and Gaussian baseline tests\n", "whether the diffusion model reproduces the morphological anisotropy of the\n", "true foreground fields — a property that scalar statistics (power spectra,\n", "moments, MFs) are blind to.\n", "\n", "**Sanity check:** A statistically isotropic Gaussian random field should give\n", "β ≈ 1 on average and θ uniformly distributed — deviations indicate genuine\n", "non-Gaussian anisotropy.\n", "\n", "**Inputs:**\n", "- Test maps: `data/low_pass/2mJy/*.npy`\n", "- DDPM samples: `data/low_pass/2mJy/new_samples_*.npy`\n", "- Norm params: `data/low_pass/2mJy/norm_params_2mJy.npy`\n", "\n", "**Outputs:**\n", "- `plots/minkowski_tensors_beta.pdf` — β(ν) curves, all tensor types, CIB and tSZ\n", "- `plots/minkowski_tensors_theta.pdf` — θ polar histograms at fixed thresholds\n", "- `plots/minkowski_tensors_residuals.pdf` — (β_Agora − β_DDPM) / σ_Agora for W012\n", "\n", "**Key module functions:**\n", "- `foregrounds_diffusion.statistics.compute_minkowski_tensors`\n", "- `foregrounds_diffusion.statistics.MINKOWSKI_TENSOR_DESCRIPTIONS`\n", "\n", "**Reference:** Schroder-Turk et al. (2013), *New J. Phys.* 15 083028" ] }, { "cell_type": "markdown", "id": "a0e6793a", "metadata": {}, "source": [ "# Load modules and set configuration.\n", "# THRESHOLDS is in the normalised [0, 1] range (after apply_maxmin_normalization);\n", "# avoid ν = 0 and ν = 1 where the excursion set is trivially full or empty.\n", "## 1 Configuration and tensor descriptions\n", "\n", "Three tensor types are computed at 30 intensity thresholds spanning ν ∈ [0.05, 0.95]\n", "of the normalised range:\n", "\n", "| Tensor | Symbol | Physical interpretation |\n", "|--------|--------|------------------------|\n", "| W012 | W^{0,2}_1 | Boundary normal tensor: sums n⊗n over boundary pixels (Sobel normals). Most sensitive to the *shape* of cluster boundaries. |\n", "| W200 | W^{2,0}_0 | Area inertia tensor: sums r⊗r over interior pixels. Measures *elongation* of the filled excursion region. |\n", "| W201 | W^{2,0}_1 | Boundary position tensor: sums r⊗r over boundary pixels. Hybrid between W012 and W200. |\n", "\n", "The eigendecomposition W = λ₁v₁v₁ᵀ + λ₂v₂v₂ᵀ (λ₁ ≤ λ₂) yields:\n", "- **β = λ₁/λ₂ ∈ [0, 1]**: anisotropy index. β = 1 means the excursion set\n", " is statistically isotropic; β → 0 means it is highly elongated.\n", "- **θ ∈ (−π/2, π/2]**: orientation of the major axis." ] }, { "cell_type": "code", "execution_count": null, "id": "d92f7f43", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from pathlib import Path\n", "\n", "from foregrounds_diffusion.preprocessing import apply_maxmin_normalization\n", "from foregrounds_diffusion.morphology import (\n", " compute_minkowski_tensors,\n", " MINKOWSKI_TENSOR_DESCRIPTIONS,\n", ")\n", "\n", "PTSRC = 2\n", "PIXEL_RES_ARCMIN = 1.40625\n", "N_MAPS = 100 # maps per source; reduce for quick tests\n", "\n", "# 30 thresholds in the normalised [0, 1] range, avoiding the extreme tails\n", "# where excursion sets are near-empty or near-full\n", "THRESHOLDS = np.linspace(0.05, 0.95, 30)\n", "\n", "TENSOR_TYPES = ('W012', 'W200', 'W201')\n", "\n", "PATCHES_DIR = Path(f'data/low_pass/{PTSRC}mJy')\n", "PROJECT_ROOT = Path('/home/apb86/cmb_foregrounds_diffusion')\n", "\n", "print('Tensor descriptions:')\n", "for tt, desc in MINKOWSKI_TENSOR_DESCRIPTIONS.items():\n", " print(f' {tt}: {desc}')\n" ] }, { "cell_type": "markdown", "id": "326c01b9", "metadata": {}, "source": [ "# Load normalised patches. The Minkowski tensor code re-normalises internally\n", "# via norm_fn, so the absolute pixel values here are less critical than for\n", "# power spectra. Using the standard z-score / min-max loading ensures\n", "# consistency with the train/test split.\n", "## 2 Load and denormalise maps\n", "\n", "The same loading pattern as notebooks 06–11. For the Minkowski tensor\n", "computation the maps are *re-normalised to [0, 1]* by `apply_maxmin_normalization`\n", "inside `compute_minkowski_tensors`, regardless of the input units. This ensures\n", "that the threshold parameter ν is dimensionless and comparable across maps with\n", "very different amplitudes (e.g. tSZ clusters of different masses)." ] }, { "cell_type": "code", "execution_count": null, "id": "9833220c", "metadata": {}, "outputs": [], "source": [ "# ---------------------------------------------------------------------------\n", "# Load and denormalise maps (identical pattern to notebooks 10 and 11)\n", "# ---------------------------------------------------------------------------\n", "norm_params = np.load(PATCHES_DIR / f'norm_params_{PTSRC}mJy.npy')\n", "cib_mean, cib_std, tsz_mean, tsz_std = norm_params\n", "print(f'CIB — mean: {cib_mean:.4f} µK std: {cib_std:.4f} µK')\n", "print(f'tSZ — mean: {tsz_mean:.4f} µK std: {tsz_std:.4f} µK')\n", "\n", "cib_maps = np.load(PATCHES_DIR / f'CIB_map_150GHz_256_st6_zscore_{PTSRC}mJy_lp.npy')\n", "tsz_maps = np.load(PATCHES_DIR / f'tSZ3_map_150GHz_256_st6_zscore_{PTSRC}mJy_lp.npy')\n", "\n", "rng = np.random.default_rng(seed=42)\n", "indices = rng.permutation(len(cib_maps))\n", "test_idx = indices[int(0.8 * len(cib_maps)):]\n", "\n", "cib_test = cib_maps[test_idx, :, :, 0] * cib_std + cib_mean\n", "tsz_test = tsz_maps[test_idx, :, :, 0] * tsz_std + tsz_mean\n", "print(f'Test set: {len(cib_test)} patches')" ] }, { "cell_type": "markdown", "id": "46ec2a0b", "metadata": {}, "source": [ "# DDPM samples shape: (N, 2, H, W) channels-first.\n", "# Index [:, 0] = CIB channel, [:, 1] = tSZ channel.\n", "# Denormalise: CIB = raw * (cib_max - cib_min) + cib_min;\n", "# tSZ = raw * tsz_std + tsz_mean.\n", "## 3 Load DDPM samples and Gaussian baseline" ] }, { "cell_type": "code", "execution_count": null, "id": "b3ea2cb0", "metadata": {}, "outputs": [], "source": [ "ddpm_raw = np.load(\n", " PROJECT_ROOT / 'data' / 'low_pass' / f'{PTSRC}mJy' /\n", " f'new_samples_cib_tsz_{PTSRC}mJy_zero_norm_6x6_w_au_lp.npy'\n", ") # (N, 2, H, W)\n", "ddpm_cib = ddpm_raw[:, 0] * cib_std + cib_mean\n", "ddpm_tsz = ddpm_raw[:, 1] * tsz_std + tsz_mean\n", "print(f'DDPM samples: {len(ddpm_cib)} patches')\n", "\n", "gauss_maps = np.load(PATCHES_DIR / f'gaussian_cib_tsz_{PTSRC}mJy_lp.npy')\n", "gauss_cib = gauss_maps[:, 0] if gauss_maps.shape[1] == 2 else gauss_maps[:, :, :, 0]\n", "gauss_tsz = gauss_maps[:, 1] if gauss_maps.shape[1] == 2 else gauss_maps[:, :, :, 1]\n", "print(f'Gaussian samples: {len(gauss_cib)} patches')\n", "\n", "N = min(N_MAPS, len(cib_test), len(ddpm_cib), len(gauss_cib))\n", "print(f'Using {N} maps per source')" ] }, { "cell_type": "markdown", "id": "414ae9a4", "metadata": {}, "source": [ "# compute_minkowski_tensors(maps_nhw, norm_fn, thresholds, tensor_types, centred)\n", "# norm_fn : normalisation applied per-map before thresholding\n", "# (use apply_maxmin_normalization so ν ∈ [0,1] is meaningful)\n", "# centred : if True, subtract excursion-set centroid from r⊗r tensors\n", "# (affects W200, W201 only; default True)\n", "# n_jobs : pass n_jobs=-1 to use all CPU cores via joblib\n", "## 4 Compute Minkowski tensors\n", "\n", "`compute_minkowski_tensors` loops over N maps and T thresholds. For each\n", "(map, threshold) pair:\n", "1. Apply `norm_fn` (here `apply_maxmin_normalization`) to map the patch to [0,1].\n", "2. Binarise: `binary = normalised_map > ν`.\n", "3. Skip trivially empty (< 2 interior pixels) or trivially full (< 2 exterior\n", " pixels) excursion sets — these are uninformative and return β = 1, θ = 0.\n", "4. Compute the requested tensor(s) from the binary mask.\n", "5. Eigendecompose to get β, θ.\n", "\n", "Returns a dict:\n", "```python\n", "{\n", " 'W012': {'beta': ndarray(N, T), 'theta': ndarray(N, T)},\n", " 'W200': {'beta': ..., 'theta': ...},\n", " 'W201': {'beta': ..., 'theta': ...},\n", "}\n", "```\n", "Runtime: ≈ 5–15 min for 100 maps × 30 thresholds × 3 tensor types on CPU.\n", "Use `n_jobs=-1` to parallelise over maps with joblib." ] }, { "cell_type": "code", "execution_count": null, "id": "3a7330e6", "metadata": {}, "outputs": [], "source": [ "# ---------------------------------------------------------------------------\n", "# Compute Minkowski tensors\n", "# NOTE: ~5-15 min for 100 maps × 30 thresholds × 3 tensor types on CPU.\n", "# Reduce N_MAPS or THRESHOLDS for a quick test.\n", "# ---------------------------------------------------------------------------\n", "norm_fn = apply_maxmin_normalization\n", "\n", "print('Computing Agora CIB tensors...')\n", "mt_agora_cib = compute_minkowski_tensors(\n", " cib_test[:N], norm_fn, THRESHOLDS, tensor_types=TENSOR_TYPES)\n", "\n", "print('Computing Agora tSZ tensors...')\n", "mt_agora_tsz = compute_minkowski_tensors(\n", " tsz_test[:N], norm_fn, THRESHOLDS, tensor_types=TENSOR_TYPES)\n", "\n", "print('Computing DDPM CIB tensors...')\n", "mt_ddpm_cib = compute_minkowski_tensors(\n", " ddpm_cib[:N], norm_fn, THRESHOLDS, tensor_types=TENSOR_TYPES)\n", "\n", "print('Computing DDPM tSZ tensors...')\n", "mt_ddpm_tsz = compute_minkowski_tensors(\n", " ddpm_tsz[:N], norm_fn, THRESHOLDS, tensor_types=TENSOR_TYPES)\n", "\n", "print('Computing Gaussian CIB tensors...')\n", "mt_gauss_cib = compute_minkowski_tensors(\n", " gauss_cib[:N], norm_fn, THRESHOLDS, tensor_types=TENSOR_TYPES)\n", "\n", "print('Computing Gaussian tSZ tensors...')\n", "mt_gauss_tsz = compute_minkowski_tensors(\n", " gauss_tsz[:N], norm_fn, THRESHOLDS, tensor_types=TENSOR_TYPES)\n", "\n", "print('Done.')" ] }, { "cell_type": "markdown", "id": "f36e246b", "metadata": {}, "source": [ "## 5 Anisotropy index β(ν): comparing CIB and tSZ morphology\n", "\n", "Plot β(ν) — the mean anisotropy index as a function of threshold ν — in a\n", "2-row (CIB, tSZ) × 3-column (W012, W200, W201) layout. Error bands show\n", "the standard deviation across the N test maps.\n", "\n", "A statistically isotropic random field (e.g. the Gaussian baseline) has\n", "β ≈ 1 at all thresholds by symmetry. Non-isotropic structure — elongated\n", "tSZ filaments, asymmetric CIB source clusters — drives β < 1. The DDPM\n", "should reproduce the Agora β(ν) curve if it has learnt the correct cluster\n", "morphology; the Gaussian baseline provides a reference for how much\n", "anisotropy is generated purely by the patch geometry." ] }, { "cell_type": "code", "execution_count": null, "id": "8413f716", "metadata": {}, "outputs": [], "source": [ "# ---------------------------------------------------------------------------\n", "# Plot β(ν): anisotropy index vs threshold\n", "# Layout: 2 rows (CIB, tSZ) × 3 cols (W012, W200, W201)\n", "# ---------------------------------------------------------------------------\n", "sources_cib = [\n", " ('Agora', mt_agora_cib, 'C0', '-'),\n", " ('DDPM', mt_ddpm_cib, 'C1', '-'),\n", " ('Gaussian', mt_gauss_cib, 'C2', '--'),\n", "]\n", "sources_tsz = [\n", " ('Agora', mt_agora_tsz, 'C0', '-'),\n", " ('DDPM', mt_ddpm_tsz, 'C1', '-'),\n", " ('Gaussian', mt_gauss_tsz, 'C2', '--'),\n", "]\n", "\n", "tensor_titles = {\n", " 'W012': r'$W^{0,2}_1$ — boundary normals',\n", " 'W200': r'$W^{2,0}_0$ — area inertia',\n", " 'W201': r'$W^{2,0}_1$ — boundary positions',\n", "}\n", "\n", "fig, axes = plt.subplots(2, 3, figsize=(15, 8), sharex=True, sharey=True)\n", "\n", "for row, (channel, sources) in enumerate([('CIB', sources_cib), ('tSZ', sources_tsz)]):\n", " for col, tt in enumerate(TENSOR_TYPES):\n", " ax = axes[row, col]\n", " for label, mt, color, ls in sources:\n", " beta = mt[tt]['beta'] # (N, T)\n", " mean = beta.mean(axis=0)\n", " std = beta.std(axis=0)\n", " ax.plot(THRESHOLDS, mean, color=color, ls=ls, lw=1.5, label=label)\n", " ax.fill_between(THRESHOLDS, mean - std, mean + std,\n", " alpha=0.2, color=color)\n", " ax.axhline(1, color='k', lw=0.6, ls=':', alpha=0.5)\n", " ax.set_ylim(0, 1.05)\n", " ax.set_xlim(THRESHOLDS[0], THRESHOLDS[-1])\n", " ax.grid(True, ls=':', lw=0.4, alpha=0.5)\n", " if row == 0:\n", " ax.set_title(tensor_titles[tt], fontsize=12)\n", " if col == 0:\n", " ax.set_ylabel(fr'{channel} — $\\beta = \\lambda_{{\\min}}/\\lambda_{{\\max}}$',\n", " fontsize=11)\n", " if row == 1:\n", " ax.set_xlabel(r'Threshold $\\nu$', fontsize=11)\n", " if row == 0 and col == 2:\n", " ax.legend(fontsize=10)\n", "\n", "plt.suptitle('Minkowski tensor anisotropy index β(ν)\\n'\n", " r'(β = 1 = isotropic; β → 0 = maximally elongated)',\n", " fontsize=13)\n", "plt.tight_layout()\n", "Path('plots').mkdir(exist_ok=True)\n", "plt.savefig('plots/minkowski_tensors_beta.pdf', dpi=200, bbox_inches='tight')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "8d6cb8e4", "metadata": {}, "source": [ "## 6 Orientation distribution θ at representative thresholds\n", "\n", "Polar histograms of the major-axis orientation θ at three representative\n", "excursion levels (low ν = 0.2, mid ν = 0.5, high ν = 0.8). A uniform\n", "distribution around the circle indicates statistical isotropy. Any preferred\n", "direction (e.g. a peak near 0° or 90°) would indicate a systematic alignment\n", "of cluster features — which should not appear in either the AGORA simulations\n", "or a well-trained DDPM unless the simulation has an inherent preferred axis\n", "(e.g. from the scanning strategy or the flat-sky projection)." ] }, { "cell_type": "code", "execution_count": null, "id": "a49be390", "metadata": {}, "outputs": [], "source": [ "# ---------------------------------------------------------------------------\n", "# Plot θ distributions: polar histograms at three representative thresholds\n", "# A statistically isotropic field should give a uniform distribution.\n", "# Deviations indicate preferred orientations in the morphology.\n", "# ---------------------------------------------------------------------------\n", "PROBE_THRESHOLDS = [0.2, 0.5, 0.8] # low / mid / high excursion level\n", "probe_idx = [int(np.argmin(np.abs(THRESHOLDS - t))) for t in PROBE_THRESHOLDS]\n", "\n", "fig, axes = plt.subplots(\n", " 2, len(PROBE_THRESHOLDS),\n", " figsize=(12, 7),\n", " subplot_kw={'projection': 'polar'},\n", ")\n", "\n", "theta_bins = np.linspace(-np.pi / 2, np.pi / 2, 19) # 18 bins × 10°\n", "bin_centres = 0.5 * (theta_bins[:-1] + theta_bins[1:])\n", "bin_width = theta_bins[1] - theta_bins[0]\n", "\n", "for row, (channel, mt_agora, mt_ddpm, mt_gauss) in enumerate([\n", " ('CIB', mt_agora_cib, mt_ddpm_cib, mt_gauss_cib),\n", " ('tSZ', mt_agora_tsz, mt_ddpm_tsz, mt_gauss_tsz),\n", "]):\n", " for col, (t_val, t_idx) in enumerate(zip(PROBE_THRESHOLDS, probe_idx)):\n", " ax = axes[row, col]\n", " for label, mt, color in [\n", " ('Agora', mt_agora, 'C0'),\n", " ('DDPM', mt_ddpm, 'C1'),\n", " ('Gaussian', mt_gauss, 'C2'),\n", " ]:\n", " theta_vals = mt['W012']['theta'][:, t_idx] # (N,)\n", " counts, _ = np.histogram(theta_vals, bins=theta_bins)\n", " counts = counts / counts.sum() # normalise to density\n", " ax.bar(bin_centres, counts, width=bin_width * 0.85,\n", " alpha=0.5, color=color, label=label)\n", " # Uniform reference\n", " ax.axhline(1 / len(bin_centres), color='k', lw=0.8, ls='--', alpha=0.6)\n", " ax.set_theta_zero_location('N')\n", " ax.set_theta_direction(-1)\n", " ax.set_thetalim(-np.pi / 2, np.pi / 2)\n", " ax.set_title(f'{channel} ν = {t_val:.1f}', fontsize=10, pad=8)\n", " if row == 0 and col == 0:\n", " ax.legend(fontsize=9, loc='upper right',\n", " bbox_to_anchor=(1.5, 1.15))\n", "\n", "plt.suptitle(r'Orientation θ of major axis (W012, $\\nu$ = threshold)'\n", " '\\nUniform = isotropic; peaked = preferred orientation',\n", " fontsize=12)\n", "plt.tight_layout()\n", "plt.savefig('plots/minkowski_tensors_theta.pdf', dpi=200, bbox_inches='tight')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "cd8a7ce1", "metadata": {}, "source": [ "## 7 Standardised residuals: (β_Agora − β_DDPM) / σ_Agora\n", "\n", "Residuals for the W012 tensor (most physically interpretable) show whether\n", "any systematic anisotropy bias persists across the threshold range, in units\n", "of the Agora sample variance. Large residuals at high thresholds (rare\n", "excursion sets) may indicate that the DDPM does not accurately reproduce\n", "the statistics of the most extreme cluster configurations." ] }, { "cell_type": "code", "execution_count": null, "id": "1c0690c5", "metadata": {}, "outputs": [], "source": [ "# ---------------------------------------------------------------------------\n", "# Residuals: (β_Agora − β_DDPM) / σ_Agora for W012\n", "# Shows where the DDPM anisotropy deviates from Agora in units of Agora σ.\n", "# ---------------------------------------------------------------------------\n", "fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)\n", "\n", "for ax, channel, mt_agora, mt_ddpm in [\n", " (axes[0], 'CIB', mt_agora_cib, mt_ddpm_cib),\n", " (axes[1], 'tSZ', mt_agora_tsz, mt_ddpm_tsz),\n", "]:\n", " beta_agora = mt_agora['W012']['beta'] # (N, T)\n", " beta_ddpm = mt_ddpm['W012']['beta']\n", " sigma = beta_agora.std(axis=0) + 1e-10\n", " resid = (beta_agora.mean(axis=0) - beta_ddpm.mean(axis=0)) / sigma\n", "\n", " ax.plot(THRESHOLDS, resid, color='C0', lw=1.5)\n", " ax.fill_between(THRESHOLDS, resid, alpha=0.3, color='C0')\n", " ax.axhline(0, color='k', lw=0.8, ls='--')\n", " ax.axhline(1, color='gray', lw=0.6, ls=':')\n", " ax.axhline(-1, color='gray', lw=0.6, ls=':')\n", " ax.set_xlabel(r'Threshold $\\nu$', fontsize=11)\n", " ax.set_title(fr'{channel} — W012 $\\beta$ residual', fontsize=12)\n", " ax.grid(True, ls=':', lw=0.4, alpha=0.5)\n", "\n", "axes[0].set_ylabel(\n", " r'$(\\bar{\\beta}^\\mathrm{Agora} - \\bar{\\beta}^\\mathrm{DDPM})\\,/\\,\\sigma^\\mathrm{Agora}$',\n", " fontsize=11)\n", "\n", "plt.suptitle('W012 anisotropy residuals (Agora − DDPM)', fontsize=13)\n", "plt.tight_layout()\n", "plt.savefig('plots/minkowski_tensors_residuals.pdf', dpi=200, bbox_inches='tight')\n", "plt.show()\n", "\n", "for channel, mt_agora, mt_ddpm in [\n", " ('CIB', mt_agora_cib, mt_ddpm_cib),\n", " ('tSZ', mt_agora_tsz, mt_ddpm_tsz),\n", "]:\n", " sigma = mt_agora['W012']['beta'].std(axis=0) + 1e-10\n", " resid = (mt_agora['W012']['beta'].mean(0) - mt_ddpm['W012']['beta'].mean(0)) / sigma\n", " print(f'{channel}: {(np.abs(resid) < 1).mean():.0%} of thresholds within 1σ')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.11.13" } }, "nbformat": 4, "nbformat_minor": 5 }