diff --git a/llms-full.txt b/llms-full.txt index 02342769..af19965d 100644 --- a/llms-full.txt +++ b/llms-full.txt @@ -121,6 +121,10 @@ AUTO-GENERATED by PyAutoHands — do not edit by hand; regenerate with generate. - Contents: Prerequisites, Mask, Dataset, Galaxy MGE Basis, Mapping Matrix, Transformed Mapping Matrix ($f$), Data Vector (D), Curvature Matrix (F), Reconstruction (Positive-Negative), Reconstruction (Positive Only), Visibilities Reconstruction, Likelihood Function, Chi Squared, Noise Normalization Term, Calculate The Log Likelihood, Fit, Galaxy Modeling, Wrap Up - [Modeling Features: Multi Gaussian Expansion (Interferometer)](scripts/interferometer/features/multi_gaussian_expansion/modeling.py): A multi-Gaussian expansion (MGE) decomposes a galaxy's light into ~15-100 Gaussians, where the `intensity` of every Gaussian is solved for via linear algebra using a process called an "inversion" (see the `linear_light_profiles` feature for a full description of this). - Contents: Advantages & Disadvantages, NUFFT (nufftax), Positive Only Solver, Model, Mask, Dataset, Over Sampling, Search, Analysis, VRAM, Run Time, Result, Wrap Up +- [Modeling Features: Operated Light Profiles (Interferometer)](scripts/interferometer/features/operated_light_profile/modeling.py): It is common for galaxies to have point-source emission, for example bright emission right at their centre due to an active galactic nuclei or a very compact knot of star formation. + - Contents: Mask, Dataset, Dataset Auto-Simulation, Over Sampling, Model, Search, Analysis, VRAM, Run Time, Model-Fit, Result +- [Simulator: Operated Light Profiles (Interferometer)](scripts/interferometer/features/operated_light_profile/simulator.py): This script simulates `Interferometer` data of a galaxy using light profiles where: + - Contents: Dataset Paths, Grid, uv-wavelengths, Simulator, Galaxies, Output, Visualize, Galaxies json - [Features: Pixelization Fit](scripts/interferometer/features/pixelization/fit.py): This script performs a single, direct fit (no non-linear search) of the `clumpy` interferometer dataset — a galaxy with a smooth central bulge plus asymmetric clumpy star formation. The fit uses a pixelization with a rectangular mesh and constant regularization scheme to reconstruct the galaxy's surface brightness on a pixel grid. - Contents: Advantages, Disadvantages, Positive Only Solver, Mask, Dataset, Dataset Auto-Simulation, Sparse Operators, Settings, Over Sampling, Mesh Shape, Pixelization, Fit, Wrap Up, Linear Objects, Grids, Reconstruction, Mapped Reconstructed Images, Linear Algebra Matrices (Advanced), Evidence Terms (Advanced), Simulated Interferometer, Future Ideas / Contributions - [Pixelization: Galaxy Reconstruction](scripts/interferometer/features/pixelization/galaxy_reconstruction.py): A common pixelization use-case is to reconstruct a galaxy’s surface brightness on a pixelization mesh, and then export this reconstruction to perform scientific analysis. diff --git a/notebooks/interferometer/features/operated_light_profile/README.md b/notebooks/interferometer/features/operated_light_profile/README.md new file mode 100644 index 00000000..b37a3d4e --- /dev/null +++ b/notebooks/interferometer/features/operated_light_profile/README.md @@ -0,0 +1,20 @@ +The `operated_light_profile` folder contains example scripts showing how to perform analysis of +`Interferometer` data using operated light profiles, which represent compact point-source emission (e.g. an +AGN) whose image-plane shape is specified directly. + +For interferometer data there is no PSF, so operated light profiles are Fourier transformed to the visibility +plane like every other light profile — the PSF-bypass behaviour of the imaging examples applies only where a +PSF exists. Using them keeps a galaxy model consistent across imaging and interferometer datasets. + +# Files + +The following example scripts illustrating galaxy modeling where: + +- `modeling`: Galaxy modeling of an `Interferometer` dataset using operated light profiles. +- `simulator`: Simulating interferometer data of a galaxy using operated light profiles. + +# Results + +These scripts only give a brief overview of how to analyse and interpret the results a galaxy model fit. + +A full guide to result analysis is given at `autogalaxy_workspace/*/results`. diff --git a/notebooks/interferometer/features/operated_light_profile/modeling.ipynb b/notebooks/interferometer/features/operated_light_profile/modeling.ipynb new file mode 100644 index 00000000..6f764e87 --- /dev/null +++ b/notebooks/interferometer/features/operated_light_profile/modeling.ipynb @@ -0,0 +1,436 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Modeling Features: Operated Light Profiles (Interferometer)\n", + "===========================================================\n", + "\n", + "It is common for galaxies to have point-source emission, for example bright emission right at their centre due\n", + "to an active galactic nuclei or a very compact knot of star formation.\n", + "\n", + "For CCD imaging data this emission is blurred by the telescope's Point Spread Function, and the\n", + "`imaging/features/operated_light_profile` example explains how operated light profiles fit it by assuming the\n", + "profile has already been convolved with the PSF.\n", + "\n", + "Interferometer data has no PSF: the visibilities are the Fourier transform of the galaxy's sky emission, and\n", + "the synthesized beam only enters when a dirty image is formed. An operated light profile therefore takes on a\n", + "simpler meaning \u2014 it is a light profile whose image-plane shape directly represents the compact emission,\n", + "with no convolution step to bypass. Its image is Fourier transformed to the visibility plane like every other\n", + "light profile.\n", + "\n", + "Using operated light profiles for this compact emission keeps a galaxy model consistent across datasets: the\n", + "same `lp_operated` / `lp_linear_operated` component fitted to imaging data can be fitted to interferometer\n", + "data of the same galaxy, with the PSF-bypass behaviour applying only where a PSF exists.\n", + "\n", + "__Advanced: Visibility-Space Overrides__\n", + "\n", + "Internally, linear operated light profiles use the inversion's `operated_mapping_matrix_override` API to\n", + "bypass PSF convolution for imaging data. Interferometer inversions also support this override for custom\n", + "linear objects: an override supplied to an interferometer inversion bypasses the NUFFT entirely and must\n", + "therefore be a complex matrix in visibility space, of shape [total_visibilities, params] (e.g. computed via\n", + "an analytic Fourier transform). See the `LinearObj.operated_mapping_matrix_override` docstring in PyAutoArray\n", + "for the full contract. The linear operated light profiles fitted in this example do not use an override for\n", + "interferometer data \u2014 their images are NUFFT'd like any other profile.\n", + "\n", + "__Model__\n", + "\n", + "This script fits an `Interferometer` dataset of a galaxy with a model where:\n", + "\n", + " - The galaxy's light is a linear parametric `Sersic` bulge.\n", + " - The galaxy includes a linear parametric operated `Gaussian` representing its compact nuclear emission.\n", + "\n", + "__Fit__\n", + "\n", + "For operated light profiles, there is no `fit.py` example found for standard light profiles, linear light\n", + "profiles and other examples.\n", + "\n", + "This is done purely to keep the number of examples in the workspace manageable. To perform a fit with operated\n", + "light profiles, simply follow one of the other `interferometer/fit.py` examples and replace the light profiles\n", + "with operated light profiles using the API described below.\n", + "\n", + "__Start Here Notebook__\n", + "\n", + "If any code in this script is unclear, refer to the `interferometer/start_here.ipynb` notebook.\n", + "\n", + "__Imaging Equivalent__\n", + "\n", + "For the CCD-imaging version of this script, see\n", + "`autogalaxy_workspace/*/imaging/features/operated_light_profile/modeling.py`.\n", + "\n", + "__Contents__\n", + "\n", + "- **Mask:** Define the `real_space_mask` which sets the grid the galaxy is evaluated on.\n", + "- **Dataset:** Loading the operated light profile interferometer dataset from FITS files.\n", + "- **Dataset Auto-Simulation:** Automatically simulating the dataset if it does not already exist.\n", + "- **Over Sampling:** Interferometer modeling does not use over-sampling (covered briefly here).\n", + "- **Model:** Composing the galaxy model with a linear Sersic bulge and operated Gaussian point source.\n", + "- **Search:** Configuring the Nautilus nested sampling non-linear search.\n", + "- **Analysis:** Creating the AnalysisInterferometer object for likelihood evaluation.\n", + "- **VRAM:** Discussion of GPU VRAM usage for operated light profiles.\n", + "- **Run Time:** Discussion of computational run times for operated light profiles.\n", + "- **Model-Fit:** Running the model-fit and monitoring output.\n", + "- **Result:** Inspecting the result object and best-fit model." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Google Colab Setup__\n", + "\n", + "This cell sets up the environment when the notebook is run on Google Colab: it installs the\n", + "required PyAuto packages, clones the workspace (configuration files and example datasets) and\n", + "points the configuration at it. If you are running the notebook elsewhere (e.g. locally via\n", + "your own installation) it does nothing, and you can run it safely.\n", + "\n", + "Colab tip: model-fits run much faster on a GPU \u2014 enable one via \"Runtime\" -> \"Change runtime\n", + "type\" -> \"Hardware accelerator\" before running the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " import google.colab\n", + "except ImportError:\n", + " from autogalaxy import setup_colab as _setup_colab\n", + "else:\n", + " import importlib\n", + " import subprocess\n", + " import sys\n", + "\n", + " subprocess.check_call(\n", + " [sys.executable, \"-m\", \"pip\", \"install\", \"autonerves\", \"--no-deps\"]\n", + " )\n", + " _setup_colab = importlib.import_module(\"autonerves.setup_colab\")\n", + "\n", + "_setup_colab.setup(\"autogalaxy\")" + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "from autogalaxy import setup_notebook; setup_notebook()\n", + "\n", + "from pathlib import Path\n", + "import autofit as af\n", + "import autogalaxy as ag\n", + "import autogalaxy.plot as aplt" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Mask__\n", + "\n", + "We define the `real_space_mask` which defines the grid the image of the galaxy is evaluated on." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "mask_radius = 3.0\n", + "\n", + "real_space_mask = ag.Mask2D.circular(\n", + " shape_native=(256, 256),\n", + " pixel_scales=0.1,\n", + " radius=mask_radius,\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dataset__\n", + "\n", + "Load and plot the galaxy `Interferometer` dataset `operated` from .fits files, using `TransformerNUFFT`\n", + "backed by `nufftax`." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dataset_name = \"operated\"\n", + "dataset_path = Path(\"dataset\") / \"interferometer\" / dataset_name" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dataset Auto-Simulation__\n", + "\n", + "If the dataset does not already exist on your system, it will be created by running the corresponding\n", + "simulator script. This ensures that all example scripts can be run without manually simulating data first." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "if ag.util.dataset.should_simulate(str(dataset_path)):\n", + " import subprocess\n", + " import sys\n", + "\n", + " subprocess.run(\n", + " [\n", + " sys.executable,\n", + " \"scripts/interferometer/features/operated_light_profile/simulator.py\",\n", + " ],\n", + " check=True,\n", + " )\n", + "\n", + "dataset = ag.Interferometer.from_fits(\n", + " data_path=dataset_path / \"data.fits\",\n", + " noise_map_path=dataset_path / \"noise_map.fits\",\n", + " uv_wavelengths_path=dataset_path / \"uv_wavelengths.fits\",\n", + " real_space_mask=real_space_mask,\n", + " transformer_class=ag.TransformerNUFFT,\n", + ")\n", + "\n", + "aplt.subplot_interferometer_dirty_images(dataset=dataset)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Over Sampling__\n", + "\n", + "If you are familiar with using imaging data, you may have seen that a numerical technique called over sampling\n", + "is used, which evaluates light profiles on a higher resolution grid than the image data to ensure the\n", + "calculation is accurate.\n", + "\n", + "Interferometer data does not observe galaxies in a way where over sampling is necessary, therefore all\n", + "interferometer calculations are performed without over sampling.\n", + "\n", + "__Model__\n", + "\n", + "We compose our model where in this example:\n", + "\n", + " - The galaxy's bulge is a linear parametric `Sersic` bulge [6 parameters].\n", + " - The galaxy's point source emission is a linear parametric operated `Gaussian` centred on the bulge\n", + " [3 parameters].\n", + "\n", + "The number of free parameters and therefore the dimensionality of non-linear parameter space is N=9.\n", + "\n", + "The prior on the operated `Gaussian`'s `sigma` value is very important, as it is often the case that this is a\n", + "very small value (e.g. ~0.1).\n", + "\n", + "By default, **PyAutoGalaxy** assumes a `UniformPrior` from 0.0 to 5.0, but the scale of this value depends on\n", + "the resolution of the data. I therefore recommend you set it manually below, using your knowledge of the\n", + "compact emission's angular size." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "bulge = af.Model(ag.lp_linear.Sersic)\n", + "psf = af.Model(ag.lp_linear_operated.Gaussian)\n", + "\n", + "psf.sigma = af.UniformPrior(lower_limit=0.0, upper_limit=5.0)\n", + "\n", + "bulge.centre = psf.centre\n", + "\n", + "galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge, psf=psf)\n", + "\n", + "model = af.Collection(galaxies=af.Collection(galaxy=galaxy))" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The `info` attribute shows the model in a readable format." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "print(model.info)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Search__\n", + "\n", + "The model is fitted to the data using a non-linear search. In this example, we use the nested sampling\n", + "algorithm Nautilus (https://nautilus.readthedocs.io/en/latest/).\n", + "\n", + "A full description of the settings below is given in the beginner modeling scripts, if anything is unclear." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "search = af.Nautilus(\n", + " path_prefix=Path(\"interferometer\") / \"features\",\n", + " name=\"operated_light_profiles\",\n", + " unique_tag=dataset_name,\n", + " n_live=100,\n", + " n_batch=20, # GPU galaxy model fits are batched and run simultaneously, see VRAM section below.\n", + " live_visual_update=False, # Set True to open a live matplotlib window (script) or refresh a Jupyter cell (notebook).\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Analysis__\n", + "\n", + "Create the `AnalysisInterferometer` object defining how the via Nautilus the model is fitted to the data." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "analysis = ag.AnalysisInterferometer(dataset=dataset, use_jax=True)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__VRAM__\n", + "\n", + "The `interferometer/modeling.py` example explains how VRAM is used during GPU-based fitting and how to print\n", + "the estimated VRAM required by a model.\n", + "\n", + "For each linear operated light profile in the model a small additional amount of VRAM is used to store its\n", + "NUFFT'd mapping matrix column. For 1-10 linear light profiles this is a tiny amount of VRAM (e.g. < 10MB per\n", + "batched likelihood). Even for large batch sizes you almost certainly will not use enough VRAM to require\n", + "monitoring.\n", + "\n", + "VRAM on interferometer datasets is driven primarily by the visibility count and the real-space mask size, not\n", + "the number of operated light profiles in the model.\n", + "\n", + "__Run Time__\n", + "\n", + "For interferometer data the likelihood evaluation time of an operated light profile is the same as that of an\n", + "ordinary light profile \u2014 both are evaluated in real space and NUFFT'd to the visibility plane (the PSF\n", + "convolution that operated profiles bypass for imaging data does not exist here).\n", + "\n", + "The overall run-time may be a little slower than a bulge-only model though, because the `psf` component adds a\n", + "few extra parameters.\n", + "\n", + "__Model-Fit__\n", + "\n", + "We begin the model-fit by passing the model and analysis object to the non-linear search (checkout the output\n", + "folder for on-the-fly visualization and results)." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "result = search.fit(model=model, analysis=analysis)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Result__\n", + "\n", + "The search returns a result object, which whose `info` attribute shows the result in a readable format:" + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "print(result.info)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We plot the maximum likelihood fit, galaxy images and posteriors inferred via Nautilus.\n", + "\n", + "The galaxy bulge and compact nuclear emission appear similar to those in the data, confirming that the\n", + "`intensity` values inferred by the inversion process are accurate." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "print(result.max_log_likelihood_instance)\n", + "\n", + "aplt.subplot_galaxies(galaxies=result.max_log_likelihood_galaxies, grid=result.grids.lp)\n", + "\n", + "aplt.subplot_fit_interferometer(fit=result.max_log_likelihood_fit)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Checkout `autogalaxy_workspace/*/guides/results/start_here.py` for a full description of the result object." + ] + } + ], + "metadata": { + "anaconda-cloud": {}, + "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.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/notebooks/interferometer/features/operated_light_profile/simulator.ipynb b/notebooks/interferometer/features/operated_light_profile/simulator.ipynb new file mode 100644 index 00000000..5282c963 --- /dev/null +++ b/notebooks/interferometer/features/operated_light_profile/simulator.ipynb @@ -0,0 +1,366 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Simulator: Operated Light Profiles (Interferometer)\n", + "===================================================\n", + "\n", + "This script simulates `Interferometer` data of a galaxy using light profiles where:\n", + "\n", + " - The galaxy's bulge is an `Sersic`.\n", + " - The galaxy has a compact component of point-source emission at its centre which is modeled as an operated\n", + " `Gaussian`.\n", + "\n", + "For interferometer data there is no Point Spread Function: the visibilities are the Fourier transform of the\n", + "galaxy's sky emission, and the synthesized beam only enters when a dirty image is formed. The operated\n", + "`Gaussian` therefore represents compact nuclear emission (e.g. an AGN or a compact knot of star formation)\n", + "whose image-plane shape is specified directly, and it is Fourier transformed to the visibility plane like\n", + "every other light profile.\n", + "\n", + "The output dataset is consumed by the companion `modeling.py` in the same folder.\n", + "\n", + "__Advanced__\n", + "\n", + "This is an advanced simulator script, meaning that detailed explanations of certain code are omitted. Refer to\n", + "simulators not in the `advanced` folder for more detailed comments.\n", + "\n", + "__Start Here Notebook__\n", + "\n", + "If any code in this script is unclear, refer to the `interferometer/simulator.ipynb` notebook.\n", + "\n", + "__Contents__\n", + "\n", + "- **Dataset Paths:** Defining the output path for the simulated dataset.\n", + "- **Grid:** Real-space grid the galaxy images are evaluated on.\n", + "- **uv-wavelengths:** Load the uv baselines used to NUFFT the image to the visibility plane.\n", + "- **Simulator:** `SimulatorInterferometer` (no PSF; uv-plane noise instead of image-plane Poisson noise).\n", + "- **Galaxies:** Defining the galaxy with a Sersic bulge and an operated Gaussian point source.\n", + "- **Output:** Saving the simulated dataset to FITS files.\n", + "- **Visualize:** Outputting subplot and dirty-image visualizations as PNG files.\n", + "- **Galaxies json:** Saving the Galaxies object as a JSON file for future reference." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Google Colab Setup__\n", + "\n", + "This cell sets up the environment when the notebook is run on Google Colab: it installs the\n", + "required PyAuto packages, clones the workspace (configuration files and example datasets) and\n", + "points the configuration at it. If you are running the notebook elsewhere (e.g. locally via\n", + "your own installation) it does nothing, and you can run it safely.\n", + "\n", + "Colab tip: model-fits run much faster on a GPU \u2014 enable one via \"Runtime\" -> \"Change runtime\n", + "type\" -> \"Hardware accelerator\" before running the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " import google.colab\n", + "except ImportError:\n", + " from autogalaxy import setup_colab as _setup_colab\n", + "else:\n", + " import importlib\n", + " import subprocess\n", + " import sys\n", + "\n", + " subprocess.check_call(\n", + " [sys.executable, \"-m\", \"pip\", \"install\", \"autonerves\", \"--no-deps\"]\n", + " )\n", + " _setup_colab = importlib.import_module(\"autonerves.setup_colab\")\n", + "\n", + "_setup_colab.setup(\"autogalaxy\")" + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "from autogalaxy import setup_notebook; setup_notebook()\n", + "\n", + "from pathlib import Path\n", + "import autogalaxy as ag\n", + "import autogalaxy.plot as aplt" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dataset Paths__\n", + "\n", + "The `dataset_type` describes the type of data being simulated and `dataset_name` gives it a descriptive name." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dataset_type = \"interferometer\"\n", + "dataset_name = \"operated\"\n", + "dataset_path = Path(\"dataset\", dataset_type, dataset_name)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Grid__\n", + "\n", + "Simulate the image using a (y,x) grid. Over-sampling is an imaging-only technique and is not used for\n", + "interferometer data." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "grid = ag.Grid2D.uniform(shape_native=(256, 256), pixel_scales=0.1)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__uv-wavelengths__\n", + "\n", + "To perform the Fourier transform we need the wavelengths of the baselines." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "uv_wavelengths_path = Path(\"dataset\", dataset_type, \"uv_wavelengths\")\n", + "uv_wavelengths = ag.ndarray_via_fits_from(\n", + " file_path=Path(uv_wavelengths_path, \"sma.fits\"), hdu=0\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Simulator__\n", + "\n", + "Create the simulator for the interferometer data, which defines the exposure time, visibility-plane\n", + "noise sigma, and transformer." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "simulator = ag.SimulatorInterferometer(\n", + " uv_wavelengths=uv_wavelengths,\n", + " exposure_time=300.0,\n", + " noise_sigma=1000.0,\n", + " transformer_class=ag.TransformerDFT,\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Galaxies__\n", + "\n", + "Setup the galaxy with a bulge (elliptical Sersic) for this simulation.\n", + "\n", + "This includes an operated `Gaussian` component which represents the compact point-source emission at the\n", + "galaxy's centre. Its image-plane shape is specified directly by the profile and is Fourier transformed to\n", + "the visibility plane like every other light profile (there is no PSF convolution for interferometer data)." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "galaxy = ag.Galaxy(\n", + " redshift=0.5,\n", + " bulge=ag.lp.Sersic(\n", + " centre=(0.0, 0.0),\n", + " ell_comps=ag.convert.ell_comps_from(axis_ratio=0.9, angle=45.0),\n", + " intensity=1.0,\n", + " effective_radius=0.6,\n", + " sersic_index=3.0,\n", + " ),\n", + " psf=ag.lp_operated.Gaussian(\n", + " centre=(0.0, 0.0), ell_comps=(0.0, 0.0), intensity=100.0, sigma=0.1\n", + " ),\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Use these galaxies to generate the image which is simulated as an `Interferometer` dataset." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "galaxies = ag.Galaxies(galaxies=[galaxy])\n", + "\n", + "aplt.plot_array(array=galaxies.image_2d_from(grid=grid), title=\"Image\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Pass the simulator galaxies, which creates the real-space image and NUFFTs it to visibilities." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dataset = simulator.via_galaxies_from(galaxies=galaxies, grid=grid)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the simulated `Interferometer` dataset before outputting it to fits." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "aplt.subplot_interferometer_dirty_images(dataset=dataset)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Output__\n", + "\n", + "Output the simulated dataset to the dataset path as .fits files." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "aplt.fits_interferometer(\n", + " dataset=dataset,\n", + " data_path=dataset_path / \"data.fits\",\n", + " noise_map_path=dataset_path / \"noise_map.fits\",\n", + " uv_wavelengths_path=dataset_path / \"uv_wavelengths.fits\",\n", + " overwrite=True,\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Visualize__\n", + "\n", + "Output a subplot of the simulated dataset and the galaxies' images to the dataset path as .png files." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "aplt.subplot_interferometer_dirty_images(\n", + " dataset=dataset, output_path=dataset_path, output_format=\"png\"\n", + ")\n", + "aplt.subplot_galaxies(\n", + " galaxies=galaxies, grid=grid, output_path=dataset_path, output_format=\"png\"\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Galaxies json__\n", + "\n", + "Save the `Galaxies` in the dataset folder as a .json file, ensuring the true light profiles and galaxies\n", + "are safely stored and available to check how the dataset was simulated in the future.\n", + "\n", + "This can be loaded via the method `galaxies = ag.from_json()`." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "ag.output_to_json(\n", + " obj=galaxies,\n", + " file_path=Path(dataset_path, \"galaxies.json\"),\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The dataset can be viewed in the folder `autogalaxy_workspace/dataset/interferometer/operated`." + ] + } + ], + "metadata": { + "anaconda-cloud": {}, + "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.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/scripts/interferometer/features/operated_light_profile/README.md b/scripts/interferometer/features/operated_light_profile/README.md new file mode 100644 index 00000000..b37a3d4e --- /dev/null +++ b/scripts/interferometer/features/operated_light_profile/README.md @@ -0,0 +1,20 @@ +The `operated_light_profile` folder contains example scripts showing how to perform analysis of +`Interferometer` data using operated light profiles, which represent compact point-source emission (e.g. an +AGN) whose image-plane shape is specified directly. + +For interferometer data there is no PSF, so operated light profiles are Fourier transformed to the visibility +plane like every other light profile — the PSF-bypass behaviour of the imaging examples applies only where a +PSF exists. Using them keeps a galaxy model consistent across imaging and interferometer datasets. + +# Files + +The following example scripts illustrating galaxy modeling where: + +- `modeling`: Galaxy modeling of an `Interferometer` dataset using operated light profiles. +- `simulator`: Simulating interferometer data of a galaxy using operated light profiles. + +# Results + +These scripts only give a brief overview of how to analyse and interpret the results a galaxy model fit. + +A full guide to result analysis is given at `autogalaxy_workspace/*/results`. diff --git a/scripts/interferometer/features/operated_light_profile/__init__.py b/scripts/interferometer/features/operated_light_profile/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/scripts/interferometer/features/operated_light_profile/modeling.py b/scripts/interferometer/features/operated_light_profile/modeling.py new file mode 100644 index 00000000..b077e1eb --- /dev/null +++ b/scripts/interferometer/features/operated_light_profile/modeling.py @@ -0,0 +1,247 @@ +""" +Modeling Features: Operated Light Profiles (Interferometer) +=========================================================== + +It is common for galaxies to have point-source emission, for example bright emission right at their centre due +to an active galactic nuclei or a very compact knot of star formation. + +For CCD imaging data this emission is blurred by the telescope's Point Spread Function, and the +`imaging/features/operated_light_profile` example explains how operated light profiles fit it by assuming the +profile has already been convolved with the PSF. + +Interferometer data has no PSF: the visibilities are the Fourier transform of the galaxy's sky emission, and +the synthesized beam only enters when a dirty image is formed. An operated light profile therefore takes on a +simpler meaning — it is a light profile whose image-plane shape directly represents the compact emission, +with no convolution step to bypass. Its image is Fourier transformed to the visibility plane like every other +light profile. + +Using operated light profiles for this compact emission keeps a galaxy model consistent across datasets: the +same `lp_operated` / `lp_linear_operated` component fitted to imaging data can be fitted to interferometer +data of the same galaxy, with the PSF-bypass behaviour applying only where a PSF exists. + +__Advanced: Visibility-Space Overrides__ + +Internally, linear operated light profiles use the inversion's `operated_mapping_matrix_override` API to +bypass PSF convolution for imaging data. Interferometer inversions also support this override for custom +linear objects: an override supplied to an interferometer inversion bypasses the NUFFT entirely and must +therefore be a complex matrix in visibility space, of shape [total_visibilities, params] (e.g. computed via +an analytic Fourier transform). See the `LinearObj.operated_mapping_matrix_override` docstring in PyAutoArray +for the full contract. The linear operated light profiles fitted in this example do not use an override for +interferometer data — their images are NUFFT'd like any other profile. + +__Model__ + +This script fits an `Interferometer` dataset of a galaxy with a model where: + + - The galaxy's light is a linear parametric `Sersic` bulge. + - The galaxy includes a linear parametric operated `Gaussian` representing its compact nuclear emission. + +__Fit__ + +For operated light profiles, there is no `fit.py` example found for standard light profiles, linear light +profiles and other examples. + +This is done purely to keep the number of examples in the workspace manageable. To perform a fit with operated +light profiles, simply follow one of the other `interferometer/fit.py` examples and replace the light profiles +with operated light profiles using the API described below. + +__Start Here Notebook__ + +If any code in this script is unclear, refer to the `interferometer/start_here.ipynb` notebook. + +__Imaging Equivalent__ + +For the CCD-imaging version of this script, see +`autogalaxy_workspace/*/imaging/features/operated_light_profile/modeling.py`. + +__Contents__ + +- **Mask:** Define the `real_space_mask` which sets the grid the galaxy is evaluated on. +- **Dataset:** Loading the operated light profile interferometer dataset from FITS files. +- **Dataset Auto-Simulation:** Automatically simulating the dataset if it does not already exist. +- **Over Sampling:** Interferometer modeling does not use over-sampling (covered briefly here). +- **Model:** Composing the galaxy model with a linear Sersic bulge and operated Gaussian point source. +- **Search:** Configuring the Nautilus nested sampling non-linear search. +- **Analysis:** Creating the AnalysisInterferometer object for likelihood evaluation. +- **VRAM:** Discussion of GPU VRAM usage for operated light profiles. +- **Run Time:** Discussion of computational run times for operated light profiles. +- **Model-Fit:** Running the model-fit and monitoring output. +- **Result:** Inspecting the result object and best-fit model. +""" + +# from autogalaxy import setup_notebook; setup_notebook() + +from pathlib import Path +import autofit as af +import autogalaxy as ag +import autogalaxy.plot as aplt + +""" +__Mask__ + +We define the `real_space_mask` which defines the grid the image of the galaxy is evaluated on. +""" +mask_radius = 3.0 + +real_space_mask = ag.Mask2D.circular( + shape_native=(256, 256), + pixel_scales=0.1, + radius=mask_radius, +) + +""" +__Dataset__ + +Load and plot the galaxy `Interferometer` dataset `operated` from .fits files, using `TransformerNUFFT` +backed by `nufftax`. +""" +dataset_name = "operated" +dataset_path = Path("dataset") / "interferometer" / dataset_name + +""" +__Dataset Auto-Simulation__ + +If the dataset does not already exist on your system, it will be created by running the corresponding +simulator script. This ensures that all example scripts can be run without manually simulating data first. +""" +if ag.util.dataset.should_simulate(str(dataset_path)): + import subprocess + import sys + + subprocess.run( + [ + sys.executable, + "scripts/interferometer/features/operated_light_profile/simulator.py", + ], + check=True, + ) + +dataset = ag.Interferometer.from_fits( + data_path=dataset_path / "data.fits", + noise_map_path=dataset_path / "noise_map.fits", + uv_wavelengths_path=dataset_path / "uv_wavelengths.fits", + real_space_mask=real_space_mask, + transformer_class=ag.TransformerNUFFT, +) + +aplt.subplot_interferometer_dirty_images(dataset=dataset) + +""" +__Over Sampling__ + +If you are familiar with using imaging data, you may have seen that a numerical technique called over sampling +is used, which evaluates light profiles on a higher resolution grid than the image data to ensure the +calculation is accurate. + +Interferometer data does not observe galaxies in a way where over sampling is necessary, therefore all +interferometer calculations are performed without over sampling. + +__Model__ + +We compose our model where in this example: + + - The galaxy's bulge is a linear parametric `Sersic` bulge [6 parameters]. + - The galaxy's point source emission is a linear parametric operated `Gaussian` centred on the bulge + [3 parameters]. + +The number of free parameters and therefore the dimensionality of non-linear parameter space is N=9. + +The prior on the operated `Gaussian`'s `sigma` value is very important, as it is often the case that this is a +very small value (e.g. ~0.1). + +By default, **PyAutoGalaxy** assumes a `UniformPrior` from 0.0 to 5.0, but the scale of this value depends on +the resolution of the data. I therefore recommend you set it manually below, using your knowledge of the +compact emission's angular size. +""" +bulge = af.Model(ag.lp_linear.Sersic) +psf = af.Model(ag.lp_linear_operated.Gaussian) + +psf.sigma = af.UniformPrior(lower_limit=0.0, upper_limit=5.0) + +bulge.centre = psf.centre + +galaxy = af.Model(ag.Galaxy, redshift=0.5, bulge=bulge, psf=psf) + +model = af.Collection(galaxies=af.Collection(galaxy=galaxy)) + +""" +The `info` attribute shows the model in a readable format. +""" +print(model.info) + +""" +__Search__ + +The model is fitted to the data using a non-linear search. In this example, we use the nested sampling +algorithm Nautilus (https://nautilus.readthedocs.io/en/latest/). + +A full description of the settings below is given in the beginner modeling scripts, if anything is unclear. +""" +search = af.Nautilus( + path_prefix=Path("interferometer") / "features", + name="operated_light_profiles", + unique_tag=dataset_name, + n_live=100, + n_batch=20, # GPU galaxy model fits are batched and run simultaneously, see VRAM section below. + live_visual_update=False, # Set True to open a live matplotlib window (script) or refresh a Jupyter cell (notebook). +) + +""" +__Analysis__ + +Create the `AnalysisInterferometer` object defining how the via Nautilus the model is fitted to the data. +""" +analysis = ag.AnalysisInterferometer(dataset=dataset, use_jax=True) + +""" +__VRAM__ + +The `interferometer/modeling.py` example explains how VRAM is used during GPU-based fitting and how to print +the estimated VRAM required by a model. + +For each linear operated light profile in the model a small additional amount of VRAM is used to store its +NUFFT'd mapping matrix column. For 1-10 linear light profiles this is a tiny amount of VRAM (e.g. < 10MB per +batched likelihood). Even for large batch sizes you almost certainly will not use enough VRAM to require +monitoring. + +VRAM on interferometer datasets is driven primarily by the visibility count and the real-space mask size, not +the number of operated light profiles in the model. + +__Run Time__ + +For interferometer data the likelihood evaluation time of an operated light profile is the same as that of an +ordinary light profile — both are evaluated in real space and NUFFT'd to the visibility plane (the PSF +convolution that operated profiles bypass for imaging data does not exist here). + +The overall run-time may be a little slower than a bulge-only model though, because the `psf` component adds a +few extra parameters. + +__Model-Fit__ + +We begin the model-fit by passing the model and analysis object to the non-linear search (checkout the output +folder for on-the-fly visualization and results). +""" +result = search.fit(model=model, analysis=analysis) + +""" +__Result__ + +The search returns a result object, which whose `info` attribute shows the result in a readable format: +""" +print(result.info) + +""" +We plot the maximum likelihood fit, galaxy images and posteriors inferred via Nautilus. + +The galaxy bulge and compact nuclear emission appear similar to those in the data, confirming that the +`intensity` values inferred by the inversion process are accurate. +""" +print(result.max_log_likelihood_instance) + +aplt.subplot_galaxies(galaxies=result.max_log_likelihood_galaxies, grid=result.grids.lp) + +aplt.subplot_fit_interferometer(fit=result.max_log_likelihood_fit) + +""" +Checkout `autogalaxy_workspace/*/guides/results/start_here.py` for a full description of the result object. +""" diff --git a/scripts/interferometer/features/operated_light_profile/simulator.py b/scripts/interferometer/features/operated_light_profile/simulator.py new file mode 100644 index 00000000..03b0fd70 --- /dev/null +++ b/scripts/interferometer/features/operated_light_profile/simulator.py @@ -0,0 +1,166 @@ +""" +Simulator: Operated Light Profiles (Interferometer) +=================================================== + +This script simulates `Interferometer` data of a galaxy using light profiles where: + + - The galaxy's bulge is an `Sersic`. + - The galaxy has a compact component of point-source emission at its centre which is modeled as an operated + `Gaussian`. + +For interferometer data there is no Point Spread Function: the visibilities are the Fourier transform of the +galaxy's sky emission, and the synthesized beam only enters when a dirty image is formed. The operated +`Gaussian` therefore represents compact nuclear emission (e.g. an AGN or a compact knot of star formation) +whose image-plane shape is specified directly, and it is Fourier transformed to the visibility plane like +every other light profile. + +The output dataset is consumed by the companion `modeling.py` in the same folder. + +__Advanced__ + +This is an advanced simulator script, meaning that detailed explanations of certain code are omitted. Refer to +simulators not in the `advanced` folder for more detailed comments. + +__Start Here Notebook__ + +If any code in this script is unclear, refer to the `interferometer/simulator.ipynb` notebook. + +__Contents__ + +- **Dataset Paths:** Defining the output path for the simulated dataset. +- **Grid:** Real-space grid the galaxy images are evaluated on. +- **uv-wavelengths:** Load the uv baselines used to NUFFT the image to the visibility plane. +- **Simulator:** `SimulatorInterferometer` (no PSF; uv-plane noise instead of image-plane Poisson noise). +- **Galaxies:** Defining the galaxy with a Sersic bulge and an operated Gaussian point source. +- **Output:** Saving the simulated dataset to FITS files. +- **Visualize:** Outputting subplot and dirty-image visualizations as PNG files. +- **Galaxies json:** Saving the Galaxies object as a JSON file for future reference. +""" + +# from autogalaxy import setup_notebook; setup_notebook() + +from pathlib import Path +import autogalaxy as ag +import autogalaxy.plot as aplt + +""" +__Dataset Paths__ + +The `dataset_type` describes the type of data being simulated and `dataset_name` gives it a descriptive name. +""" +dataset_type = "interferometer" +dataset_name = "operated" +dataset_path = Path("dataset", dataset_type, dataset_name) + +""" +__Grid__ + +Simulate the image using a (y,x) grid. Over-sampling is an imaging-only technique and is not used for +interferometer data. +""" +grid = ag.Grid2D.uniform(shape_native=(256, 256), pixel_scales=0.1) + +""" +__uv-wavelengths__ + +To perform the Fourier transform we need the wavelengths of the baselines. +""" +uv_wavelengths_path = Path("dataset", dataset_type, "uv_wavelengths") +uv_wavelengths = ag.ndarray_via_fits_from( + file_path=Path(uv_wavelengths_path, "sma.fits"), hdu=0 +) + +""" +__Simulator__ + +Create the simulator for the interferometer data, which defines the exposure time, visibility-plane +noise sigma, and transformer. +""" +simulator = ag.SimulatorInterferometer( + uv_wavelengths=uv_wavelengths, + exposure_time=300.0, + noise_sigma=1000.0, + transformer_class=ag.TransformerDFT, +) + +""" +__Galaxies__ + +Setup the galaxy with a bulge (elliptical Sersic) for this simulation. + +This includes an operated `Gaussian` component which represents the compact point-source emission at the +galaxy's centre. Its image-plane shape is specified directly by the profile and is Fourier transformed to +the visibility plane like every other light profile (there is no PSF convolution for interferometer data). +""" +galaxy = ag.Galaxy( + redshift=0.5, + bulge=ag.lp.Sersic( + centre=(0.0, 0.0), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.9, angle=45.0), + intensity=1.0, + effective_radius=0.6, + sersic_index=3.0, + ), + psf=ag.lp_operated.Gaussian( + centre=(0.0, 0.0), ell_comps=(0.0, 0.0), intensity=100.0, sigma=0.1 + ), +) + +""" +Use these galaxies to generate the image which is simulated as an `Interferometer` dataset. +""" +galaxies = ag.Galaxies(galaxies=[galaxy]) + +aplt.plot_array(array=galaxies.image_2d_from(grid=grid), title="Image") + +""" +Pass the simulator galaxies, which creates the real-space image and NUFFTs it to visibilities. +""" +dataset = simulator.via_galaxies_from(galaxies=galaxies, grid=grid) + +""" +Plot the simulated `Interferometer` dataset before outputting it to fits. +""" +aplt.subplot_interferometer_dirty_images(dataset=dataset) + +""" +__Output__ + +Output the simulated dataset to the dataset path as .fits files. +""" +aplt.fits_interferometer( + dataset=dataset, + data_path=dataset_path / "data.fits", + noise_map_path=dataset_path / "noise_map.fits", + uv_wavelengths_path=dataset_path / "uv_wavelengths.fits", + overwrite=True, +) + +""" +__Visualize__ + +Output a subplot of the simulated dataset and the galaxies' images to the dataset path as .png files. +""" +aplt.subplot_interferometer_dirty_images( + dataset=dataset, output_path=dataset_path, output_format="png" +) +aplt.subplot_galaxies( + galaxies=galaxies, grid=grid, output_path=dataset_path, output_format="png" +) + +""" +__Galaxies json__ + +Save the `Galaxies` in the dataset folder as a .json file, ensuring the true light profiles and galaxies +are safely stored and available to check how the dataset was simulated in the future. + +This can be loaded via the method `galaxies = ag.from_json()`. +""" +ag.output_to_json( + obj=galaxies, + file_path=Path(dataset_path, "galaxies.json"), +) + +""" +The dataset can be viewed in the folder `autogalaxy_workspace/dataset/interferometer/operated`. +""" diff --git a/workspace_index.json b/workspace_index.json index 3b704653..357344af 100644 --- a/workspace_index.json +++ b/workspace_index.json @@ -2187,6 +2187,53 @@ "summary": "A multi-Gaussian expansion (MGE) decomposes a galaxy's light into ~15-100 Gaussians, where the `intensity` of every Gaussian is solved for via linear algebra using a process called an \"inversion\" (see the `linear_light_profiles` feature for a full description of this).", "title": "Modeling Features: Multi Gaussian Expansion (Interferometer)" }, + { + "contents": [ + "Mask", + "Dataset", + "Dataset Auto-Simulation", + "Over Sampling", + "Model", + "Search", + "Analysis", + "VRAM", + "Run Time", + "Model-Fit", + "Result" + ], + "cross_refs": [ + "/guides/results/start_here.py", + "/imaging/features/operated_light_profile/modeling.py", + "fit.py", + "interferometer/fit.py", + "interferometer/modeling.py", + "interferometer/start_here.ipynb" + ], + "notebook": "notebooks/interferometer/features/operated_light_profile/modeling.ipynb", + "path": "scripts/interferometer/features/operated_light_profile/modeling.py", + "summary": "It is common for galaxies to have point-source emission, for example bright emission right at their centre due to an active galactic nuclei or a very compact knot of star formation.", + "title": "Modeling Features: Operated Light Profiles (Interferometer)" + }, + { + "contents": [ + "Dataset Paths", + "Grid", + "uv-wavelengths", + "Simulator", + "Galaxies", + "Output", + "Visualize", + "Galaxies json" + ], + "cross_refs": [ + "interferometer/simulator.ipynb", + "modeling.py" + ], + "notebook": "notebooks/interferometer/features/operated_light_profile/simulator.ipynb", + "path": "scripts/interferometer/features/operated_light_profile/simulator.py", + "summary": "This script simulates `Interferometer` data of a galaxy using light profiles where:", + "title": "Simulator: Operated Light Profiles (Interferometer)" + }, { "contents": [ "Advantages",