tutorials/spectrum_models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectral models in Gammapy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "This notebook explains how to use the functions and classes in [gammapy.spectrum.models](https://docs.gammapy.org/dev/spectrum/#module-gammapy.spectrum.models) in order to work with spectral models.\n",
    "\n",
    "The following clases will be used:\n",
    "\n",
    "* [gammapy.spectrum.models.PowerLaw](https://docs.gammapy.org/dev/api/gammapy.spectrum.models.PowerLaw.html)\n",
    "* [gammapy.utils.fitting.Parameter](https://docs.gammapy.org/dev/api/gammapy.utils.fitting.Parameter.html)\n",
    "* [gammapy.utils.fitting.Parameters](https://docs.gammapy.org/dev/api/gammapy.utils.fitting.Parameters.html)\n",
    "* [gammapy.spectrum.models.SpectralModel](https://docs.gammapy.org/dev/api/gammapy.spectrum.models.SpectralModel.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Same procedure as in every script ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from gammapy.spectrum import models\n",
    "from gammapy.utils.fitting import Parameter, Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a model\n",
    "\n",
    "To create a spectral model, instantiate an object of the spectral model class you're interested in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwl = models.PowerLaw()\n",
    "print(pwl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will use default values for the model parameters, which is rarely what you want.\n",
    "\n",
    "Usually you will want to specify the parameters on object creation.\n",
    "One way to do this is to pass `astropy.utils.Quantity` objects like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwl = models.PowerLaw(\n",
    "    index=2.3, amplitude=1e-12 * u.Unit(\"cm-2 s-1 TeV-1\"), reference=1 * u.TeV\n",
    ")\n",
    "print(pwl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you see, some of the parameters have default ``min`` and ``values`` as well as a ``frozen`` flag. This is only relevant in the context of spectral fitting and thus covered in [spectrum_analysis.ipynb](https://github.com/gammapy/gammapy/blob/master/tutorials/spectrum_analysis.ipynb). Also, the parameter errors are not set. This will be covered later in this tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get and set model parameters\n",
    "\n",
    "The model parameters are stored in the ``Parameters`` object on the spectal model. Each model parameter is a ``Parameter`` instance. It has a ``value`` and a ``unit`` attribute, as well as a ``quantity`` property for convenience."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pwl.parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pwl.parameters[\"index\"])\n",
    "pwl.parameters[\"index\"].value = 2.6\n",
    "print(pwl.parameters[\"index\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pwl.parameters[\"amplitude\"])\n",
    "pwl.parameters[\"amplitude\"].quantity = 2e-12 * u.Unit(\"m-2 TeV-1 s-1\")\n",
    "print(pwl.parameters[\"amplitude\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## List available models\n",
    "\n",
    "All spectral models in gammapy are subclasses of ``SpectralModel``. The list of available models is shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "models.SpectralModel.__subclasses__()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting\n",
    "\n",
    "In order to plot a model you can use the ``plot`` function. It expects an energy range as argument. You can also chose flux and energy units as well as an energy power for the plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_range = [0.1, 10] * u.TeV\n",
    "pwl.plot(energy_range, energy_power=2, energy_unit=\"GeV\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameter errors\n",
    "\n",
    "Parameters are stored internally as covariance matrix. There are, however, convenience methods to set individual parameter errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwl.parameters.set_parameter_errors(\n",
    "    {\"index\": 0.2, \"amplitude\": 0.1 * pwl.parameters[\"amplitude\"].quantity}\n",
    ")\n",
    "print(pwl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can access the parameter errors like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwl.parameters.covariance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwl.parameters.error(\"index\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can plot the butterfly using the ``plot_error`` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = pwl.plot_error(energy_range, color=\"blue\", alpha=0.2)\n",
    "pwl.plot(energy_range, ax=ax, color=\"blue\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integral fluxes\n",
    "\n",
    "You've probably asked yourself already, if it's possible to integrated models. Yes, it is. Where analytical solutions are available, these are used by default. Otherwise, a numerical integration is performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwl.integral(emin=1 * u.TeV, emax=10 * u.TeV)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User-defined model\n",
    "\n",
    "Now we'll see how you can define a custom model. To do that you need to subclass ``SpectralModel``. All ``SpectralModel`` subclasses need to have an ``__init__`` function, which sets up the ``Parameters`` of the model and a ``static`` function called ``evaluate`` where the mathematical expression for the model is defined.\n",
    "\n",
    "As an example we will use a PowerLaw plus a Gaussian (with fixed width)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class UserModel(models.SpectralModel):\n",
    "    def __init__(self, index, amplitude, reference, mean, width):\n",
    "        self.parameters = Parameters(\n",
    "            [\n",
    "                Parameter(\"index\", index, min=0),\n",
    "                Parameter(\"amplitude\", amplitude, min=0),\n",
    "                Parameter(\"reference\", reference, frozen=True),\n",
    "                Parameter(\"mean\", mean, min=0),\n",
    "                Parameter(\"width\", width, min=0, frozen=True),\n",
    "            ]\n",
    "        )\n",
    "\n",
    "    @staticmethod\n",
    "    def evaluate(energy, index, amplitude, reference, mean, width):\n",
    "        pwl = models.PowerLaw.evaluate(\n",
    "            energy=energy,\n",
    "            index=index,\n",
    "            amplitude=amplitude,\n",
    "            reference=reference,\n",
    "        )\n",
    "        gauss = amplitude * np.exp(-(energy - mean) ** 2 / (2 * width ** 2))\n",
    "        return pwl + gauss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = UserModel(\n",
    "    index=2,\n",
    "    amplitude=1e-12 * u.Unit(\"cm-2 s-1 TeV-1\"),\n",
    "    reference=1 * u.TeV,\n",
    "    mean=5 * u.TeV,\n",
    "    width=0.2 * u.TeV,\n",
    ")\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_range = [1, 10] * u.TeV\n",
    "model.plot(energy_range=energy_range);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What's next?\n",
    "\n",
    "In this tutorial we learnd how to work with spectral models.\n",
    "\n",
    "Go to [gammapy.spectrum](https://docs.gammapy.org/dev/spectrum/index.html) to learn more."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}