Spectral Models

Spectral models are provided via astromodels. For details, visit the astromodels documentation.

The important points are breifly covered below.

Building Custom Models

One of the most powerful aspects of astromodels and 3ML is the ability to quickly build custom models on the fly. The source code for a model can be pure python, FORTRAN linked via f2py, C++ linked via cython, etc. Anything that provides a python function can be used to fit data.

To build a custom spectral model in 3ML, we need to import a few things that will allow astromodels to recognize your model.

from astromodels.functions.function import Function1D, FunctionMeta, ModelAssertionViolation

Function1D is the base class for 1D spectral models and FunctionMeta is ABC class that ensures all the needed parts of a model are in the class as well as making the class function as it should.

There are three basic parts to declaring a model:

• the docstring

• the units setter

• the evaluate function

Let’s look at the simple case of the power law already define in astromodels.

class Powerlaw(Function1D):
        r"""
        description :
            A simple power-law
        latex : $ K~\frac{x}{piv}^{index} $
        parameters :
            K :
                desc : Normalization (differential flux at the pivot value)
                initial value : 1.0
                is_normalization : True
                transformation : log10
                min : 1e-30
                max : 1e3
                delta : 0.1
            piv :
                desc : Pivot value
                initial value : 1
                fix : yes
            index :
                desc : Photon index
                initial value : -2
                min : -10
                max : 10

        """

        __metaclass__ = FunctionMeta

        def _set_units(self, x_unit, y_unit):
            # The index is always dimensionless
            self.index.unit = astropy_units.dimensionless_unscaled

            # The pivot energy has always the same dimension as the x variable
            self.piv.unit = x_unit

            # The normalization has the same units as the y

            self.K.unit = y_unit

        # noinspection PyPep8Naming
        def evaluate(self, x, K, piv, index):

            xx = np.divide(x, piv)

            return K * np.power(xx, index)

The docstring

We have used the docstring interface to provide a YAML description of the model. This sets up the important information used in the fitting process and record keeping. The docstring has three parts:

• description

  • The description is a text string that provides readable info about the model. Nothing fancy, but good descriptions help to inform the user.

• latex

  • If the model is analytic, a latex formula can be included

• parameters

  • For each parameter, a description and initial value must be included. Transformations for fitting, min/max values and fixing the parameter can also be described here.

Keep in mind that this is in YAML format.

Set units

3ML and astromodels keep track of units for you. However, a model must be set up to properly describe the units with astropy’s unit system. Keep in mind that models are fit with a differential photon flux,

\[\frac{d N_p}{dA dt dE}\]

so your units should reflect this convention. Therefore, proper normalizations should be taken into account.

Evaluate

This is where the function is evaluated. The first argumument must be called x and the parameter names and ordering must reflect what is in the docstring. Any number of operations can take place inside the evaluate call, but remember that the return must be in the form of a differential photon flux.

A model is defined in a python session. If you save the results of a fit to an AnalysisResults file and try to load this file without loading this model, you will get a error,

Custom models with Cython

What if your model is built from a C++ function and you want to fit that directly to the data? Cython provides and interface for both importing C++ functions as well as typing python code in order to provide speed. However, there are some caveats that we must consider when using Cython (and other interfaces that do not play nice with units). First, we need to talk about how 3ML performs a fit. Since astromodels keeps track of units, spectral models called with and array of energies/wavelengths will return an astropy Quantity. Astropy Quantities are very slow to evaluate, and we want fits to be fast. Thus, models can be called with (slow) and without (fast) units. Cython calls will crash if an array of Quantities are passed to the function. So, we will have to provide a work around.

Let’s create a quick Cython function. The same will apply to Cython functions that are wrapping C++ function calls.

%load_ext Cython

%%cython --annotate

cpdef cython_function(a):
    # we could wrap a c++ function here
    return a

Generated by Cython 0.25.2

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

 1: 

+2: cpdef cython_function(a):

