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.
[1]:
from astromodels.functions.function import (
Function1D,
FunctionMeta,
ModelAssertionViolation,
)
19:55:25 WARNING The naima package is not available. Models that depend on it will not be functions.py:48 available
WARNING The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it functions.py:69 will not be available.
WARNING The ebltable package is not available. Models that depend on it will not be absorption.py:36 available
19:55:26 WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704
WARNING We have set the min_value of F to 1e-99 because there was a postive transform parameter.py:704
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704
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.
[2]:
class Powerlaw(Function1D, metaclass=FunctionMeta):
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
"""
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,
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 in other langauges
What if your model is built from a C++ function and you want to fit that directly to the data? using Cython, pybind, f2py, etc, you can wrap these models and call them easily.
[3]:
def cpp_function_wrapper(a):
# we could wrap a c++ function here
# with cython, pybind11, etc
return a
[4]:
cpp_function_wrapper(2.0)
[4]:
2.0
Now we will define a spectral model that will handle both the unit and non-unit call.
[5]:
import astropy.units as astropy_units
class CppModel(Function1D, metaclass=FunctionMeta):
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
"""
def _set_units(self, x_unit, y_unit):
# The normalization has the same units as the y
self.a.unit = y_unit
def evaluate(self, x, a):
# check is the function is being called 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.0
# call the cython function
flux = cpp_function_wrapper(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
[6]:
cpp_spectrum = CppModel()
from astromodels import PointSource
point_source = PointSource("ps", 0, 0, spectral_shape=cpp_spectrum)
print(point_source(10.0))
point_source(10.0 * astropy_units.keV)
1.0
[6]:
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.
[7]:
from astromodels import Band
import numpy as np
model = Band()
# we won't need to modify the normalization
model.K = 1.0
# 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.
[8]:
from astromodels import TemplateModelFactory
tmf = TemplateModelFactory(
"my_template", "A test template", energies, ["alpha", "xp", "beta"]
)
19:55:28 WARNING Energy unit is not a Quantity instance, so units has not been provided. template_model.py:106 Using keV.
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.
[9]:
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!
[10]:
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.
[11]:
tmf.save_data(overwrite=True)
from astromodels import TemplateModel
reloaded_table_model = TemplateModel("my_template")
[12]:
reloaded_table_model(energies)
[12]:
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])
[ ]: