diff --git a/notebooks/NIRCam/NIRCam_wisp_subtraction/nircam_wisp_subtraction.ipynb b/notebooks/NIRCam/NIRCam_wisp_subtraction/nircam_wisp_subtraction.ipynb new file mode 100644 index 000000000..573968867 --- /dev/null +++ b/notebooks/NIRCam/NIRCam_wisp_subtraction/nircam_wisp_subtraction.ipynb @@ -0,0 +1,553 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "00417c2b-3e15-4ebd-bd0a-43307a5e1708", + "metadata": {}, + "source": [ + "# NIRCam Wisp Removal\n", + "---\n", + "**Author**: Ben Sunnquist (bsunnquist@stsci.edu)
\n", + "**Latest Update**: 24 July 2024
\n", + "**Use case**: NIRCam Imaging detectors A3, A4, B3, and B4.
\n", + "**Data**: PID 1063 Obs 6 Imaging Flats
\n", + "**Test Pipeline Version**: 1.15.1" + ] + }, + { + "cell_type": "markdown", + "id": "338e0783-ee8e-4291-9574-755faaa73610", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "\n", + "## Introduction\n", + "\n", + "This notebook demonstrates how to remove \n", + "[wisps](https://jwst-docs.stsci.edu/known-issues-with-jwst-data/nircam-known-issues/nircam-scattered-light-artifacts#NIRCamScatteredLightArtifacts-wispsWisps) from NIRCam imaging data. Wisps are a scattered light feature affecting detectors A3, A4, B3, and B4. For a given filter, wisps appear in the same detector location with only their brightness varying between exposures; therefore, they can be removed from science data by scaling and subtracting a wisp template (i.e. a median combination of all wisp appearances).\n", + "\n", + "The wisp templates used in this notebook will be downloaded in the [Data](#data) section, but they are also available in the \"version3\" folder within the [NIRCam wisp template Box folder](https://stsci.box.com/s/1bymvf1lkrqbdn9rnkluzqk30e8o2bne). Future updates to these templates will be added to this same Box folder, and users are encouraged to check this folder regularly to ensure they have the latest version.\n", + "\n", + "This notebook uses the `subtract_wisp.py` code to scale and subtract the wisps. That code can be used by itself within python, and is preferred if calibrating a large number of files in parallel, but this notebook will be used to demonstrate the various parameters available to optimize wisp removal. For each notebook cell, we'll also show the corresponding command to run the equivalent in python." + ] + }, + { + "cell_type": "markdown", + "id": "3f29f4b4-f889-4cbb-89be-4dc5d8c7e5d1", + "metadata": {}, + "source": [ + "\n", + "## Table of Contents\n", + "* [Imports](#imports)\n", + "* [Data](#data)\n", + "* [Run the wisp subtraction code](#Run-the-wisp-subtraction-code)\n", + "* [Run the wisp subtraction code using custom settings](#Run-the-wisp-subtraction-code-using-custom-settings)\n", + " * [New wisp scaling settings](#Example-1:-New-wisp-scaling-method)\n", + " * [Data quality flagging](#Example-2:-Data-quality-flagging)\n", + " * [Other misc settings](#Example-3:-Other-misc-settings)\n", + "* [Using custom wisp templates](#Using-custom-wisp-templates)\n", + "* [Running the code in pieces](#Runing-the-code-in-pieces)" + ] + }, + { + "cell_type": "markdown", + "id": "0a824094-c56b-44b3-80e6-9466a6f4f57a", + "metadata": {}, + "source": [ + "\n", + "## Imports\n", + "\n", + "How to make an environment to run this notebook:\n", + "\n", + "\n", + "conda create -n jwst1140 python=3.11 notebook\n", + "source activate jwst1140\n", + "pip install -r requirements.txt\n", + "\n", + "\n", + "
\n", + "The version of the pipeline used by this notebook can be updated in the `requirments.txt` file before installation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fa348563-e007-4135-9481-7fed3e6c09e6", + "metadata": {}, + "outputs": [], + "source": [ + "# load functions from the subtract_wisp.py code \n", + "from subtract_wisp import make_segmap, process_file, subtract_wisp\n", + "\n", + "# load important packages\n", + "from astropy.io import fits\n", + "from astroquery.mast import Mast, Observations\n", + "import glob\n", + "import os\n", + "import urllib\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "id": "38e477be-b001-49be-9b8f-7fb503686199", + "metadata": {}, + "source": [ + "## Data\n", + "\n", + "We'll use F200W imaging data from program 1063 observation 6 taken on\n", + "detectors A3, A4, B3, and B4. These are 6 images of sparse fields that\n", + "are affected by bright wisps. We'll also use the corresponding\n", + "longwave images for source detection since the longwave is not impacted by wisps.\n", + "\n", + "Download the relevant data from MAST using\n", + "[astroquery](https://astroquery.readthedocs.io/en/latest/mast/mast.html).
\n", + "\n", + "
\n", + "Warning: If this notebook is halted during this step the downloaded\n", + "file may be incomplete, and cause crashes later on!\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "499b9d8d-1f8a-4258-925d-3b60d44c7de7", + "metadata": {}, + "outputs": [], + "source": [ + "# Select the data relevant data from MAST\n", + "params = {\"columns\": \"dataURI, filename, exp_type\",\n", + " \"filters\": [{\"paramName\": \"program\", \"values\": ['1063']},\n", + " {\"paramName\": \"observtn\", \"values\": ['6']},\n", + " {\"paramName\": \"exposure\", \"values\": ['00005']},\n", + " {\"paramName\": \"visit\", \"values\": ['004']},\n", + " {\"paramName\": \"detector\", \"values\": ['NRCA3', 'NRCA4', 'NRCB3', 'NRCB4', 'NRCALONG', 'NRCBLONG']},\n", + " {\"paramName\": \"productLevel\", \"values\": ['2b']}]}\n", + "t = Mast().service_request('Mast.Jwst.Filtered.Nircam', params)\n", + "for row in t:\n", + " if '_cal' in row['filename']: # only want cal files\n", + " result = Observations().download_file(row['dataURI'], cache=False)" + ] + }, + { + "cell_type": "markdown", + "id": "ff808e6b-ddbf-4e58-ac0b-26b1bf24606a", + "metadata": {}, + "source": [ + "Next, we'll download the wisp templates. These templates are equivalent to those in the \"version3\" folder within the [NIRCam wisp template Box folder](https://stsci.box.com/s/1bymvf1lkrqbdn9rnkluzqk30e8o2bne). Users are encouraged to check this folder regularly for any updates. For this notebook, only the F200W templates are needed:\n", + "\n", + "* WISP_NRCA3_F200W_CLEAR.fits\n", + "* WISP_NRCA4_F200W_CLEAR.fits\n", + "* WISP_NRCB3_F200W_CLEAR.fits\n", + "* WISP_NRCB4_F200W_CLEAR.fits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef409ef8-280c-4156-bd78-f1e06e5820e9", + "metadata": {}, + "outputs": [], + "source": [ + "# Download the wisp templates\n", + "boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/nircam_wisp_templates/'\n", + "\n", + "for detector in ['NRCA3', 'NRCA4', 'NRCB3', 'NRCB4']:\n", + " boxfile = os.path.join(boxlink, 'WISP_{}_F200W_CLEAR.fits'.format(detector))\n", + " urllib.request.urlretrieve(boxfile, os.path.basename(boxfile))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ad7f9846-e6b2-4d94-9d96-94057487c980", + "metadata": {}, + "outputs": [], + "source": [ + "# Check if template files are downloaded\n", + "template_files = ['WISP_NRCA3_F200W_CLEAR.fits', 'WISP_NRCA4_F200W_CLEAR.fits',\n", + " 'WISP_NRCB3_F200W_CLEAR.fits', 'WISP_NRCB4_F200W_CLEAR.fits']\n", + "\n", + "for tem_file in template_files:\n", + " if not os.path.isfile(tem_file):\n", + " print(f'file {tem_file} does not exist, please downoad it from Box')" + ] + }, + { + "cell_type": "markdown", + "id": "32af21b1-b980-433c-907a-578d34629339", + "metadata": {}, + "source": [ + "## Run the wisp subtraction code\n", + "\n", + "In this section, we show how to run the wisp subtraction code on all of the shortwave images using the default parameters. The exception is we set `show_plot=True` to allow the diagnostic plots to display in the notebook.\n", + "\n", + "The main processing function `process_file()` combines the segmap creation and wisp scaling/subtraction steps together. The basic process the code uses is the following:\n", + "\n", + "1. Generate a source mask using the corresponding longwave image, since the longwave isn't affected by wisps. Ignore these pixels going forward. Occurs in the `make_segmap()` function.\n", + "2. Load the relevant wisp template.\n", + "3. Apply a series of scale factors to the template. For each scale factor, multiply it onto the wisp template and subtract the result from the input science data. Record the residual noise within the wisp region for each of the scale factors tested. A correction to residual 1/f (horizontal) noise is applied during this step to help with the fitting. Fit a polynomial to the residuals, and choose the scale factor with the lowest noise based on this fit. Occurs in the `subtract_wisp()` function.\n", + "4. Multiply the chosen wisp factor onto the wisp template to generate the final wisp model. Occurs in the `subtract_wisp()` function.\n", + "5. Subtract the final wisp model from the input data, and write out the wisp-subtracted data, the wisp model, and diagnostic plots summarizing the results. Occurs in the `subtract_wisp()` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24f58c3a-1f7b-4fdc-aa7f-fbe824c2358f", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Collect all the science data files \n", + "files = glob.glob('*_cal.fits')\n", + "\n", + "files = [cal_file for cal_file in files if 'long' not in cal_file] # only want shortwave files" + ] + }, + { + "cell_type": "markdown", + "id": "3b16a76b-75cf-4861-81f0-6b249d0533f7", + "metadata": {}, + "source": [ + "Run the `process_file()` function in each dataset. The code generates the following:\n", + "\n", + "* A wisp-corrected version of the input file, with the same name as the input file besides a _wisp suffix added (note: the suffix used here can be changed using the `suffix` argument, which can also be set to an empty string to overwrite the input file).\n", + "* The wisp model subtracted from the input image, with the same name as the input file besides a _wisp_model suffix added.\n", + "* A diagnostic plot showing the original image, the wisp-corrected image, the wisp model subtracted, and a plot showing all of the wisp scale factors tested along with their corresponding residual noise. The factor with the lowest noise is the one applied to the wisp template to generate the wisp model.\n", + "\n", + "Note: calling `process_files()` rather than `process_file()` will process the files in parallel. This will save time, but the plots won't display inline in the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c196dbd7-902a-42b1-9f85-5402b5220fcb", + "metadata": {}, + "outputs": [], + "source": [ + "# Remove wisps from the science data\n", + "for file in files:\n", + " results = process_file(file, show_plot=True)" + ] + }, + { + "cell_type": "markdown", + "id": "7e5c7f06-a38f-4a37-85d1-5a4cdbc24dc5", + "metadata": {}, + "source": [ + "Python equivalent: \n", + "`python subtract_wisp.py --files *_cal.fits`" + ] + }, + { + "cell_type": "markdown", + "id": "eb65ba67-65c1-4623-ae03-a3267f0e2d37", + "metadata": {}, + "source": [ + "\n", + "## Run the wisp subtraction code using custom settings" + ] + }, + { + "cell_type": "markdown", + "id": "3c9cd47f-7d44-4817-b2a4-e754db0e0888", + "metadata": {}, + "source": [ + "In this section, we'll run the wisp subtraction code using custom settings that may be helpful to optimize the wisp removal. The different settings shown here can be used in combination with those from other examples or as a standalone." + ] + }, + { + "cell_type": "markdown", + "id": "996ff57a-cd4e-44e5-8f74-9fe21dd90023", + "metadata": {}, + "source": [ + "\n", + "### Example 1: New wisp scaling method" + ] + }, + { + "cell_type": "markdown", + "id": "752133eb-2271-4e22-b1c1-8a605f8ad807", + "metadata": {}, + "source": [ + "By default, the wisp factor is decided based on that which results in the lowest residual noise. However, the scaling can instead be based on minimizing the overall signal level differences inside and outside the wisp-affected detector regions by setting `scale_method='median'`. For this scaling method, it's recommended to also turn off the polynomial fitting on the residuals (`poly_degree=0`) to simply choose the factor that results in the lowest difference.\n", + "\n", + "During the scaling, residual 1/f noise is corrected by default (`correct_rows=True`); however, for datasets that have strong e.g. amp offsets or odd-even column residuals, this same correction can be applied in the vertical direction as well (`correct_cols=True`). Note that these corrections aren't subtracted in the output file, they are only used to help the wisp scaling procedure.\n", + "\n", + "Below, we perform a new run using these settings." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "67c3d9ae-c62e-4eb7-a550-dc181183c159", + "metadata": {}, + "outputs": [], + "source": [ + "# Remove wisps using the scaling method\n", + "results = process_file('jw01063006004_02101_00005_nrcb4_cal.fits',\n", + " scale_method='median', poly_degree=0, \n", + " correct_cols=True, show_plot=True)" + ] + }, + { + "cell_type": "markdown", + "id": "e5f0466b-afed-4df9-a82f-bb6cb39eb9f0", + "metadata": {}, + "source": [ + "Python equivalent: `python subtract_wisp.py --files jw01063006004_02101_00005_nrcb4_cal.fits --scale_method median --poly_degree 0 --correct_cols`" + ] + }, + { + "cell_type": "markdown", + "id": "e69a261c-00ed-4b78-b5a0-fa1a38ee736d", + "metadata": {}, + "source": [ + "Sometimes, e.g. for crowded fields, there may not be enough empty regions to scale the wisp properly. In this case, simply subtracting the wisp template with no scaling may provide the best results (`scale_wisp=False`, or similarly, in python `--no-scale_wisp`)." + ] + }, + { + "cell_type": "markdown", + "id": "2bbf490a-45f4-4475-b9dd-aa5c89408126", + "metadata": {}, + "source": [ + "\n", + "### Example 2: Data quality flagging" + ] + }, + { + "cell_type": "markdown", + "id": "e1817da3-aa13-471c-aa3f-ff2af64e0f10", + "metadata": {}, + "source": [ + "Another option is to flag any pixels with wisp signals above a certain threshold in the image's data quality array, but not subtract the wisp itself from the science data. This may be useful if you want to ignore these pixels when combining several images into a mosaic or flag them while performing photometry. This can be accomplished by setting `flag_wisp_thresh` to the chosen wisp signal threshold (units are assumed to be the same as the input file) and `sub_wisp=False`. The data quality value used here is set by the `dq_val` parameter.\n", + "\n", + "The default value of `dq_val` is 1, which is equivalent to the JWST Pipeline DO_NOT_USE flag, and it is ignored by default during image3 drizzling. If you would prefer the pixels to be flagged, but not ignored during that step, another option is using `dq_val=1073741824`, which is equivalent to the JWST Pipeine OTHER_BAD_PIXEL flag. More details on the various data quality flags can be found [here](https://jwst-pipeline.readthedocs.io/en/latest/jwst/references_general/references_general.html#data-quality-flags)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e69fcf4f-a289-48f4-88a3-21609e87b770", + "metadata": {}, + "outputs": [], + "source": [ + "# Do not remove wisps, just flag them\n", + "results = process_file('jw01063006004_02101_00005_nrcb4_cal.fits',\n", + " flag_wisp_thresh=0.03, dq_val=1073741824,\n", + " sub_wisp=False, show_plot=True)" + ] + }, + { + "cell_type": "markdown", + "id": "98429059-9038-40a3-be3c-ca5edb81b0d8", + "metadata": {}, + "source": [ + "Python equivalent: `python subtract_wisp.py --files jw01063006004_02101_00005_nrcb4_cal.fits --flag_wisp_thresh 0.03 --dq_val 1073741824 --no-sub_wisp`\n" + ] + }, + { + "cell_type": "markdown", + "id": "20506cbb-6f6a-41a0-b6cf-4492b7ce4859", + "metadata": {}, + "source": [ + "Check that the data quality array in this file was updated appropriately" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e5736e6e-1e40-4e60-8490-c5740a94d13d", + "metadata": {}, + "outputs": [], + "source": [ + "# Check data quality (DQ) array \n", + "dq = fits.getdata('jw01063006004_02101_00005_nrcb4_cal_wisp.fits', 'DQ')\n", + "\n", + "# Show pixels flagged as OTHER_BAD_PIXEL\n", + "dq = (dq & 1073741824 != 0).astype(int) # only want to see pixels with dq_val 1073741824\n", + "plt.imshow(dq, cmap='gray', origin='lower', vmin=0, vmax=0.1)" + ] + }, + { + "cell_type": "markdown", + "id": "275c048c-f642-46c0-9f0b-68b9bcda0d53", + "metadata": {}, + "source": [ + "\n", + "### Example 3: Other misc settings" + ] + }, + { + "cell_type": "markdown", + "id": "6fbfb5c2-0b40-4fa4-8890-3fdb1a66e284", + "metadata": {}, + "source": [ + "This example will show a variety of other settings that may be\n", + "useful either used alone or in combination with any of the\n", + "other settings shown above. The full list of optional settings\n", + "can be seen in the `make_segmap()` and `subtract_wisp()` docstrings.\n", + "\n", + "In this example, we'll do the source-finding on the input image\n", + "itself (rather than the corresponding longwave image) by\n", + "setting `seg_from_lw=False`. We'll increase the source\n", + "detection `sigma` to help avoid flagging the wisp itself, and\n", + "we'll save the resulting segmentation mask by setting\n", + "`save_segmap=True` for inspection. The segmentation mask\n", + "filename is the same as the input file, just with the\n", + "_seg suffix added.\n", + "\n", + "We'll also only subtract wisp signals above 0.01 MJy/sr\n", + "(`min_wisp=0.01`). To save a little time, we'll only test wispscale factors between 0.5 and 1.5 in 0.05 steps." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8eaf492c-bff4-48e9-9220-e14e56e43e9d", + "metadata": {}, + "outputs": [], + "source": [ + "# Get full list of optional settings for make_segmap() using help(make_segmap)\n", + "#\n", + "# Use special settings\n", + "results = process_file('jw01063006004_02101_00005_nrcb4_cal.fits',\n", + " seg_from_lw=False, sigma=1.5, save_segmap=True, \n", + " min_wisp=0.01, factor_min=0.5, factor_max=1.5,\n", + " factor_step=0.05, show_plot=True)" + ] + }, + { + "cell_type": "markdown", + "id": "8195caf8-ca91-497a-b904-e54d4a990e8d", + "metadata": {}, + "source": [ + "Python equivalent: `python subtract_wisp.py --files jw01063006004_02101_00005_nrcb4_cal.fits --no-seg_from_lw --sigma 1.5 --save_segmap --min_wisp 0.01 --factor_min 0.5 --factor_max 1.5 --factor_step 0.05`\n" + ] + }, + { + "cell_type": "markdown", + "id": "793c22b0-9ad0-4692-91bb-d96792696ab1", + "metadata": {}, + "source": [ + "As seen above, these custom settings resulted in an oversubtraction. This is likely due to poor source finding in the wisp region (e.g. the wisp itself being flagged as a source), and why source finding using the longwave image is recommended. We can confirm this by inspecting the generated segmentation map." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "acaaaaf0-9285-42ee-a6b8-3bf84f4cc826", + "metadata": {}, + "outputs": [], + "source": [ + "data = fits.getdata('jw01063006004_02101_00005_nrcb4_cal.fits')\n", + "segmap = fits.getdata('jw01063006004_02101_00005_nrcb4_cal_seg.fits')\n", + "\n", + "# Show the images\n", + "fig, axes = plt.subplots(1, 2, figsize=(20, 10))\n", + "axes[0].imshow(data, origin='lower', cmap='gray', vmin=0.18, vmax=0.3)\n", + "axes[1].imshow(segmap, origin='lower', cmap='gray', vmin=0, vmax=0.1)" + ] + }, + { + "cell_type": "markdown", + "id": "b57fa256-94df-4577-a8d7-404e2824a1f4", + "metadata": {}, + "source": [ + "\n", + "## Using custom wisp templates" + ] + }, + { + "cell_type": "markdown", + "id": "a01cb24d-13af-4745-b89b-3ccc7d0e1dc4", + "metadata": {}, + "source": [ + "This section shows how you can incorporate custom wisp templates while still using the wisp subtraction procedure contained in the main code. This may be useful if you have your own wisp templates or want to apply a manual scale factor to the existing templates.\n", + "\n", + "In this example, we'll use the existing template, but multiply it by our own scale factor and turn the scaling procedure off (`scale_wisp=False`) so the code simply applies the wisp model as given. Note that turning off the scaling procedure is not required when using a custom template; users can both provide a custom template and still keep the scaling procedure on. Since the scaling procedure is off in this example, the source finding is irrelevant, and so we can skip the segmap creation function `make_segmap()` and run `subtract_wisp()` directly." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17dbbd90-a406-416e-a70d-2fc45b0a7b38", + "metadata": {}, + "outputs": [], + "source": [ + "# Get the wisp template and apply a manual scale factor\n", + "wisp_data = fits.getdata('WISP_NRCB4_F200W_CLEAR.fits') * 1.05\n", + "\n", + "# Process the file with the custom wisp template\n", + "results = subtract_wisp('jw01063006004_02101_00005_nrcb4_cal.fits', wisp_data, scale_wisp=False, show_plot=True)" + ] + }, + { + "cell_type": "markdown", + "id": "5aedb011-4c13-49ce-bd18-c5d60417d0fc", + "metadata": {}, + "source": [ + "\n", + "## Runing the code in pieces" + ] + }, + { + "cell_type": "markdown", + "id": "a712be10-d672-4d67-be18-b90a94393b34", + "metadata": {}, + "source": [ + "In this section, we'll show how to run the full wisp removal process in pieces. This may be useful to provide more flexibility or when finding the optimal settings to use for your data. We'll use a variety of custom settings - see the function docstrings for more details." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aaf41481-d1f1-4edc-9b8a-9e2d0f06527f", + "metadata": {}, + "outputs": [], + "source": [ + "# The file to process\n", + "cal_file = 'jw01063006004_02101_00005_nrcb4_cal.fits'\n", + "\n", + "# The wisp template to use. We'll use a \"custom\" template again.\n", + "wisp_data = fits.getdata('WISP_NRCB4_F200W_CLEAR.fits') * 1.05\n", + "\n", + "# Generate the segmentation map\n", + "segmap_data = make_segmap(cal_file, sigma=1.0, npixels=8, dilate_segmap=10)\n", + "\n", + "# Scale and subtract the wisp template\n", + "results = subtract_wisp(cal_file, wisp_data=wisp_data, segmap_data=segmap_data, show_plot=True)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/NIRCam/NIRCam_wisp_subtraction/requirements.txt b/notebooks/NIRCam/NIRCam_wisp_subtraction/requirements.txt new file mode 100644 index 000000000..161385ab4 --- /dev/null +++ b/notebooks/NIRCam/NIRCam_wisp_subtraction/requirements.txt @@ -0,0 +1,2 @@ +astroquery==0.4.7 +jwst==1.15.1 \ No newline at end of file diff --git a/notebooks/NIRCam/NIRCam_wisp_subtraction/subtract_wisp.py b/notebooks/NIRCam/NIRCam_wisp_subtraction/subtract_wisp.py new file mode 100644 index 000000000..a5840d31b --- /dev/null +++ b/notebooks/NIRCam/NIRCam_wisp_subtraction/subtract_wisp.py @@ -0,0 +1,570 @@ +#! /usr/bin/env python + +""" +Scales and subtracts wisp templates from full-frame NIRCam images. + +Authors +------- + Ben Sunnquist, 2024 + +Use +--- + To run from the command line using all default arguments (assumes input files + and wisp templates are in the current working directory): + >>> python subtract_wisp.py + + To run on all NRCB4 files with no wisp scaling or segmap creation, and flagging + wisp values above 0.01 in the input images' data quality arrays: + >>> python subtract_wisp.py --files ./*nrcb4_cal.fits --no-scale_wisp --no-create_segmap --flag_wisp_thresh 0.01 + + To run the same example above within python e.g. a Jupyter notebook: + >>> from subtract_wisp import process_files + >>> process_files(files='./*nrcb4_cal.fits', create_segmap=False, scale_wisp=False, flag_wisp_thresh=0.01) + +Notes +----- + Wisp templates designed to work with this code are available at the following STScI Box link: + https://stsci.app.box.com/s/1bymvf1lkrqbdn9rnkluzqk30e8o2bne + + This code was developed using the JWST Pipeline v1.14.0 environment: + https://github.com/spacetelescope/jwst + > conda create -n python=3.11 + > conda activate + > pip install jwst==1.14.0 + Other environments may work fine, but be sure to use Python 3.9+. +""" + +import argparse +from functools import partial +import multiprocessing +import os +import warnings + +from astropy.convolution import convolve, Gaussian2DKernel +from astropy.io import fits +from astropy.stats import median_absolute_deviation, sigma_clipped_stats +from astropy.wcs import WCS +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from photutils.segmentation import detect_sources, detect_threshold +from scipy.ndimage import binary_dilation, generate_binary_structure +matplotlib.use('Agg') + +warnings.filterwarnings('ignore', message="Input data contains invalid values*") # nan values expected throughout code +warnings.filterwarnings('ignore', message="All-NaN slice encountered*") +warnings.filterwarnings('ignore', message="'obsfix' made the change*") # from astropy wcs during lw segmap blotting +warnings.filterwarnings('ignore', message="'datfix' made the change*") # from astropy wcs during lw segmap blotting + + +# ----------------------------------------------------------------------------- + + +def make_segmap(f, seg_from_lw=True, sigma=0.8, npixels=10, dilate_segmap=5, save_segmap=False): + """ + Make a segmentation map for the input file. + + Parameters + ---------- + f : str + The filename to make a segmentation map for, e.g. + jw01063006004_02101_00005_nrca3_cal.fits. Source-finding is only performed + on cal files, so the cal file must be present if a rate file is input here. + + seg_from_lw : bool + Option to perform the source-finding on the corresponding + longwave image. This is generally preferred since the longwave is not + affected by wisps. + + sigma: float + The source detection sigma threshold. + + npixels: int + The number of connected pixels above sigma to include as a source. + + dilate_segmap: int + The number of pixels to grow the segmap outwards by. + + save_segmap : bool + Option to save the generated segmentation map. + + Returns + ------- + segmap : numpy.ndarray + The generated segmentation map. + + Outputs + ------- + {f}_seg.fits : fits image + The generated segmentation map. + """ + + print('Making segmap for {}'.format(f)) + outfile = f.replace('.fits', '_seg.fits') + + # Get the input data; always source-find on the cal image + detector = os.path.basename(f).split('_')[-2].lower() + f_sw = f.replace('_rate.fits', '_cal.fits') + if (seg_from_lw) & ('long' not in detector): + if 'a' in detector: + f_lw = f.replace(detector.lower(), 'nrcalong').replace('_rate.fits', '_cal.fits') + if 'b' in detector: + f_lw = f.replace(detector.lower(), 'nrcblong').replace('_rate.fits', '_cal.fits') + data = fits.getdata(f_lw, 'SCI') + dq = fits.getdata(f_lw, 'DQ') + else: + data = fits.getdata(f_sw, 'SCI') + dq = fits.getdata(f_sw, 'DQ') + + # Make the segmentation map + threshold = detect_threshold(data, sigma) + g = Gaussian2DKernel(x_stddev=3) + data_conv = convolve(data, g, mask=dq & 1 != 0) # Smooth input image before detecting sources + seg = detect_sources(data_conv, threshold, npixels=npixels, mask=dq & 1 != 0) # avoid bad pixels as sources + segmap_data = seg.data + segmap_data[segmap_data != 0] = 1 + + # Dilate the segmap outwards + if dilate_segmap != 0: + segmap_data = binary_dilation(segmap_data, iterations=dilate_segmap, structure=generate_binary_structure(2, 2)) + segmap_data[segmap_data != 0] = 1 + + # Blot LW segmap back onto SW detector space + if (seg_from_lw) & ('long' not in detector): + segmap_tmp = np.zeros(segmap_data.shape).astype(int) + wcs = WCS(fits.getheader(f_lw, 'SCI')) # lw cal wcs + seg_y, seg_x = np.where(segmap_data != 0) + sky_coords = wcs.pixel_to_world(seg_x, seg_y) + wcs = WCS(fits.getheader(f_sw, 'SCI')) # sw cal wcs + coords = wcs.world_to_pixel(sky_coords) + for i in np.arange(len(coords[0])): + y, x = int(coords[1][i]), int(coords[0][i]) + if (y < 2048) & (x < 2048) & (y >= 0) & (x >= 0): + segmap_tmp[y, x] = 1 + # Dilate to compensate for y,x rounding due to different lw/sw pixel scales + segmap_data = binary_dilation(segmap_tmp, iterations=1, structure=generate_binary_structure(2, 2)).astype(int) + + # Save the segmap + segmap_data = segmap_data.astype(int) + if save_segmap: + fits.writeto(outfile, segmap_data, overwrite=True) + + return segmap_data + + +# ----------------------------------------------------------------------------- + + +def process_file(f, wisp_dir='./', create_segmap=True, seg_from_lw=True, sigma=0.8, npixels=10, dilate_segmap=5, + save_segmap=False, sub_wisp=True, gauss_smooth_wisp=False, gauss_stddev=3.0, scale_wisp=True, + scale_method='mad', poly_degree=5, factor_min=0.0, factor_max=2.0, factor_step=0.01, min_wisp=None, + flag_wisp_thresh=None, dq_val=1, correct_rows=True, correct_cols=False, save_data=True, save_model=True, + plot=True, show_plot=False, suffix='_wisp'): + """ + The main processing function. Combines the segmap creation and wisp scaling/subtraction steps together. + + See make_segmap() and subtract_wisp() docstrings for more details on all of the input arguments. + + Parameters + ---------- + f : str + The input filename to subtract the wisp from. Must be a full-frame NIRCam image, e.g. + jw01063006004_02101_00005_nrca3_cal.fits. + + wisp_dir : str + The directory containing the wisp templates. The templates are assumed to have the + form WISP_{DETECTOR}_{FILTER}_{PUPIL}.fits. + """ + + # Get the relevant wisp template + print('Processing {}'.format(f)) + header = fits.getheader(f) + det, fltr, pupil = header['DETECTOR'], header['FILTER'], header['PUPIL'] + file_type = f.split('_')[-1].replace('.fits', '').upper() # CAL or RATE + wisp_data = fits.getdata(os.path.join(wisp_dir, 'WISP_{}_{}_{}.fits'.format(det, fltr, pupil)), file_type) + + # Make the segmentation map + if create_segmap: + segmap_data = make_segmap(f, seg_from_lw=seg_from_lw, sigma=sigma, npixels=npixels, + dilate_segmap=dilate_segmap, save_segmap=save_segmap) + else: + segmap_data = np.zeros(wisp_data.shape).astype(int) + + # Scale and subtract wisp template + _ = subtract_wisp(f, wisp_data=wisp_data, segmap_data=segmap_data, sub_wisp=sub_wisp, + gauss_smooth_wisp=gauss_smooth_wisp, gauss_stddev=gauss_stddev, + scale_wisp=scale_wisp, scale_method=scale_method, poly_degree=poly_degree, + factor_min=factor_min, factor_max=factor_max, factor_step=factor_step, + min_wisp=min_wisp, flag_wisp_thresh=flag_wisp_thresh, dq_val=dq_val, + correct_rows=correct_rows, correct_cols=correct_cols, save_data=save_data, + save_model=save_model, plot=plot, show_plot=show_plot, suffix=suffix) + print('Processing complete for {}'.format(f)) + +# ----------------------------------------------------------------------------- + + +def process_files(files, nproc=6, **kwargs): + """"Wrapper around the process_file() function to allow for multiprocessing.""" + + # Remove any files that are not in a detector impacted by wisps + files = [f for f in files if any(substring in f for substring in ['nrca3', 'nrca4', 'nrcb3', 'nrcb4'])] + print('Found {} relevant input files.'.format(len(files))) + + # Proess the files + process_file_partial = partial(process_file, **kwargs) + p = multiprocessing.Pool(nproc) + _ = p.map(process_file_partial, files) + p.close() + p.join() + + +# ----------------------------------------------------------------------------- + + +def subtract_wisp(f, wisp_data, segmap_data=None, sub_wisp=True, gauss_smooth_wisp=False, gauss_stddev=3.0, scale_wisp=True, + scale_method='mad', poly_degree=5, factor_min=0.0, factor_max=2.0, factor_step=0.01, min_wisp=None, + flag_wisp_thresh=None, dq_val=1, correct_rows=True, correct_cols=False, save_data=True, save_model=True, plot=True, + show_plot=False, suffix='_wisp'): + """Scales and subtracts a wisp template from the input file. + + Parameters + ---------- + f : str + The file to subtract the wisp from. Must be a full-frame NIRCam image, e.g. + jw01063006004_02101_00005_nrca3_cal.fits. + + wisp_data : numpy.ndarray + The wisp template data. Must be the same shape as input file data. + + segmap_data : numpy.ndarray + A mask of the sources in the input file. Must be the same shape as input file data. + + sub_wisp : bool + Option to subtract the wisp template from the input file. Generally this is always True + except for cases where you only want to flag high wisp values in the input file's data + quality array, i.e. using flag_wisp_thresh. + + gauss_smooth_wisp: bool + Option to smooth the wisp template with a Gaussian filter before scaling/subtracting it. + + gauss_stddev: float + The standard deviation of the Gaussian filter to apply to the wisp if gauss_smooth_wisp is True. + + scale_wisp : bool + Option to scale the wisp templates before subtracting it from the input file. + + scale_method : str + The method used to scale the wisp template. Options are 'mad' and 'median'. The former + chooses the scale factor by minimizing the median absolute deviation of the wisp-corrected + image. The latter chooses the scale factor by minimizing the overall signal level + difference between the wisp-affected region and the remaining detector region. + + poly_degree : int + The degree of the fitting polynomial applied to the residuals of the chosen scale_method. + The minimum of this polynomial fit decides the wisp scale factor. Set to 0 to ignore, + which is generally recommended when scale_method is 'median'. + + factor_min: float + The lowest factor to scale the wisp template by. + + factor_max: float + The highest factor to scale the wisp template by. + + factor_step: float + The step distance between factor_min and factor_max to try scaling the wisp template by. + Scale factors from factor_min to factor_max in iterations equal to the step distance + specified here will all be tested. + + min_wisp : float + The minimum wisp value to perform wisp subtraction. Everything below this value will be set + to zero in the wisp template before subtracting from the input file. The units of this + value should match the units in the input file. Set to None to ignore. + + flag_wisp_thresh : float + Flag values in the wisp template above this threshold in the input file's data quality + array. The units of this value should match the units in the input file. Set to None to ignore. + The DQ value used here is set by the dq_val parameter. + + dq_val : int + The data quality value applied to the input file's DQ array to pixels above flag_wisp_thresh. + The default of 1, i.e. DO_NOT_USE in the JWST Pipeline, results in these pixels being ignored + in the JWST Pipeline image3 drizzling step by default. 1073741824, i.e. OTHER_BAD_PIXEL, is another + option to use if you want these pixels flagged in the input file's DQ array, but still used in the + default drizzling step. Only relevant if flag_wisp_thresh is True. + + correct_rows : bool + Option to correct horizontal noise in the input file (e.g. 1/f noise residuals). This won't be + subtracted from the input file; it is only used to help the wisp scaling step. + + correct_cols : bool + Option to correct vertical noise in the input file (e.g. residual odd-even column effects + and amp offsets). This noise is generally much lower than the horizontal noise in correct_rows. + This won't be subtracted from the input file; it is only used to help the wisp scaling step. + + save_data : bool + Option to save the wisp-subtracted data. The filename will be {f}{suffix}.fits. + + save_model : bool + Option to save the wisp model subtracted fom the input file. The filename will + be {f}{suffix}_model.fits. + + plot : bool + Option to save a diagnostic plot showing the original data, wisp-subtracted data, final wisp model, + and a plot of the various wisp scale factors tested with their corresponding residual noise. + The filename will be {f}{suffix}_plots.png. + + show_plot : bool + Option to show the diagnostic plot in e.g. jupyter notebook (note you'll have to set + %matplotlib inline in the notebook as well for this to work). + + suffix : str + The suffix added to all of the output products. If left blank, the original input file will be + overwritten if save_data is True. + + Returns + ------- + new_data : numpy.ndarray + The wisp-subtracted data. + + wisp_model : numpy.ndarray + The wisp model subtracted from the input file. + + factors : numpy.ndarray + The wisp scale factors tested. + + residuals : numpy.ndarray + The residuals that resulted when applying each of the wisp scale factors tested. The exact + meaning of the residuals depends on the chosen scale_method. + + factor : float + The chosen scale factor applied to the wisp template to create wisp_model. + + Outputs + ------- + {f}{suffix}.fits : fits image + The wisp-subtracted input file. Only generated if save_data is True. + + {f}{suffix}_model.fits : fits image + The wisp model that was subtracted from the input file. Only generated if save_model is True. + + {f}{suffix}_plots.png : png image + A diagnostic plot showing the original data, wisp-subtracted data, final wisp model, and + a plot of the various wisp scale factors tested with their corresponding residuals. + Only generated if plot is True. + """ + + # Get the input data + print('Applying wisp template to {}'.format(f)) + h = fits.open(f) + data = h['SCI'].data + dq = h['DQ'].data + if segmap_data is None: + print('Warning: No segmap data given for {}. Assuming no sources.'.format(f)) + segmap_data = np.zeros(data.shape) + + # Smooth wisp template + if gauss_smooth_wisp: + g = Gaussian2DKernel(x_stddev=gauss_stddev) + wisp_data = convolve(wisp_data, g) + + # Scale the wisp template + if scale_wisp: + # Make a mask of the brightest wisp region + mean, med, stddev = sigma_clipped_stats(wisp_data) + wisp_mask = (wisp_data > med + 2 * stddev).astype(int) + + # Make versions of the original data and wisp model where only good + # pixels in the wisp region are unmasked. + data_masked = np.copy(data) + wisp_data_masked = np.copy(wisp_data) + data_masked[(dq & 1 != 0) | (segmap_data != 0) | (wisp_mask == 0)] = np.nan + wisp_data_masked[(dq & 1 != 0) | (segmap_data != 0) | (wisp_mask == 0)] = np.nan + + # Make a version of the original data where only good pixels outside + # the wisp region are unmasked. + data_masked_ff = np.copy(data) + data_masked_ff[(dq & 1 != 0) | (segmap_data != 0) | (wisp_mask != 0)] = np.nan + + # Correct median-collapsed row/column offsets, representing the 1/f residuals + # and odd-even column residuals and amp offsets, respectively. + med = np.nanmedian(data_masked_ff) + if correct_rows: + collapsed_rows = np.nanmedian(data_masked_ff - med, axis=1) + else: + collapsed_rows = np.zeros(2048) + if correct_cols: + collapsed_cols = np.nanmedian(data_masked_ff - med, axis=0) + else: + collapsed_cols = np.zeros(2048) + correction_image = np.tile(collapsed_cols, (2048, 1)) + np.swapaxes(np.tile(collapsed_rows, (2048, 1)), 0, 1) + data_masked = data_masked - correction_image + data_masked_ff = data_masked_ff - correction_image + med = np.nanmedian(data_masked_ff) + + # Scale wisp template and record residuals + factors = np.arange(factor_min, factor_max, factor_step) + residuals = [] + for factor in factors: + wisp_model = wisp_data_masked * factor + new_data = data_masked - wisp_model + if scale_method == 'mad': + residual = median_absolute_deviation(new_data, ignore_nan=True) + if scale_method == 'median': + residual = abs(med - np.nanmedian(new_data)) + residuals.append(residual) + residuals = np.array(residuals) + + # Choose the wisp template with the lowest noise + if poly_degree != 0: + fn = np.poly1d(np.polyfit(factors, residuals, poly_degree)) # smooth results + factor = factors[np.argmin(fn(factors))] + else: + factor = factors[np.argmin(residuals)] + wisp_model = wisp_data * factor + else: + # Don't scale the wisp data + factors, residuals = np.ones(1), np.zeros(1) # dummy data + factor = 1.0 + wisp_model = wisp_data * factor + + # Only subtract wisp values above the specified threshold + if min_wisp is not None: + wisp_model[wisp_model < min_wisp] = 0 + + # Flag wisp values above the specified thereshold in DQ array + if flag_wisp_thresh is not None: + new_dq = np.copy(dq) + new_dq[(dq & dq_val == 0) & (wisp_model > flag_wisp_thresh)] += dq_val + else: + new_dq = dq + + # Subtract the final wisp template + if sub_wisp: + new_data = data - wisp_model + else: + new_data = data + + # Save the wisp-subtracted data and model + if save_data: + h['SCI'].data = new_data.astype('float32') + h['DQ'].data = new_dq + h.writeto(f.replace('.fits', '{}.fits'.format(suffix)), overwrite=True) + if save_model: + fits.writeto(f.replace('.fits', '{}_model.fits'.format(suffix)), wisp_model, overwrite=True) + h.close() + + # Make diagnostic plots + if plot: + fig, axes = plt.subplots(1, 4, figsize=(60, 10)) + for ax in axes: + ax.tick_params(axis='both', which='major', labelsize=20) + # Plot original image + mean, med, stddev = sigma_clipped_stats(data) + vmin, vmax = med-3*stddev, med+3*stddev + im = axes[0].imshow(data, origin='lower', vmin=vmin, vmax=vmax, cmap='gray') + cbar = fig.colorbar(im, ax=axes[0]) + cbar.ax.tick_params(labelsize=20) + axes[0].set_title('Original', fontsize=35) + # Plot wisp-corrected image + im = axes[1].imshow(new_data, origin='lower', vmin=vmin, vmax=vmax, cmap='gray') + cbar = fig.colorbar(im, ax=axes[1]) + cbar.ax.tick_params(labelsize=20) + axes[1].set_title('Wisp Corrected', fontsize=35) + # Plot wisp model + mean, med, stddev = sigma_clipped_stats(wisp_model) + im = axes[2].imshow(wisp_model, origin='lower', vmin=0, vmax=med+5*stddev, cmap='gray') + cbar = fig.colorbar(im, ax=axes[2]) + cbar.ax.tick_params(labelsize=20) + axes[2].set_title('Wisp Model', fontsize=35) + # Plot scatterplot of wisp scale factors tested vs residuals + if (scale_wisp) & (poly_degree != 0): + axes[3].plot(factors, fn(factors), color='black') + axes[3].scatter(factors, residuals) + axes[3].axvline(factor, color='black') + y_label = 'Residual Median Absolute Deviation' if scale_method == 'mad' else 'Residual Median Difference' + axes[3].set_ylabel(y_label, fontsize=25) + axes[3].set_xlabel('Wisp Factor', fontsize=25) + axes[3].set_title('Chosen Wisp Factor: {:.3f}'.format(factor), fontsize=35) + axes[3].grid(ls='--', alpha=0.5) + fig.savefig(f.replace('.fits', '{}_plots.png'.format(suffix)), dpi=100, bbox_inches='tight') + if not show_plot: + plt.close() + + return new_data, wisp_model, factors, residuals, factor + + +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- + + +def parse_args(): + """ + Parses command line arguments. + + See make_segmap() and subtract_wisp() docstrings for more details on all of the input arguments. + + Returns + ------- + args : object + Contains the input arguments. + """ + + # Make the help strings + files_help = 'The files to subtract the wisp templates from. Wildcards are supported, e.g. ./data/*nrca3*_cal.fits' + nproc_help = 'The number of processes to use during multiprocessing.' + wisp_dir_help = 'The directory containing the wisp templates. The templates are assumed to have the form WISP_{DETECTOR}_{FILTER}_{PUPIL}.fits.' + create_segmap_help = 'Option to make a source segmentation map to help with scaling the wisp template.' + + # Add the potential arguments + parser = argparse.ArgumentParser() + parser.add_argument('--files', dest='files', action='store', nargs='+', type=str, required=False, help=files_help, default='./*_cal.fits') + parser.add_argument('--nproc', dest='nproc', action='store', type=int, required=False, help=nproc_help, default=6) + parser.add_argument('--wisp_dir', dest='wisp_dir', action='store', type=str, required=False, help=wisp_dir_help, default='./') + parser.add_argument('--create_segmap', dest='create_segmap', action=argparse.BooleanOptionalAction, required=False, help=create_segmap_help, default=True) + + # Add arguments for make_segmap() + parser.add_argument('--seg_from_lw', dest='seg_from_lw', action=argparse.BooleanOptionalAction, required=False, default=True) + parser.add_argument('--sigma', dest='sigma', action='store', type=float, required=False, default=0.8) + parser.add_argument('--npixels', dest='npixels', action='store', type=int, required=False, default=10) + parser.add_argument('--dilate_segmap', dest='dilate_segmap', action='store', type=int, required=False, default=5) + parser.add_argument('--save_segmap', dest='save_segmap', action=argparse.BooleanOptionalAction, required=False, default=False) + + # Add arguments for subtract_wisp() + parser.add_argument('--sub_wisp', dest='sub_wisp', action=argparse.BooleanOptionalAction, required=False, default=True) + parser.add_argument('--gauss_smooth_wisp', dest='gauss_smooth_wisp', action=argparse.BooleanOptionalAction, required=False, default=False) + parser.add_argument('--gauss_stddev', dest='gauss_stddev', action='store', type=float, required=False, default=3.0) + parser.add_argument('--scale_wisp', dest='scale_wisp', action=argparse.BooleanOptionalAction, required=False, default=True) + parser.add_argument('--scale_method', dest='scale_method', action='store', type=str, required=False, default='mad') + parser.add_argument('--poly_degree', dest='poly_degree', action='store', type=int, required=False, default=5) + parser.add_argument('--factor_min', dest='factor_min', action='store', type=float, required=False, default=0.0) + parser.add_argument('--factor_max', dest='factor_max', action='store', type=float, required=False, default=2.0) + parser.add_argument('--factor_step', dest='factor_step', action='store', type=float, required=False, default=0.01) + parser.add_argument('--min_wisp', dest='min_wisp', action='store', type=float, required=False, default=None) + parser.add_argument('--flag_wisp_thresh', dest='flag_wisp_thresh', action='store', type=float, required=False, default=None) + parser.add_argument('--dq_val', dest='dq_val', action='store', type=int, required=False, default=1) + parser.add_argument('--correct_rows', dest='correct_rows', action=argparse.BooleanOptionalAction, required=False, default=True) + parser.add_argument('--correct_cols', dest='correct_cols', action=argparse.BooleanOptionalAction, required=False, default=False) + parser.add_argument('--save_data', dest='save_data', action=argparse.BooleanOptionalAction, required=False, default=True) + parser.add_argument('--save_model', dest='save_model', action=argparse.BooleanOptionalAction, required=False, default=True) + parser.add_argument('--plot', dest='plot', action=argparse.BooleanOptionalAction, required=False, default=True) + parser.add_argument('--show_plot', dest='show_plot', action=argparse.BooleanOptionalAction, required=False, default=False) + parser.add_argument('--suffix', dest='suffix', action='store', type=str, required=False, default='_wisp') + + # Get the arguments + args = parser.parse_args() + + return args + + +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- + + +if __name__ == '__main__': + + # Get the command line arguments + args = parse_args() + + # Process the input files + results = process_files(**vars(args)) + print('subtract_wisp.py complete.')