static PyObject *__pyx_pw_46_cython_magic_fe51d6c672d899f52c26bbea574d1b98_1cython_function(PyObject *__pyx_self, PyObject *__pyx_v_a); /*proto*/
static PyObject *__pyx_f_46_cython_magic_fe51d6c672d899f52c26bbea574d1b98_cython_function(PyObject *__pyx_v_a, CYTHON_UNUSED int __pyx_skip_dispatch) {
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("cython_function", 0);
/* … */
  /* function exit code */
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_fe51d6c672d899f52c26bbea574d1b98_1cython_function(PyObject *__pyx_self, PyObject *__pyx_v_a); /*proto*/
static PyObject *__pyx_pw_46_cython_magic_fe51d6c672d899f52c26bbea574d1b98_1cython_function(PyObject *__pyx_self, PyObject *__pyx_v_a) {
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("cython_function (wrapper)", 0);
  __pyx_r = __pyx_pf_46_cython_magic_fe51d6c672d899f52c26bbea574d1b98_cython_function(__pyx_self, ((PyObject *)__pyx_v_a));

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_fe51d6c672d899f52c26bbea574d1b98_cython_function(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_a) {
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("cython_function", 0);
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_1 = __pyx_f_46_cython_magic_fe51d6c672d899f52c26bbea574d1b98_cython_function(__pyx_v_a, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_r = __pyx_t_1;
  __pyx_t_1 = 0;
  goto __pyx_L0;

  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_AddTraceback("_cython_magic_fe51d6c672d899f52c26bbea574d1b98.cython_function", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 3:     # we could wrap a c++ function here

+4:     return a

  __Pyx_XDECREF(__pyx_r);
  __Pyx_INCREF(__pyx_v_a);
  __pyx_r = __pyx_v_a;
  goto __pyx_L0;

cython_function(2.)

2.0

Now we will define a spectral model that will handle both the unit and non-unit call.

import astropy.units as astropy_units

class CythonModel(Function1D):
        r"""
        description :
            A spectral model wrapping a cython function
        latex : $$
        parameters :
            a :
                desc : Normalization (differential flux)
                initial value : 1.0
                is_normalization : True
                min : 1e-30
                max : 1e3
                delta : 0.1
        """

        __metaclass__ = FunctionMeta

        def _set_units(self, x_unit, y_unit):

            # The normalization has the same units as the y

            self.a.unit = y_unit

        # noinspection PyPep8Naming
        def evaluate(self, x, a):

            # check is the function is being alled with units

            if isinstance(a, astropy_units.Quantity):

                # get the values
                a_ = a.value

                # save the unit
                unit_ = self.y_unit

            else:

                # we do not need to do anything here
                a_ = a

                # this will basically be ignored
                unit_ = 1.

            # call the cython function
            flux = cython_function(a_)

            # add back the unit if needed
            return flux * unit_

We can check the unit and non-unit call by making a point source and evaluating it

cython_spectrum = CythonModel()

from astromodels import PointSource

point_source = PointSource('ps',0,0,spectral_shape=cython_spectrum)

print(point_source(10.))
point_source(10. * astropy_units.keV)

1.0

$1 \; \mathrm{\frac{1}{keV\,s\,cm^{2}}}$

Template (Table) Models

3ML (via astromodels) provides the ability to construct models in tabluated form. This is very useful for models that are from numerical simualtions. While in other software special care must be taken to construct table models into FITS files, in 3ML the construction of the table model is taken care of for you. Here is an example of how to build a template model from a pre-existing function.

Constructing a template

First, we grab a function and make an energy grid.

from astromodels import Band
import numpy as np

model = Band()

# we won't need to modify the normalization
model.K = 1.


# if no units are provided for the energy grid, keV will be assumed!
energies = np.logspace(1, 3, 50)

Now we define a template model factory. This takes a name, a description, the energy grid and an array of parameter names as input.

from astromodels import TemplateModelFactory
tmf = TemplateModelFactory('my_template', 'A test template', energies, ['alpha', 'xp', 'beta'])

Now, define our grid in parameter space. While we are using a function here, this grid could be from a text file, a database of simulations, etc. We then assign these grid points to the template model factory.

alpha_grid = np.linspace(-1.5, 1, 15)
beta_grid = np.linspace(-3.5, -1.6, 15)
xp_grid = np.logspace(1, 3, 20)



tmf.define_parameter_grid('alpha', alpha_grid)
tmf.define_parameter_grid('beta', beta_grid)
tmf.define_parameter_grid('xp', xp_grid)

Finally, we loop over our grid and set the interpolation data to the template model factory. The units of the fluxes must be a differential photon flux!

for a in alpha_grid:

    for b in beta_grid:

        for xp in xp_grid:

            # change our model parameters
            model.alpha = a
            model.beta = b
            model.xp = xp

            tmf.add_interpolation_data(model(energies), alpha=a, xp=xp, beta=b)

We can now save our model to disk. The formal is an HDF5 file which is saved to the astromodel data directory (~/.astromodels/data). The HDF5 file can easily be passed around to other users as all information defining the model is stored in the file. The other user would place the file in thier astromodels data directory.

tmf.save_data(overwrite=True)

from astromodels import TemplateModel

reloaded_table_model = TemplateModel('my_template')

reloaded_table_model(energies)

array([  1.49086966e+00,   1.43117010e+00,   1.37151134e+00,
         1.31187019e+00,   1.25223305e+00,   1.19259749e+00,
         1.13297396e+00,   1.07338763e+00,   1.01388024e+00,
         9.54512018e-01,   8.95363460e-01,   8.36537051e-01,
         7.78158661e-01,   7.20378565e-01,   6.63371907e-01,
         6.07338417e-01,   5.52501210e-01,   4.99104446e-01,
         4.47409687e-01,   3.97690776e-01,   3.50227162e-01,
         3.05295643e-01,   2.63160667e-01,   2.24063443e-01,
         1.88210330e-01,   1.55761153e-01,   1.26818273e-01,
         1.01829853e-01,   8.10538848e-02,   6.39354690e-02,
         5.03106841e-02,   3.95893700e-02,   3.11527907e-02,
         2.45140645e-02,   1.92900651e-02,   1.51793110e-02,
         1.19445674e-02,   9.39915463e-03,   7.39617472e-03,
         5.82003410e-03,   4.57977241e-03,   3.60381313e-03,
         2.83583286e-03,   2.23151083e-03,   1.75597111e-03,
         1.38176992e-03,   1.08731180e-03,   8.55603333e-04,
         6.73272437e-04,   5.29796643e-04])