{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9e4ab6ac",
   "metadata": {},
   "source": [
    "# Results Analysis\n",
    "\n",
    "We have successfully run a sampling on some source with a given photometric system using the terminal commands (i.e. ``galapy-fit``).\n",
    "We now want to analyse the results obtained. To this end, we have to first load the results object that has been stored at the end of the sampling run.\n",
    "\n",
    "This can be done with the dedicated function ``galapy.sampling.Results.load_results``.\n",
    "First of all we load the function from the ``galapy.sampling`` sub-package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "886aa766",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.sampling.Results import load_results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3192fa04",
   "metadata": {},
   "source": [
    "We can then simply call the function with the path to the output file as an argument.\n",
    "\n",
    "\n",
    "> Since we have selected ``store_lightweigth = True`` in the parameter file for this run,  \n",
    "> loading the results object will require some tens of seconds.\n",
    "> The ``load_results`` function automathically understands what type of file is being passed \n",
    "> (from the extension and from internal meta-data)\n",
    "> if not specified differently in the optional arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41e8ca6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = load_results('data/sampling_4D+noise_dynesty_results_light.galapy.hdf5')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f6ee0b9",
   "metadata": {},
   "source": [
    "``res`` is an instance of type ``Results``, these objects contain:\n",
    "\n",
    "* attributes and functions to re-construct the characteristics of the run, such as\n",
    "    - the model\n",
    "    - the observation\n",
    "    - the handler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4d7a92d",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = res.get_model()\n",
    "observation = res.get_observation()\n",
    "handler = res.get_handler()\n",
    "print('The run has ', res.Ndof, 'degrees of freedom')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2fcada1",
   "metadata": {},
   "source": [
    "* tables of relevant physical quantities pre-computed at each sampled position in the parameter space, \n",
    "  a list can be shown by calling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89b5a5de",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.get_stored_quantities()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1241b023",
   "metadata": {},
   "source": [
    "* a set of useful functions to analyse the results, such as ``get_mean``, ``get_std``, ``get_bestfit``, \n",
    "  all of which take, as first argument a string with one of the quantities stored (i.e. ``res.get_mean('Mstar')`` \n",
    "  will return the **weighted mean** of all the sampled positions).\n",
    "  Weighting is done automathically by these functions.\n",
    "  Some other functions, such as ``res.get_chi2`` and ``res.get_residuals``, \n",
    "  compute quantities related to some chosen statistics \n",
    "  (i.e. ``res.get_chi2('bestfit')`` will return the reduced $\\chi^2$ computed at the \n",
    "  best-fitting values of parameters).\n",
    "  Note that all these functions also account for the noise treatment chosen when running the sampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f3fa343",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Best-fitting stellar mass: {res.get_bestfit('Mstar'):.2e} Msol\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55471573",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Median stellar mass: {res.get_median('Mstar'):.2e} Msol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8421934",
   "metadata": {},
   "source": [
    "> **Tip:** The ``handler`` object contains information on the parameterisation chosen and provides a \n",
    "> function to easily return the argument that has to be passed to the function ``model.set_parameters`` \n",
    "> (i.e. ``handler.return_nested()``).\n",
    "> For instance, it contains the priors limits, the list of free-parameters (by keyword) and which of them has\n",
    "> been sampled in logarithmic space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7bc36a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "handler.par_prior, handler.par_free, handler.par_log"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f32e2ee",
   "metadata": {},
   "source": [
    "## Plots\n",
    "\n",
    "The ``galapy.analysis`` sub-package contains sub-modules to ease the analysis and plotting of a sampling run results. To have access to plotting functions, we can import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e95cbefd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import galapy.analysis.plot as gplot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d31819e",
   "metadata": {},
   "source": [
    "Let us also import the formatted version of ``matplotlib.pyplot`` present inside the module, for consistency: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae6db6ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.analysis.plot import plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad8c3fde",
   "metadata": {},
   "source": [
    "The function ``gplot.sed_flux_res`` allows to plot the sampled model and compare it to the input dataset, by passing the results object as input argument (it also provides some further optional arguments to customise the plot, i.e. axes limits, colours, legend keyword arguments). The returned value is an object of type ``matplotlib.axes.Axes``, that can be accessed by the user to further customize the plot (or overplot further lines and artists to it) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "feaf4c01",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = gplot.sed_flux_res( \n",
    "    res, plot_components=True, plot_observation=True, plot_contours=True,\n",
    "    ax_kwargs = {\n",
    "       'xlim':(1.e+3, 2.e+8),\n",
    "       'ylim':(2.e-6,1.e+3),\n",
    "    },\n",
    "    legend_kwargs = {\n",
    "        'l1': {'loc': 'upper left', 'fontsize':12},\n",
    "        'l2': {'loc': 'upper right', 'ncol': 2, 'fontsize': 12}\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eef15df1",
   "metadata": {},
   "source": [
    "The function ``gplot.sed_residuals_res``, equivalently and consistently, plots the **standardised residuals** (i.e. the residuals in units of the error in the measurement):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ebb03b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can also build matplotlib figure and axes before calling the function\n",
    "# e.g. if we want to overwrite the default figure size:\n",
    "fig, ax = plt.subplots(1,1,figsize=(6,2)) \n",
    "\n",
    "# it is then sufficient to pass the generated axes to the galapy function:\n",
    "_ = gplot.sed_residuals_res( \n",
    "    res, frame='obs', plot_contours=True, plot_chi2 = True,\n",
    "    ax = ax, # <------ here!\n",
    "    ax_kwargs={\n",
    "        'xlim':(1.e+3, 2.e+8),\n",
    "        'ylim':(-3., +3.)\n",
    "    },\n",
    "    text_kwargs={'loc':'lower right'}\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "336437de",
   "metadata": {},
   "source": [
    "The library leaves users freedom to customise the plots as they prefer, by combining the different functions we can, e.g., obtain a two panel plot of the information above: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b73490c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2,1,figsize=(6,5),tight_layout=True, \n",
    "                                      sharex = True, \n",
    "                         gridspec_kw={'height_ratios':(4.5,1.5), 'hspace':0.0})\n",
    "\n",
    "###################################################\n",
    "# Plot the SED:\n",
    "\n",
    "_ = gplot.sed_flux_res( \n",
    "    res, plot_components=True, plot_observation=True, plot_contours=True,\n",
    "    ax = axes[0], # passing the first Axis to the ax argument\n",
    "    ax_kwargs = {\n",
    "       'xlim':(1.e+3, 2.e+8),\n",
    "       'ylim':(2.e-6,1.e+4),\n",
    "    },\n",
    "    legend_kwargs = {\n",
    "        'l1': {'loc': 'upper left', 'fontsize':12},\n",
    "        'l2': {'loc': 'upper right', 'ncol': 2, 'fontsize': 12}\n",
    "    }\n",
    ")\n",
    "\n",
    "###################################################\n",
    "# Plot the residuals\n",
    "\n",
    "_ = gplot.sed_residuals_res( \n",
    "    res, frame='obs', plot_contours=True, plot_chi2 = True,\n",
    "    ax = axes[1], # passing the second Axis to the ax argument\n",
    "    ax_kwargs={\n",
    "        'xlim':(1.e+3, 2.e+8),\n",
    "        'ylim':(-3., +3.)\n",
    "    },\n",
    "    text_kwargs={'loc':'lower right'}\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aeacf3c7",
   "metadata": {},
   "source": [
    "### Plot the derived attenuation curve\n",
    "\n",
    "With GalaPy is possible to compute the derived attenuation curve for a given model among those sampled by the algorithm.\n",
    "\n",
    "Let us first extract the element in the wavelength grid that we will assume to approximate emission in the visible band (so that to normalize the derived attenuation curve)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09d0c1f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from galapy.internal.utils import find_nearest\n",
    "ll = model.wl()\n",
    "w5500 = find_nearest(ll, 5500) # index in the wavelength grid corresponding to 5500 Angstrom"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b55eb1c7",
   "metadata": {},
   "source": [
    "We will plot the average attenuation for the best-fit model, to that purpose we set the parameters to the position in the parameter space that maximazes the likelihood:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db9f0b76",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.set_parameters(**res.get_bestfit('params'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18f0429e",
   "metadata": {},
   "source": [
    "Before actually computing the attenuation for this model we have to compute the emission for this model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbd7d39d",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = model.get_emission()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75901e10",
   "metadata": {},
   "source": [
    "We can now compute the average total attenuation and normalize it to the its value in the visible band:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd0d32cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "Abest = model.get_avgAtt()\n",
    "Abest /= Abest[w5500]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb98357f",
   "metadata": {},
   "source": [
    "Before plotting let's also define a function that computes Calzetti-like ([Calzetti et al., 2020](https://iopscience.iop.org/article/10.1086/308692)) attenuation curves for reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25f3baf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "def att_calzetti ( ll, Rv ) :\n",
    "    \"\"\"from Calzetti et al., 2000\"\"\"\n",
    "    # convert angstrom to micron:\n",
    "    ll = numpy.array(ll) * 1.e-4 \n",
    "    kp = 2.659 * ( - 2.156 + 1.509 / ll - 0.198 / ll**2. + 0.011 / ll**3. ) + Rv\n",
    "    wl = numpy.where( ll >= 0.63 )\n",
    "    kp[wl] = 2.659 * ( - 1.857 + 1.04 / ll[wl] ) + Rv\n",
    "    return kp / Rv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16a4e424",
   "metadata": {},
   "source": [
    "With a bit of cosmetics, we can plot (in black in the plot below) the derived average total attenuation of the best-fitting galaxy model: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d1a99a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot set-up\n",
    "fig, ax = plt.subplots(1,1,figsize=(6,3),constrained_layout=True)\n",
    "ax.set(\n",
    "    xscale='log', yscale='log',\n",
    "    xlim=(9.12e+2,5.e+4), ylim=(2.e-1,1.e+1), \n",
    "    xlabel='$\\\\lambda_\\\\mathrm{rest}\\\\ [\\\\AA]$',\n",
    "    ylabel='$A(\\\\lambda)/A_V$',\n",
    ")\n",
    "\n",
    "#########################\n",
    "# Plot Calzetti reference\n",
    "\n",
    "# Get a sequence of colors for values of the Rv parameter\n",
    "Rvs = [ 2.,4.,6.,8.,10.,20. ]\n",
    "cmaplist = plt.cm.viridis(numpy.linspace(0.1,0.9,len(Rvs)))\n",
    "\n",
    "# Plot for each Rv value\n",
    "_ = [ ax.plot(ll, att_calzetti(ll, Rv), \n",
    "              color=clr, ls=':', lw=2.5, label=f'$R_v = ${Rv:.0f}') \n",
    "      for Rv, clr in zip(Rvs,cmaplist) ]\n",
    "\n",
    "##########################\n",
    "# Plot derived Attenuation\n",
    "\n",
    "ax.plot( ll, Abest, color = 'k', lw=2.5, label='derived' )\n",
    "\n",
    "# legend\n",
    "_ = ax.legend(ncol=2, fontsize=11)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "406f8d11",
   "metadata": {},
   "source": [
    "### Posteriors\n",
    "\n",
    "Another useful diagnostic plot one might want to check, is the triangular 2D and 1D marginalisation of the parameters posteriors (i.e. the contour plot, corner plot or triangle plot).\n",
    "\n",
    "This can be done with the ``gplot.corner_res`` function which also accepts an argument ``which_params`` to filter what sub-sets of the free parameters for which to compute the marginal probabilities.\n",
    "So, if our free parameters are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ace8eca3",
   "metadata": {},
   "outputs": [],
   "source": [
    "handler.par_free"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ceed2223",
   "metadata": {},
   "source": [
    "and we want to plot only the free-parameters of the galaxy model (thus marginalising over the noise parameters), we can tell the function to only select parameters with the word ``'galaxy'`` in the name (note the usage of the magic character ``'*'``): ``which_params='galaxy*'``\n",
    "\n",
    "> **Tip:** Internally, GalaPy uses ``getdist`` \n",
    "> (see [online documentation](https://getdist.readthedocs.io/en/latest/index.html) of the package) \n",
    "> to generate corner plots, by passing a dictionary of keyword arguments to the argument ``getdist_settings`` \n",
    "> one can modify the internal defaults of getdist.\n",
    "\n",
    "Once again, conveniently, the function accepts the ``res`` object as first argument: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b1b0243",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, axes = gplot.corner_res(\n",
    "    res, # results object\n",
    "    which_params='galaxy*', # what parameters to plot (use None for all the free parameters)\n",
    "    mark = 'median', # what statistics to plot dashed lines ('median','mean' or 'bestfit')\n",
    "    getdist_settings={\n",
    "        'contours' : [0.68, 0.95],\n",
    "        'smooth_scale_2D' : 0.5,\n",
    "        'fine_bins': 1024\n",
    "    },\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7468e14e",
   "metadata": {},
   "source": [
    "## Tables and TeX\n",
    "\n",
    "It is possible to print values already prepared for TeX math mode with the results of the sampling run.\n",
    "\n",
    "Import the sub-module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "726dc075",
   "metadata": {},
   "outputs": [],
   "source": [
    "import galapy.analysis.funcs as gtabs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20e1055f",
   "metadata": {},
   "source": [
    "The function ``get_parameters_label_strings`` returns a dictionary with the keyword identifying a specific free parameter as key and the corresponding TeX math-mode symbol as value (takes the handler as input): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "933c8410",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtabs.get_parameters_label_strings(handler)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26581f02",
   "metadata": {},
   "source": [
    "to get a collection of summary statistics computed on the free-parameters, already converted to TeX math-mode, use function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cae14d03",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtabs.get_parameters_summary_strings(res, stat_type='quantiles', quantile=(0.025, 0.5, 0.975))"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}