Spectral models are provided via astromodels. For details, visit the astromodels documentation.
The important points are breifly covered below.
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.
In [18]:
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:
Let's look at the simple case of the power law already define in astromodels.
In [2]:
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)
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:
Keep in mind that this is in YAML format.
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.
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,
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.
In [2]:
%load_ext Cython
In [15]:
%%cython --annotate
cpdef cython_function(a):
# we could wrap a c++ function here
return a
Out[15]:
In [16]:
cython_function(2.)
Out[16]:
Now we will define a spectral model that will handle both the unit and non-unit call.
In [34]:
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
In [33]:
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)
Out[33]:
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.
First, we grab a function and make an energy grid.
In [8]:
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.
In [10]:
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.
In [11]:
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!
In [13]:
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.
In [16]:
tmf.save_data(overwrite=True)
from astromodels import TemplateModel
reloaded_table_model = TemplateModel('my_template')
In [17]:
reloaded_table_model(energies)
Out[17]: