{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Bandpass Sampling\n", "\n", "This notebook demonstrates the bandpass sampling capability in PySM, which enables generating realistic detector-specific bandpass variations for large-scale simulation studies.\n", "\n", "## What This Is Used For\n", "\n", "Bandpass sampling is used when you want many plausible detector bandpasses that are consistent with a measured (or nominal) bandpass. Typical uses include:\n", "\n", "1. Propagating bandpass uncertainty into simulated sky maps and time-ordered data.\n", "2. Monte Carlo studies of systematic effects that depend on bandpass shape (for example, foreground mixing and calibration).\n", "3. Building realistic per-detector or per-wafer bandpasses for end-to-end pipeline validation.\n", "\n", "\n", "## Introduction\n", "\n", "Modern Cosmic Microwave Background (CMB) experiments like Simons Observatory have arrays with hundreds of detectors. Each detector has a unique bandpass due to manufacturing variations. For accurate systematic uncertainty studies, we need to:\n", "\n", "1. Generate detector-specific bandpasses that preserve statistical properties\n", "2. Ensure bandpasses are smooth and physically realistic\n", "3. Control the level of variation between detectors\n", "\n", "PySM implements the bandpass sampling methodology used in Map-Based Simulations 16 (MBS 16) for Simons Observatory (https://github.com/simonsobs/map_based_simulations/tree/main/mbs-s0016-20241111#readme), based on kernel density estimation (KDE) with bootstrap resampling.\n", "\n", "## Methodology\n", "\n", "The bandpass sampling process follows these steps:\n", "\n", "### 1. Normalization\n", "Normalize the input bandpass so that $\\int b(\\nu) d\\nu = 1$. This treats the bandpass as a probability density function (PDF) over frequencies.\n", "\n", "### 2. Cumulative Distribution Function (CDF)\n", "Compute the cumulative distribution function (CDF) by integrating the normalized bandpass:\n", "$$P(\\nu) = \\int_{\\nu_{min}}^{\\nu} b(\\nu') d\\nu'$$\n", "\n", "### 3. Bootstrap Resampling\n", "Draw N samples from a uniform distribution U(0,1) and map them to frequencies using the inverse CDF. This creates a discrete frequency sample that statistically represents the bandpass.\n", "\n", "### 4. Kernel Density Estimation\n", "Apply Gaussian KDE to the frequency samples to create a smooth, continuous bandpass:\n", "$$\\hat{b}(\\nu) = \\frac{1}{Nh} \\sum_{i=1}^{N} K\\left(\\frac{\\nu - \\nu_i}{h}\\right)$$\n", "\n", "where K is a Gaussian kernel and h is the bandwidth optimized via leave-one-out cross-validation.\n", "\n", "### 5. Moment Calculation\n", "For each resampled bandpass, compute:\n", "- **Centroid (first moment)**: $\\nu_c = \\int \\nu b(\\nu) d\\nu$\n", "- **Bandwidth (second moment)**: $\\sigma = \\sqrt{\\int (\\nu - \\nu_c)^2 b(\\nu) d\\nu}$\n", "\n", "## Example Usage\n", "\n", "Let's demonstrate with a synthetic bandpass mimicking a realistic CMB detector.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2026-02-05T21:43:04.265553Z", "iopub.status.busy": "2026-02-05T21:43:04.265065Z", "iopub.status.idle": "2026-02-05T21:43:09.276826Z", "shell.execute_reply": "2026-02-05T21:43:09.276221Z" }, "tags": [ "nbval-ignore-output" ] }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import pysm3\n", "\n", "%matplotlib inline\n", "import warnings\n", "warnings.filterwarnings('ignore')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create a Synthetic Bandpass\n", "\n", "We create a realistic bandpass with:\n", "- **Primary response**: Gaussian centered at 100 GHz with 8 GHz width\n", "- **Secondary lobe**: Smaller Gaussian at 110 GHz (common in real detectors)\n", "- **Noise**: Small random fluctuations to simulate measurement uncertainty\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2026-02-05T21:43:09.280869Z", "iopub.status.busy": "2026-02-05T21:43:09.280446Z", "iopub.status.idle": "2026-02-05T21:43:09.478467Z", "shell.execute_reply": "2026-02-05T21:43:09.477758Z" }, "tags": [ "nbval-ignore-output" ] }, "outputs": [], "source": [ "# Create frequency grid from 80 to 120 GHz\n", "nu = np.linspace(80, 120, 200)\n", "\n", "# Primary Gaussian response at 100 GHz\n", "primary = np.exp(-0.5 * ((nu - 100) / 8)**2)\n", "\n", "# Secondary lobe at 110 GHz (10% amplitude)\n", "secondary = 0.1 * np.exp(-0.5 * ((nu - 110) / 3)**2)\n", "\n", "# Add realistic noise\n", "np.random.seed(42)\n", "noise = 0.02 * np.random.randn(len(nu))\n", "bnu = np.maximum(primary + secondary + noise, 0) # Ensure non-negative\n", "\n", "# Visualize the input bandpass\n", "plt.figure(figsize=(10, 5))\n", "plt.plot(nu, bnu, 'k-', linewidth=2, label='Input Bandpass')\n", "plt.xlabel('Frequency [GHz]', fontsize=12)\n", "plt.ylabel('Transmission', fontsize=12)\n", "plt.title('Synthetic Input Bandpass', fontsize=14)\n", "plt.legend(fontsize=11)\n", "plt.grid(True, alpha=0.3)\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute Bandpass Moments\n", "\n", "First, normalize the bandpass and compute its statistical properties:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2026-02-05T21:43:09.513863Z", "iopub.status.busy": "2026-02-05T21:43:09.513702Z", "iopub.status.idle": "2026-02-05T21:43:09.517423Z", "shell.execute_reply": "2026-02-05T21:43:09.516732Z" }, "tags": [ "nbval-ignore-output" ] }, "outputs": [], "source": [ "try:\n", " from numpy import trapezoid\n", "except ImportError:\n", " from numpy import trapz as trapezoid\n", "\n", "# Normalize bandpass\n", "bnu_norm = bnu / trapezoid(bnu, nu)\n", "\n", "# Compute moments\n", "orig_centroid, orig_bandwidth = pysm3.compute_moments(nu, bnu_norm)\n", "\n", "print(f\"Original Bandpass Properties:\")\n", "print(f\" Centroid (\u03bd_c): {orig_centroid:.4f} GHz\")\n", "print(f\" Bandwidth (\u03c3): {orig_bandwidth:.4f} GHz\")\n", "print(f\" Normalization: {trapezoid(bnu_norm, nu):.6f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Resample for Multiple Detectors\n", "\n", "Now generate unique bandpasses for 6 different detectors (representing 6 wafers in a detector array).\n", "\n", "### Parameters:\n", "- **num_wafers**: Number of detector bandpasses to generate\n", "- **bootstrap_size**: Number of frequency samples drawn (controls smoothness vs variation)\n", "- **random_seed**: For reproducibility\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2026-02-05T21:43:09.519356Z", "iopub.status.busy": "2026-02-05T21:43:09.519214Z", "iopub.status.idle": "2026-02-05T21:43:09.528560Z", "shell.execute_reply": "2026-02-05T21:43:09.527886Z" }, "tags": [ "nbval-ignore-output" ] }, "outputs": [], "source": [ "num_wafers = 6\n", "results = pysm3.resample_bandpass(\n", " nu, bnu,\n", " num_wafers=num_wafers,\n", " bootstrap_size=128,\n", " random_seed=1929\n", ")\n", "\n", "print(f\"Generated {num_wafers} unique detector bandpasses\\n\")\n", "print(f\"{'Detector':<10} {'Centroid [GHz]':<18} {'Bandwidth [GHz]':<18} {'\u0394\u03bd_c [GHz]':<15}\")\n", "print(\"-\" * 65)\n", "for i, r in enumerate(results):\n", " delta_cent = r['centroid'] - orig_centroid\n", " print(f\"Wafer {i:<4} {r['centroid']:>16.4f} {r['bandwidth']:>16.4f} {delta_cent:>13.4f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualize Resampled Bandpasses\n", "\n", "Compare each resampled bandpass against the original. Notice how:\n", "- The overall shape is preserved\n", "- Peak locations vary slightly (realistic detector variation)\n", "- Secondary features are maintained\n", "- All curves are smooth and artifact-free\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2026-02-05T21:43:09.530338Z", "iopub.status.busy": "2026-02-05T21:43:09.530171Z", "iopub.status.idle": "2026-02-05T21:43:10.095415Z", "shell.execute_reply": "2026-02-05T21:43:10.094683Z" }, "tags": [ "nbval-ignore-output" ] }, "outputs": [], "source": [ "fig, axes = plt.subplots(2, 3, figsize=(15, 10))\n", "axes = axes.flatten()\n", "\n", "for i, (ax, result) in enumerate(zip(axes, results)):\n", " # Plot original\n", " ax.plot(nu, bnu_norm, 'k--', linewidth=1.5, alpha=0.5, label='Original')\n", " \n", " # Plot resampled\n", " ax.plot(result['frequency'], result['weights'], \n", " 'C1-', linewidth=2, alpha=0.8, label=f\"Wafer {i}\")\n", " \n", " # Formatting\n", " ax.set_xlabel('Frequency [GHz]', fontsize=11)\n", " ax.set_ylabel('Normalized Transmission', fontsize=11)\n", " ax.set_title(f\"Wafer {i}: \u03bd_c = {result['centroid']:.2f} GHz, \"\n", " f\"\u03c3 = {result['bandwidth']:.2f} GHz\", fontsize=11)\n", " ax.legend(fontsize=10)\n", " ax.grid(True, alpha=0.3)\n", "\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Statistical Validation\n", "\n", "To verify the resampling produces statistically valid results, generate many bandpasses and examine the distribution of moments.\n", "\n", "We expect:\n", "- **Centroid distribution**: Centered near the original value with realistic spread\n", "- **Bandwidth distribution**: Similar mean to original, with controlled variation\n", "- **No outliers**: All resampled bandpasses should be physically reasonable\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2026-02-05T21:43:10.098188Z", "iopub.status.busy": "2026-02-05T21:43:10.097956Z", "iopub.status.idle": "2026-02-05T21:43:10.412712Z", "shell.execute_reply": "2026-02-05T21:43:10.411420Z" }, "tags": [ "nbval-ignore-output" ] }, "outputs": [], "source": [ "n_samples = 50\n", "results_many = pysm3.resample_bandpass(\n", " nu, bnu, \n", " num_wafers=n_samples, \n", " bootstrap_size=128, \n", " random_seed=1929\n", ")\n", "\n", "centroids = np.array([r['centroid'] for r in results_many])\n", "bandwidths = np.array([r['bandwidth'] for r in results_many])\n", "\n", "# Create statistical plots\n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))\n", "\n", "# Centroid distribution\n", "ax1.hist(centroids, bins=20, alpha=0.7, color='C0', edgecolor='black', density=False)\n", "ax1.axvline(orig_centroid, color='red', linestyle='--', linewidth=2, \n", " label=f'Original: {orig_centroid:.2f} GHz')\n", "ax1.axvline(np.mean(centroids), color='green', linestyle='--', linewidth=2,\n", " label=f'Mean: {np.mean(centroids):.2f} GHz')\n", "ax1.set_xlabel('Centroid [GHz]', fontsize=12)\n", "ax1.set_ylabel('Count', fontsize=12)\n", "ax1.set_title(f'Centroid Distribution (N={n_samples})', fontsize=13)\n", "ax1.legend(fontsize=11)\n", "ax1.grid(True, alpha=0.3)\n", "\n", "# Bandwidth distribution\n", "ax2.hist(bandwidths, bins=20, alpha=0.7, color='C1', edgecolor='black', density=False)\n", "ax2.axvline(orig_bandwidth, color='red', linestyle='--', linewidth=2,\n", " label=f'Original: {orig_bandwidth:.2f} GHz')\n", "ax2.axvline(np.mean(bandwidths), color='green', linestyle='--', linewidth=2,\n", " label=f'Mean: {np.mean(bandwidths):.2f} GHz')\n", "ax2.set_xlabel('Bandwidth [GHz]', fontsize=12)\n", "ax2.set_ylabel('Count', fontsize=12)\n", "ax2.set_title(f'Bandwidth Distribution (N={n_samples})', fontsize=13)\n", "ax2.legend(fontsize=11)\n", "ax2.grid(True, alpha=0.3)\n", "\n", "plt.tight_layout()\n", "\n", "# Print statistical summary\n", "print(\"\\nStatistical Summary:\")\n", "print(f\"\\nCentroids:\")\n", "print(f\" Mean: {np.mean(centroids):.4f} GHz\")\n", "print(f\" Std Dev: {np.std(centroids):.4f} GHz\")\n", "print(f\" Bias: {np.mean(centroids) - orig_centroid:+.4f} GHz\")\n", "print(f\"\\nBandwidths:\")\n", "print(f\" Mean: {np.mean(bandwidths):.4f} GHz\")\n", "print(f\" Std Dev: {np.std(bandwidths):.4f} GHz\")\n", "print(f\" Bias: {np.mean(bandwidths) - orig_bandwidth:+.4f} GHz\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Overlay Visualization\n", "\n", "Show 10 resampled bandpasses overlaid on the original to visualize the ensemble of variations:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "execution": { "iopub.execute_input": "2026-02-05T21:43:10.414661Z", "iopub.status.busy": "2026-02-05T21:43:10.414487Z", "iopub.status.idle": "2026-02-05T21:43:10.596595Z", "shell.execute_reply": "2026-02-05T21:43:10.595745Z" }, "tags": [ "nbval-ignore-output" ] }, "outputs": [], "source": [ "# Generate 10 samples for overlay\n", "results_overlay = pysm3.resample_bandpass(\n", " nu, bnu, \n", " num_wafers=10, \n", " bootstrap_size=128, \n", " random_seed=2024\n", ")\n", "\n", "plt.figure(figsize=(12, 6))\n", "\n", "# Plot original\n", "plt.plot(nu, bnu_norm, 'k-', linewidth=3, label='Original', zorder=10)\n", "\n", "# Plot resampled with transparency\n", "for i, r in enumerate(results_overlay):\n", " plt.plot(r['frequency'], r['weights'], \n", " 'C1-', linewidth=1.5, alpha=0.3)\n", "\n", "# Add single legend entry for resampled\n", "plt.plot([], [], 'C1-', linewidth=1.5, alpha=0.5, label='Resampled (N=10)')\n", "\n", "plt.xlabel('Frequency [GHz]', fontsize=13)\n", "plt.ylabel('Normalized Transmission', fontsize=13)\n", "plt.title('Bandpass Variation Envelope', fontsize=14)\n", "plt.legend(fontsize=12)\n", "plt.grid(True, alpha=0.3)\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Command Line Usage\n", "\n", "PySM also provides a command-line tool for batch processing:\n", "\n", "```bash\n", "# Resample a bandpass file for 6 detectors\n", "pysm_bandpass_sampler input_bandpass.txt \\\n", " --num-wafers 6 \\\n", " --bootstrap-size 128 \\\n", " --seed 42 \\\n", " --output-dir ./resampled_bandpasses/\n", "```\n", "\n", "Input file format (ASCII, 2 columns):\n", "```\n", "# frequency transmission\n", "80.0 0.01\n", "81.0 0.05\n", "...\n", "```\n", "\n", "Output: IPAC table format files, one per wafer.\n", "\n", "## Applications\n", "\n", "This bandpass sampling capability enables:\n", "\n", "1. **Large-scale simulations**: Generate unique bandpasses for hundreds of detectors\n", "2. **Systematic studies**: Propagate bandpass uncertainties through analysis pipelines \n", "3. **Monte Carlo ensembles**: Create realistic detector variation ensembles\n", "4. **Instrument modeling**: Model detector arrays (e.g., Simons Observatory's 6 wafers per frequency band)\n", "\n", "## Further Reading\n", "\n", "- See the [bandpass sampling comparison notebook](bandpass_sampling_comparison.ipynb) for validation against reference Simons Observatory data\n", "- Thorne et al. 2017, \"The Python Sky Model\" ([arXiv:1608.02841](https://arxiv.org/abs/1608.02841))\n", "- Zonca et al. 2021, \"The Python Sky Model 3\" ([arXiv:2108.01444](https://arxiv.org/abs/2108.01444))\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.2" }, "nbsphinx": { "execute": "always" } }, "nbformat": 4, "nbformat_minor": 4 }