{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "09eb9c87",
   "metadata": {},
   "source": [
    "# Correction Grids\n",
    "\n",
    "Correction grids can be thought of as lookup tables that provide corrected Petrosian derived values. \n",
    "These grids are typically generated through simulations or theoretical models to offer \n",
    "the best correction for Petrosian profiles. This section discusses the concept of correction grids: how they're generated,\n",
    "when and why they should be used, and how they enhance accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb595677",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "# Hidden cell\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d48236c6",
   "metadata": {},
   "source": [
    "## Necessity and Implications\n",
    "\n",
    "**Are Correction Grids Always Needed?**\n",
    "\n",
    "The short answer is no, but they are very useful in some cases. \n",
    "Petrosian corrections offer more precise measurements of galaxies. \n",
    "The necessity hinges on your data quality, research goals, and desired precision. \n",
    "Here are key considerations on why the default settings are usually good enough:\n",
    "\n",
    "1. **Standard Practices**: Numerous studies, such as those from the Sloan Digital Sky Survey (SDSS), \n",
    "   rely on default settings. The parameters $\\eta = 0.2$ and $\\epsilon = 2$ are widely accepted because\n",
    "   they typically produce reliable results across many scenarios.\n",
    "\n",
    "2. **Noise Constraints**: If the fainter parts of galaxies in your images are submerged below the noise level, \n",
    "   detailed corrections might not significantly impact the outcome since that portion of the galaxy is not meaningfully\n",
    "   measurable. Here, the additional steps and computational demands of correction grids might be unwarranted.\n",
    "\n",
    "3. **Practical Implications**: While correction grids can produce improved results, \n",
    "   they require optimal background subtraction. Not all datasets might meet this standard.\n",
    " \n",
    "Correction grids promise enhanced precision but aren't always essential. \n",
    "For a majority of cases, default settings suffice.\n",
    "\n",
    "**When are Correction Grids Needed?**\n",
    "\n",
    "While the default settings serve most general purposes, there are specific scenarios where \n",
    "correction grids are beneficial:\n",
    "\n",
    "1. **Galaxy Size Estimation**: If your research aims to measure the size of a galaxy which suffers due to \n",
    "   poor signal-to-noise, you can estimate its total flux radius based on the corrected epsilon. Remember, \n",
    "   this isn't about recovering the total flux which includes regions submerged under noise, but \n",
    "   about the estimating the full extent of the galaxy based on the measurable components of the light profile. \n",
    "\n",
    "2. **Estimating Sérsic Index**: Correction grids can be used to estimate the Sérsic index of single component galaxies.\n",
    "\n",
    "3. **High-Precision Studies**: For projects that require a granular study of galaxy structures, \n",
    "   or when analyzing galaxies with distinct features that might get overlooked with standard parameters, \n",
    "   correction grids provide the extra precision that may be needed.\n",
    "\n",
    "4. **Superior Data Quality**: If you have high-quality images with excellent background subtraction \n",
    "   and reasonable signal-to-noise ratios, the total flux can be estimated to higher accuracy. That is,\n",
    "   high quality data allows you to measure the flux up to the corrected total flux radius. \n",
    "\n",
    "In essence, when the goal shifts from general measurements to intricate, high-resolution studies of galactic structures and sizes, correction grids become a useful tool.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b119cb69",
   "metadata": {},
   "source": [
    "## Generating a Correction Grid\n",
    "\n",
    "Correction grids are essentially lookup tables that map certain observational properties of a galaxy \n",
    "(like Petrosian radius, uncorrected half light radius, and concentration indices) to more intrinsic \n",
    "properties (like Sérsic index and the corrected epsilon value). In this section, we'll walk through the steps of creating a correction grid. \n",
    "\n",
    "To start with `PetroFit`, simply import it as follows:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc5cca42",
   "metadata": {},
   "outputs": [],
   "source": [
    "import petrofit as pf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef5216be",
   "metadata": {},
   "source": [
    "### Load PSF (Optional)\n",
    "\n",
    "We load an HST `F105W` PSF and normalize it to sum to 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "573a17f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.rcParams['figure.figsize'] = [6, 6]\n",
    "plt.rcParams['image.origin'] = 'lower'\n",
    "plt.rcParams['font.size'] = 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "potential-princeton",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "\n",
    "# Load PSF image (2D array)\n",
    "PSF = fits.getdata('data/f105w_psf.fits.gz')\n",
    "\n",
    "# Normalize PSF \n",
    "PSF = PSF / PSF.sum()\n",
    "\n",
    "# Note that the PSF shape is odd on all sides\n",
    "print(\"PSF Shape = {}\".format(PSF.shape))\n",
    "\n",
    "# Plot PSF and use vmax and vmin to show difraction spikes\n",
    "plt.imshow(PSF, vmin=0, vmax=5e-4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c20bce2",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "### Define Sampling Points \n",
    "\n",
    "The correction grid is generated along with a set of half-light radii and Sersic indices. The generator loops through the Sersic indices for each half-light radius in the radius list. For this documentation, we define a small set of values, with the Sersic indices of a Gaussian and de Vaucouleurs’ profile.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcec33cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "\n",
    "r_eff_list = np.array([7, 10, 15])\n",
    "n_list = np.array([0.5, 1., 4.])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f5cf7ae",
   "metadata": {},
   "source": [
    "### Grid Simulation\n",
    "\n",
    "We do a simple API call to generate the correction grid. We provide a PSF and oversampling rule as well. Oversampling becomes important for small half-light radii since the model needs to be sampled well in the center. \n",
    "\n",
    "The `generate_petrosian_sersic_correction` grid follows the following steps to generate the correction grid:\n",
    "\n",
    "- Computes the total flux (`total_flux`) of an ideal Sersic profile with the sampling points using the ` petrofit.modeling.models.sersic_enclosed` function and setting the radius to `np.inf`.\n",
    "- Computes the radius equal to `total_flux * 0.99`.\n",
    "- Makes a PSF convolved Sersic Model image and measures the photometry. \n",
    "- Measures the uncorrected Petrosian radius. \n",
    "- Computes the correct `epsilon` value as ` corrected_epsilon = r(total_flux) / r_petrosian`.\n",
    "- It also saves other values, such as the uncorrected `C2080` which can be used to map to the Sersic index.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dedicated-helen",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_file_name='temp/example_correction_gid.ecsv'\n",
    "\n",
    "petrosian_grid = pf.generate_petrosian_sersic_correction(\n",
    "    output_file_name=output_file_name, # Output file name\n",
    "    psf=PSF, # PSF (optional)\n",
    "    r_eff_list=r_eff_list, # List of r_e to sample\n",
    "    n_list=n_list, # List of n to sample\n",
    "    oversample=4, # Oversample factor or region, see fitting docs\n",
    "    out_format='ascii.ecsv', # Output format, should match file name\n",
    "    overwrite=True, # Overwrite output\n",
    "    ipython_widget=True, # Progress bar\n",
    "    n_cpu=None, # int value for number of mp threads\n",
    "    plot=False, # Plot each step, can not be True if n_cpu > 1\n",
    "\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd6daece",
   "metadata": {},
   "outputs": [],
   "source": [
    "petrosian_grid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85317016",
   "metadata": {},
   "source": [
    "## PetrosianCorrection\n",
    "\n",
    "\n",
    "The `PetrosianCorrection` class is designed to compute corrections for Petrosian based on \n",
    "default Petrosian measurements. \n",
    "\n",
    "### Loading and Navigating\n",
    "\n",
    "1. **Grid Input**: To initialize the `PetrosianCorrection` class, you'll need a correction grid. \n",
    "   This grid is typically an Astropy Table that has been generated with the function \n",
    "   `petrofit.correction.generate_petrosian_sersic_correction`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9177ef8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input generated grid table\n",
    "pc = pf.PetrosianCorrection(petrosian_grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb64600f",
   "metadata": {},
   "source": [
    "2. **Reading and Writing Grids**: \n",
    "    - Use the `PetrosianCorrection.read(file_path)` class method to read a correction grid from a file.\n",
    "    - Similarly, to save your correction grid to a file, utilize the `write(grid_file, file_format=None)` method.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53e324c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read F105W grid:\n",
    "pc = pf.PetrosianCorrection.read(output_file_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0899985b",
   "metadata": {},
   "source": [
    "\n",
    "3. **Working with Grid Data**:\n",
    "    - The `grid_keys` property will provide you with the column names of your grid.\n",
    "    - Use `unique_grid_values(key)` to obtain unique values from a specific column.\n",
    "    - The `filter_grid(key, value)` method helps filter the grid based on a specific key and value.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8972898",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read ideal Sersic grid:\n",
    "pc = pf.PetrosianCorrection.read('data/no_psf_corr.ecsv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4332cd16",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pc.grid_keys)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7d3bf50",
   "metadata": {},
   "outputs": [],
   "source": [
    "# List all Sersic indices available \n",
    "pc.unique_grid_values('n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15b85c14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Provide rows for a grid key's value\n",
    "# For example all rows with n=1\n",
    "pc.filter_grid('n', 1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a0dfb85",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting relations \n",
    "pc_c_c2080_list = np.sort(pc.unique_grid_values('c_c2080')) # Get corrected c_c2080 \n",
    "approx_n_values = pf.PetroApprox.c2080_to_n(pc_c_c2080_list) # Get n from PetroApprox for comparison\n",
    "\n",
    "# Plot\n",
    "plt.scatter( pc.grid['c_c2080'], pc.grid['n'], color='black')\n",
    "plt.plot(pc_c_c2080_list, approx_n_values, color='red')\n",
    "\n",
    "plt.xlabel('Corrected $C_{2080}$')\n",
    "plt.ylabel('n')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0864e222",
   "metadata": {},
   "source": [
    "### Correcting a Petrosian Profile\n",
    "\n",
    "In the [Photometry Chapter](./photometry.ipynb#Photometry) we constructed a curve of growth for the football shaped galaxy:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3fbce43",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy.table import Table\n",
    "phot_table = Table.read('data/abell_2744_galaxy_f105w_photometry.ecsv') # Read table \n",
    "\n",
    "# Load data\n",
    "r_list = np.array(phot_table['r_list'])\n",
    "flux_arr = np.array(phot_table['flux_arr'])\n",
    "area_arr = np.array(phot_table['area_arr'])\n",
    "error_arr = np.array(phot_table['error_arr'])\n",
    "\n",
    "\n",
    "# Make Petrosian profile\n",
    "p = pf.Petrosian(r_list, area_arr, flux_arr, flux_err=error_arr)\n",
    "\n",
    "p.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78f75f7f",
   "metadata": {},
   "source": [
    "1. **Correct Profile**:\n",
    "    - Simply use the `correct(p)` function to correct the profile. \n",
    "    - Use `plot_correction` function to see the closest value.\n",
    "    - The parameters used for correction are derived directly from the grid. Specifically:\n",
    "        - `'p02'`: Represents the Petrosian radius $P_{02} = R(\\eta_{0.2})$.\n",
    "        - `'u_r_50'`: Stands for the uncorrected half-light radius.\n",
    "        - `'u_c2080'`: Corresponds to the uncorrected concentration index.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7ab8441",
   "metadata": {},
   "outputs": [],
   "source": [
    "pc = pf.PetrosianCorrection.read('data/f105w_psf_corr.ecsv')\n",
    "\n",
    "# Plot grid and uncorrected profile on grid (p)\n",
    "pc.plot_correction(p)\n",
    "plt.show()\n",
    "\n",
    "# Pass uncorrected p to the correct function\n",
    "p_corrected = pc.correct(p)\n",
    "\n",
    "# Plot corrected profile\n",
    "p_corrected.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e775a2d7",
   "metadata": {},
   "source": [
    "2. **Estimating Parameters**:\n",
    "    - The method `estimate_n(p)` will give an estimated Sérsic index `n` based on the half-light radius and `c2080` computed using the default epsilon value.\n",
    "    - The method `estimate_epsilon(p)` will provide a corrected epsilon value given the half-light radius and `c2080` computed with the default epsilon.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67ab440d",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = pc.estimate_n(p)\n",
    "n "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f4538ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_epsilon = pc.estimate_epsilon(p)\n",
    "corr_epsilon"
   ]
  }
 ],
 "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}