notebooks/overview/overview_6_interferometer.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Overview: Interferometer\n",
        "------------------------\n",
        "\n",
        "**PyAutoLens** supports modeling of interferometer data from submillimeter and radio observatories such as ALMA\n",
        "or LOFAR.\n",
        "\n",
        "The visibilities of the interferometer dataset are fitted directly in the uv-plane, circumventing issues that arise\n",
        "when fitting a dirty image produced via the visibilities.\n",
        "\n",
        "The most important issue this addresses is removing correlated noise from impacting the fit."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "%matplotlib inline\n",
        "from pyprojroot import here\n",
        "workspace_path = str(here())\n",
        "%cd $workspace_path\n",
        "print(f\"Working Directory has been set to `{workspace_path}`\")\n",
        "\n",
        "import autofit as af\n",
        "import autolens as al\n",
        "import autolens.plot as aplt\n",
        "from os import path\n",
        "import numpy as np"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Real Space Mask__\n",
        "\n",
        "To begin, we define a real-space mask. Although interferometer lens modeling is performed in the uv-plane and \n",
        "therefore Fourier space, we still need to define the grid of coordinates in real-space from which the lensed source's \n",
        "images are computed. It is this image that is mapped to Fourier space to compare to the uv-plane data.\n",
        "\n",
        "The size and resolution of this mask depend on the baselines of your interferometer dataset. datasets with longer\n",
        "baselines (i.e. higher resolution data) require higher resolution and larger masks."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "real_space_mask = al.Mask2D.circular(\n",
        "    shape_native=(400, 400), pixel_scales=0.025, radius=3.0\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Dataset__\n",
        "\n",
        "We next load an interferometer dataset from fits files, which follows the same API that we have seen for an `Imaging`\n",
        "object."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dataset_name = \"simple\"\n",
        "dataset_path = path.join(\"dataset\", \"interferometer\", dataset_name)\n",
        "\n",
        "dataset = al.Interferometer.from_fits(\n",
        "    data_path=path.join(dataset_path, \"data.fits\"),\n",
        "    noise_map_path=path.join(dataset_path, \"noise_map.fits\"),\n",
        "    uv_wavelengths_path=path.join(dataset_path, \"uv_wavelengths.fits\"),\n",
        "    real_space_mask=real_space_mask,\n",
        "    transformer_class=al.TransformerNUFFT,\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The **PyAutoLens** plot module has tools for plotting interferometer datasets, including the visibilities, noise-map\n",
        "and uv wavelength which represent the interferometer`s baselines. \n",
        "\n",
        "The data used in this tutorial contains only ~300 visibilities and is representative of a low resolution \n",
        "Square-Mile Array (SMA) dataset. \n",
        "\n",
        "We made this choice so the script runs fast, and we discuss below how **PyAutoLens** can scale up to large visibilities\n",
        "datasets from an instrument like ALMA."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dataset_plotter = aplt.InterferometerPlotter(dataset=dataset)\n",
        "dataset_plotter.figures_2d(data=True, uv_wavelengths=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "It can also plot dirty images of the dataset in real-space, using the fast Fourier transforms described below."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dataset_plotter = aplt.InterferometerPlotter(dataset=dataset)\n",
        "dataset_plotter.figures_2d(dirty_image=True, dirty_signal_to_noise_map=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Tracer__\n",
        "\n",
        "To perform uv-plane modeling, **PyAutoLens** generates an image of the strong lens system in real-space via a tracer.\n",
        "\n",
        "Lets quickly set up the `Tracer` we'll use in this example."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "lens_galaxy = al.Galaxy(\n",
        "    redshift=0.5,\n",
        "    mass=al.mp.Isothermal(\n",
        "        centre=(0.0, 0.0),\n",
        "        einstein_radius=1.6,\n",
        "        ell_comps=al.convert.ell_comps_from(axis_ratio=0.9, angle=45.0),\n",
        "    ),\n",
        "    shear=al.mp.ExternalShear(gamma_1=0.05, gamma_2=0.05),\n",
        ")\n",
        "\n",
        "source_galaxy = al.Galaxy(\n",
        "    redshift=1.0,\n",
        "    bulge=al.lp.Sersic(\n",
        "        centre=(0.0, 0.0),\n",
        "        ell_comps=al.convert.ell_comps_from(axis_ratio=0.8, angle=60.0),\n",
        "        intensity=0.3,\n",
        "        effective_radius=1.0,\n",
        "        sersic_index=2.5,\n",
        "    ),\n",
        ")\n",
        "\n",
        "tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])\n",
        "\n",
        "tracer_plotter = aplt.TracerPlotter(\n",
        "    tracer=tracer, grid=real_space_mask.derive_grid.unmasked\n",
        ")\n",
        "tracer_plotter.figures_2d(image=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__UV-Plane__\n",
        "\n",
        "To perform uv-plane modeling, **PyAutoLens** next Fourier transforms this image from real-space to the uv-plane.\n",
        "\n",
        "This operation uses a `Transformer` object, of which there are multiple available in **PyAutoLens**. This includes\n",
        "a direct Fourier transform which performs the exact Fourier transform without approximation."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "transformer_class = al.TransformerDFT"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "However, the direct Fourier transform is inefficient. For ~10 million visibilities, it requires thousands of seconds\n",
        "to perform a single transform. This approach is therefore unfeasible for high quality ALMA and radio datasets.\n",
        "\n",
        "For this reason, **PyAutoLens** supports the non-uniform fast fourier transform algorithm\n",
        "**PyNUFFT** (https://github.com/jyhmiinlin/pynufft), which is significantly faster, being able to perform a Fourier\n",
        "transform of ~10 million in less than a second!\n",
        "\n",
        "This was passed input the `from_fits` method when loading the dataset above."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "transformer_class = al.TransformerNUFFT"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Fitting__\n",
        "\n",
        "The interferometer can now be used with a `FitInterferometer` object to fit it to a dataset:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "fit = al.FitInterferometer(dataset=dataset, tracer=tracer)\n",
        "\n",
        "fit_plotter = aplt.FitInterferometerPlotter(fit=fit)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualization of the fit is provided both in the uv-plane and in real-space. \n",
        "\n",
        "Note that the fit is not performed in real-space, but plotting it in real-space is often more informative."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "fit_plotter.subplot_fit()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Pixelized Sources__\n",
        "\n",
        "Interferometer data can also be modeled using pixelized source's, which again performs the source reconstruction by\n",
        "directly fitting the visibilities in the uv-plane."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "pixelization = al.Pixelization(\n",
        "    image_mesh=al.image_mesh.Overlay(shape=(30, 30)),\n",
        "    mesh=al.mesh.Delaunay(),\n",
        "    regularization=al.reg.Constant(coefficient=1.0),\n",
        ")\n",
        "\n",
        "source_galaxy = al.Galaxy(redshift=1.0, pixelization=pixelization)\n",
        "\n",
        "tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])\n",
        "\n",
        "fit = al.FitInterferometer(\n",
        "    dataset=dataset,\n",
        "    tracer=tracer,\n",
        "    settings_inversion=al.SettingsInversion(use_linear_operators=True),\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The source reconstruction is visualized in real space:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "fit_plotter = aplt.FitInterferometerPlotter(fit=fit)\n",
        "fit_plotter.figures_2d_of_planes(plane_index=1, plane_image=True)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Efficiency__\n",
        "\n",
        "Computing this source reconstruction would be extremely inefficient if **PyAutoLens** used a traditional approach to\n",
        "linear algebra which explicitly stored in memory the values required to solve for the source fluxes. In fact, for an\n",
        "interferometer dataset of ~10 million visibilities this would require **hundreds of GB of memory**!\n",
        "\n",
        "**PyAutoLens** uses the library **PyLops** (https://pylops.readthedocs.io/en/latest/) to represent this calculation as\n",
        "a sequence of memory-light linear operators.\n",
        "\n",
        "The combination of **PyNUFFT** and **PyLops** makes the analysis of ~10 million visibilities from observatories such as\n",
        "ALMA and JVLA feasible in **PyAutoLens**.\n",
        "\"\n",
        "__Lens Modeling__\n",
        "\n",
        "It is straight forward to fit a lens model to an interferometer dataset, using the same API that we saw for imaging\n",
        "data.\n",
        "\n",
        "We first compose the model, omitted the lens light components given that most strong lenses observed at submm /\n",
        "radio wavelengths do not have visible lens galaxy emission."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Lens:\n",
        "\n",
        "mass = af.Model(al.mp.Isothermal)\n",
        "\n",
        "lens = af.Model(al.Galaxy, redshift=0.5, mass=mass)\n",
        "\n",
        "# Source:\n",
        "\n",
        "disk = af.Model(al.lp.Exponential)\n",
        "\n",
        "source = af.Model(al.Galaxy, redshift=1.0, disk=disk)\n",
        "\n",
        "# Overall Lens Model:\n",
        "\n",
        "model = af.Collection(galaxies=af.Collection(lens=lens, source=source))"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Non-linear Search__\n",
        "\n",
        "We again choose the non-linear search `Nautilus` (https://github.com/joshspeagle/Nautilus)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "search = af.Nautilus(path_prefix=\"overview\", name=\"interferometer\")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Analysis__\n",
        "\n",
        "Whereas we previously used an `AnalysisImaging` object, we instead use an `AnalysisInterferometer` object which fits \n",
        "the lens model in the correct way for an interferometer dataset. \n",
        "\n",
        "This includes mapping the lens model from real-space to the uv-plane via the Fourier transform discussed above."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "analysis = al.AnalysisInterferometer(dataset=dataset)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Model-Fit__\n",
        "\n",
        "We can now begin the model-fit by passing the model and analysis object to the search, which performs a non-linear\n",
        "search to find which models fit the data with the highest likelihood.\n",
        "\n",
        "The results can be found in the `output/overview_interferometer` folder in the `autolens_workspace`."
      ]
    },
    {
      "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 **PyAutoLens** visualization library and `FitInterferometer` object includes specific methods for plotting the \n",
        "results, for example the maximum log likelihood fit:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "fit_plotter = aplt.FitInterferometerPlotter(fit=result.max_log_likelihood_fit)\n",
        "fit_plotter.subplot_fit()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Simulation__\n",
        "\n",
        "Simulated interferometer datasets can be generated using the ``SimulatorInterferometer`` object, which includes adding\n",
        "Gaussian noise to the visibilities:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "simulator = al.SimulatorInterferometer(\n",
        "    uv_wavelengths=dataset.uv_wavelengths, exposure_time=300.0, noise_sigma=0.01\n",
        ")\n",
        "\n",
        "real_space_grid = al.Grid2D.uniform(\n",
        "    shape_native=real_space_mask.shape_native,\n",
        "    pixel_scales=real_space_mask.pixel_scales,\n",
        ")\n",
        "\n",
        "dataset = simulator.via_tracer_from(tracer=tracer, grid=real_space_grid)\n",
        "\n",
        "dataset_plotter = aplt.InterferometerPlotter(dataset=dataset)\n",
        "dataset_plotter.subplot_dataset()"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Wrap Up__\n",
        "\n",
        "The `interferometer` package of the `autolens_workspace` contains numerous example scripts for performing \n",
        "interferometer modeling and simulating strong lens interferometer datasets."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    }
  ],
  "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
}