{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Double 2D Sersic example\n",
    "\n",
    "In this example we create a (simplified) synthetic galaxy image consisting of two Sersic components, add some \"realism\" to it (PSF + noise), and then run statmorph in order to recover the parameters of the two components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.visualization import simple_norm\n",
    "from astropy.modeling.models import Sersic2D\n",
    "from astropy.convolution import convolve, Gaussian2DKernel\n",
    "from photutils.segmentation import detect_threshold, detect_sources\n",
    "import time\n",
    "import statmorph\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up\n",
    "\n",
    "**Creating the model galaxy image**\n",
    "\n",
    "We assume that the image size is 240x240 pixels and that the \"true\" light distribution corresponds to a *double* 2D Sersic model with the following parameters (note that the two components share the same center, by construction):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ny, nx = 240, 240\n",
    "y, x = np.mgrid[0:ny, 0:nx]\n",
    "doublesersic_model = statmorph.DoubleSersic2D(\n",
    "    x_0=120.5, y_0=96.5,\n",
    "    amplitude_1=1, r_eff_1=10, n_1=5.0, ellip_1=0.6, theta_1=2.0,\n",
    "    amplitude_2=2, r_eff_2=20, n_2=1.0, ellip_2=0.4, theta_2=0.5)\n",
    "image = doublesersic_model(x, y)\n",
    "\n",
    "# Visualize \"idealized\" image\n",
    "plt.imshow(image, cmap='gray', origin='lower',\n",
    "           norm=simple_norm(image, stretch='log', log_a=10000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Applying realism**\n",
    "\n",
    "We now apply some \"realism\" (PSF + noise) to the idealized image (see the [tutorial](./tutorial.html) for more details):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convolve with PSF\n",
    "kernel = Gaussian2DKernel(2.0)\n",
    "kernel.normalize()  # make sure kernel adds up to 1\n",
    "psf = kernel.array  # we only need the numpy array\n",
    "image = convolve(image, psf)\n",
    "\n",
    "# Apply shot noise\n",
    "np.random.seed(3)\n",
    "gain = 1e5\n",
    "image = np.random.poisson(image * gain) / gain\n",
    "\n",
    "# Apply background noise\n",
    "sky_sigma = 0.01\n",
    "image += sky_sigma * np.random.standard_normal(size=(ny, nx))\n",
    "\n",
    "# Visualize \"realistic\" image\n",
    "plt.imshow(image, cmap='gray', origin='lower',\n",
    "           norm=simple_norm(image, stretch='log', log_a=10000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Creating a segmentation map**\n",
    "\n",
    "We also need to create a segmentation image (see the [tutorial](./tutorial.html) for more details):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = detect_threshold(image, 1.5)\n",
    "npixels = 5  # minimum number of connected pixels\n",
    "convolved_image = convolve(image, psf)\n",
    "segmap = detect_sources(convolved_image, threshold, npixels)\n",
    "plt.imshow(segmap, origin='lower', cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Running statmorph\n",
    "\n",
    "We now have all the input necessary to run statmorph. However, unlike in the [tutorial](./tutorial.html), this time we include the option ``include_doublesersic = True``, which is necessary in order to carry out the double Sersic fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = time.time()\n",
    "source_morphs = statmorph.source_morphology(\n",
    "    image, segmap, gain=gain, psf=psf, include_doublesersic=True)\n",
    "print('Time: %g s.' % (time.time() - start))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examining the output\n",
    "\n",
    "We focus on the first (and only) source labeled in the segmap:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "morph = source_morphs[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We print some of the morphological properties just calculated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('BASIC MEASUREMENTS (NON-PARAMETRIC)')\n",
    "print('xc_centroid =', morph.xc_centroid)\n",
    "print('yc_centroid =', morph.yc_centroid)\n",
    "print('ellipticity_centroid =', morph.ellipticity_centroid)\n",
    "print('elongation_centroid =', morph.elongation_centroid)\n",
    "print('orientation_centroid =', morph.orientation_centroid)\n",
    "print('xc_asymmetry =', morph.xc_asymmetry)\n",
    "print('yc_asymmetry =', morph.yc_asymmetry)\n",
    "print('ellipticity_asymmetry =', morph.ellipticity_asymmetry)\n",
    "print('elongation_asymmetry =', morph.elongation_asymmetry)\n",
    "print('orientation_asymmetry =', morph.orientation_asymmetry)\n",
    "print('rpetro_circ =', morph.rpetro_circ)\n",
    "print('rpetro_ellip =', morph.rpetro_ellip)\n",
    "print('rhalf_circ =', morph.rhalf_circ)\n",
    "print('rhalf_ellip =', morph.rhalf_ellip)\n",
    "print('r20 =', morph.r20)\n",
    "print('r80 =', morph.r80)\n",
    "print('Gini =', morph.gini)\n",
    "print('M20 =', morph.m20)\n",
    "print('F(G, M20) =', morph.gini_m20_bulge)\n",
    "print('S(G, M20) =', morph.gini_m20_merger)\n",
    "print('sn_per_pixel =', morph.sn_per_pixel)\n",
    "print('C =', morph.concentration)\n",
    "print('A =', morph.asymmetry)\n",
    "print('S =', morph.smoothness)\n",
    "print()\n",
    "print('SINGLE SERSIC')\n",
    "print('sersic_amplitude =', morph.sersic_amplitude)\n",
    "print('sersic_rhalf =', morph.sersic_rhalf)\n",
    "print('sersic_n =', morph.sersic_n)\n",
    "print('sersic_xc =', morph.sersic_xc)\n",
    "print('sersic_yc =', morph.sersic_yc)\n",
    "print('sersic_ellip =', morph.sersic_ellip)\n",
    "print('sersic_theta =', morph.sersic_theta)\n",
    "print('sersic_chi2_dof =', morph.sersic_chi2_dof)\n",
    "print('sersic_aic =', morph.sersic_aic)\n",
    "print('sersic_bic =', morph.sersic_bic)\n",
    "print()\n",
    "print('DOUBLE SERSIC')\n",
    "print('doublesersic_xc =', morph.doublesersic_xc)\n",
    "print('doublesersic_yc =', morph.doublesersic_yc)\n",
    "print('doublesersic_amplitude1 =', morph.doublesersic_amplitude1)\n",
    "print('doublesersic_rhalf1 =', morph.doublesersic_rhalf1)\n",
    "print('doublesersic_n1 =', morph.doublesersic_n1)\n",
    "print('doublesersic_ellip1 =', morph.doublesersic_ellip1)\n",
    "print('doublesersic_theta1 =', morph.doublesersic_theta1)\n",
    "print('doublesersic_amplitude2 =', morph.doublesersic_amplitude2)\n",
    "print('doublesersic_rhalf2 =', morph.doublesersic_rhalf2)\n",
    "print('doublesersic_n2 =', morph.doublesersic_n2)\n",
    "print('doublesersic_ellip2 =', morph.doublesersic_ellip2)\n",
    "print('doublesersic_theta2 =', morph.doublesersic_theta2)\n",
    "print('doublesersic_chi2_dof =', morph.doublesersic_chi2_dof)\n",
    "print('doublesersic_aic =', morph.doublesersic_aic)\n",
    "print('doublesersic_bic =', morph.doublesersic_bic)\n",
    "print()\n",
    "print('OTHER')\n",
    "print('sky_mean =', morph.sky_mean)\n",
    "print('sky_median =', morph.sky_median)\n",
    "print('sky_sigma =', morph.sky_sigma)\n",
    "print('flag =', morph.flag)\n",
    "print('flag_sersic =', morph.flag_sersic)\n",
    "print('flag_doublesersic =', morph.flag_doublesersic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, the fitted double Sersic model is consistent with the \"true\" light distribution that we originally created (n1 = 5, n2 = 1, etc.) and the reduced chi-squared statistic (doublesersic_chi2_dof) is close to 1, indicating a good fit without overfitting. On the other hand, the *single* Sersic fit has a reduced chi-squared statistic much larger than 1, indicating a poor fit (as expected).\n",
    "\n",
    "We also calculate the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) for the two models, which again favor the *double* Sersic model as the statistically preferred one, since it returns much lower AIC and BIC values.\n",
    "\n",
    "Also note that statmorph now returns an additional quality flag:\n",
    "\n",
    "- ``flag_doublesersic`` : indicates the quality of the double Sersic fit. Like ``flag`` and ``flag_sersic``, it can take the following values: 0 (good), 1 (suspect), 2 (bad), and 4 (catastrophic)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Visualizing the individual components**\n",
    "\n",
    "For some applications (e.g. bulge/disk decompositions) it might be useful to analyze the two fitted components separately, as we do below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ny, nx = image.shape\n",
    "y, x = np.mgrid[0:ny, 0:nx]\n",
    "sersic1 = Sersic2D(morph.doublesersic_amplitude1,\n",
    "                   morph.doublesersic_rhalf1,\n",
    "                   morph.doublesersic_n1,\n",
    "                   morph.doublesersic_xc,\n",
    "                   morph.doublesersic_yc,\n",
    "                   morph.doublesersic_ellip1,\n",
    "                   morph.doublesersic_theta1)\n",
    "sersic2 = Sersic2D(morph.doublesersic_amplitude2,\n",
    "                   morph.doublesersic_rhalf2,\n",
    "                   morph.doublesersic_n2,\n",
    "                   morph.doublesersic_xc,\n",
    "                   morph.doublesersic_yc,\n",
    "                   morph.doublesersic_ellip2,\n",
    "                   morph.doublesersic_theta2)\n",
    "image1 = sersic1(x, y)\n",
    "image2 = sersic2(x, y)\n",
    "image_total = image1 + image2\n",
    "\n",
    "fig = plt.figure(figsize=(15,5))\n",
    "ax = fig.add_subplot(131)\n",
    "ax.imshow(image1, cmap='gray', origin='lower',\n",
    "           norm=simple_norm(image_total, stretch='log', log_a=10000))\n",
    "ax.set_title('First component')\n",
    "ax.text(0.04, 0.93, 'n1 = %.4f' % (morph.doublesersic_n1,),\n",
    "        bbox=dict(facecolor='white'), transform=ax.transAxes)\n",
    "ax = fig.add_subplot(132)\n",
    "ax.imshow(image2, cmap='gray', origin='lower',\n",
    "           norm=simple_norm(image_total, stretch='log', log_a=10000))\n",
    "ax.set_title('Second component')\n",
    "ax.text(0.04, 0.93, 'n2 = %.4f' % (morph.doublesersic_n2,),\n",
    "        bbox=dict(facecolor='white'), transform=ax.transAxes)\n",
    "ax = fig.add_subplot(133)\n",
    "ax.imshow(image_total, cmap='gray', origin='lower',\n",
    "           norm=simple_norm(image_total, stretch='log', log_a=10000))\n",
    "ax.set_title('Composite model')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the two Sersic components shown above are *not* convolved with the PSF, since they are meant to recover the \"true\" light distributions of the two components of the galaxy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Examining the single Sersic fit**\n",
    "\n",
    "For illustration puposes, below we compare the original (realistic) image to the *single* Sersic fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bg_noise = sky_sigma * np.random.standard_normal(size=(ny, nx))\n",
    "model = statmorph.ConvolvedSersic2D(\n",
    "    morph.sersic_amplitude,\n",
    "    morph.sersic_rhalf,\n",
    "    morph.sersic_n,\n",
    "    morph.sersic_xc,\n",
    "    morph.sersic_yc,\n",
    "    morph.sersic_ellip,\n",
    "    morph.sersic_theta)\n",
    "model.set_psf(psf)  # must set PSF by hand\n",
    "image_model = model(x, y)\n",
    "\n",
    "fig = plt.figure(figsize=(15,5))\n",
    "ax = fig.add_subplot(131)\n",
    "ax.imshow(image, cmap='gray', origin='lower',\n",
    "           norm=simple_norm(image, stretch='log', log_a=10000))\n",
    "ax.set_title('Original image')\n",
    "ax = fig.add_subplot(132)\n",
    "ax.imshow(image_model + bg_noise, cmap='gray', origin='lower',\n",
    "           norm=simple_norm(image, stretch='log', log_a=10000))\n",
    "ax.set_title('Single Sersic fit')\n",
    "ax = fig.add_subplot(133)\n",
    "residual = image - image_model\n",
    "ax.imshow(residual, cmap='gray', origin='lower',\n",
    "           norm=simple_norm(residual, stretch='linear'))\n",
    "ax.set_title('Single Sersic residual')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Examining the double Sersic fit**\n",
    "\n",
    "Similarly, below we compare the original (realistic) image to the *double* Sersic fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = statmorph.ConvolvedDoubleSersic2D(\n",
    "    morph.doublesersic_xc,\n",
    "    morph.doublesersic_yc,\n",
    "    morph.doublesersic_amplitude1,\n",
    "    morph.doublesersic_rhalf1,\n",
    "    morph.doublesersic_n1,\n",
    "    morph.doublesersic_ellip1,\n",
    "    morph.doublesersic_theta1,\n",
    "    morph.doublesersic_amplitude2,\n",
    "    morph.doublesersic_rhalf2,\n",
    "    morph.doublesersic_n2,\n",
    "    morph.doublesersic_ellip2,\n",
    "    morph.doublesersic_theta2)\n",
    "model.set_psf(psf)  # must set PSF by hand\n",
    "image_model = model(x, y)\n",
    "\n",
    "fig = plt.figure(figsize=(15,5))\n",
    "ax = fig.add_subplot(131)\n",
    "ax.imshow(image, cmap='gray', origin='lower',\n",
    "           norm=simple_norm(image, stretch='log', log_a=10000))\n",
    "ax.set_title('Original image')\n",
    "ax = fig.add_subplot(132)\n",
    "ax.imshow(image_model + bg_noise, cmap='gray', origin='lower',\n",
    "           norm=simple_norm(image, stretch='log', log_a=10000))\n",
    "ax.set_title('Double Sersic fit')\n",
    "ax = fig.add_subplot(133)\n",
    "residual = image - image_model\n",
    "ax.imshow(residual, cmap='gray', origin='lower',\n",
    "           norm=simple_norm(residual, stretch='linear'))\n",
    "ax.set_title('Double Sersic residual')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig('doublesersic.png', dpi=150)\n",
    "plt.close(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Concluding remarks\n",
    "\n",
    "The fact that statmorph uses Astropy's modeling utility behind the scenes provides a great deal of flexibility. For example, if one is interested in fitting a de Vaucouleurs + exponential model (these components are, of course, special cases of the Sersic model with `n = 4` and `n = 1`, respectively), one simply has to add the following option when calling statmorph:\n",
    "\n",
    "    doublesersic_model_args = {\n",
    "        'n_1': 4, 'n_2': 1, 'fixed': {'n_1': True, 'n_2': True}}\n",
    "\n",
    "Furthermore, in some applications it might make sense to \"tie\" the ellipticity and position angle of the two Sersic components. This can also be accomplished using ``doublesersic_model_args`` in combination with the ``tied`` property of Astropy parameters, although the syntax is slightly more involved (more details [here](https://docs.astropy.org/en/stable/modeling/parameters.html)). Alternatively, statmorph provides the following option for this purpose, which achieves the same effect:\n",
    "\n",
    "    doublesersic_tied_ellip = 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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}