Functions tutorial

In astromodels functions can be used as spectral shapes for sources, or to describe time-dependence, phase-dependence, or links among parameters.

To get the list of available functions just do:

%%capture
from astromodels import *

list_functions()

name	Description
Asymm_Gaussian_on_sphere	A bidimensional Gaussian function on a sphere (in spherical coordinates)\nsee https://en.wikipedia.org/wiki/Gaussian_function#Two-dimensional_Gaussian_function
Band	Band model from Band et al., 1993, parametrized with the peak energy
Band_Calderone	The Band model from Band et al. 1993, implemented however in a way which reduces the covariances between the parameters (Calderone et al., MNRAS, 448, 403C, 2015)
Band_grbm	Band model from Band et al., 1993, parametrized with the cutoff energy
Beta	A beta distribution function
Blackbody	A blackbody function
Broken_powerlaw	A broken power law function
Cauchy	The Cauchy distribution
Constant	Return k
Continuous_injection_diffusion	Positron and electrons diffusing away from the accelerator
...	...
SpatialTemplate_2D	User input Spatial Template. Expected to be normalized to 1/sr
StepFunction	A function which is constant on the interval lower_bound - upper_bound and 0 outside the interval. The extremes of the interval are counted as part of the interval.
StepFunctionUpper	A function which is constant on the interval lower_bound - upper_bound and 0 outside the interval. The upper interval is open.
Super_cutoff_powerlaw	A power law with a super-exponential cutoff
TbAbs	Photometric absorption (Tbabs implementation), f(E) = exp(- NH * sigma(E)) contributed by Dominique Eckert
TemplateModel	A template model
Truncated_gaussian	A truncated Gaussian function defined on the interval between the lower_bound (a) and upper_bound (b)
Uniform_prior	A function which is constant on the interval lower_bound - upper_bound and 0 outside the interval. The extremes of the interval are counted as part of the interval.
WAbs	Photometric absorption (Wabs implementation), f(E) = exp(- NH * sigma(E)) contributed by Dominique Eckert
ZDust	Extinction by dust grains from Pei (1992), suitable for IR, optical and UV energy bands, including the full energy ranges of the Swift UVOT and XMM-Newton OM detectors. Three models are included which characterize the extinction curves of (1) the Milky Way, (2) the LMC and (3) the SMC. The models can be modified by redshift and can therefore be applied to extragalactic sources. The transmission is set to unity shortward of 912 Angstroms in the rest frame of the dust. This is incorrect physically but does allow the model to be used in combination with an X-ray photoelectric absorption model such as phabs. Parameter 1 (method) describes which extinction curve (MW, LMC or SMC) will be constructed and should never be allowed to float during a fit. The extinction at V, A(V) = E(B-V) x Rv. Rv should typically remain frozen for a fit. Standard values for Rv are MW = 3.08, LMC = 3.16 and SMC = 2.93 (from table 2 of Pei 1992), although these may not be applicable to more distant dusty sources.
_ComplexTestFunction	A useless function to be used during automatic tests

If you need more info about a function, you can obtain it by using:

Gaussian.info()

• description: A Gaussian function
• formula: $ K \frac{1}{\sigma \sqrt{2 \pi}}\exp{\frac{(x-\mu)^2}{2~(\sigma)^2}} $
• default parameters:
  • F:
    • value: 1.0
    • desc: Integral between -inf and +inf. Fix this to 1 to obtain a Normal distribution
    • min_value: None
    • max_value: None
    • unit:
    • is_normalization: False
    • delta: 0.1
    • free: True
  • mu:
    • value: 0.0
    • desc: Central value
    • min_value: None
    • max_value: None
    • unit:
    • is_normalization: False
    • delta: 0.1
    • free: True
  • sigma:
    • value: 1.0
    • desc: standard deviation
    • min_value: 1e-12
    • max_value: None
    • unit:
    • is_normalization: False
    • delta: 0.1
    • free: True

Note that you don’t need to create an instance in order to call the info() method.

Creating functions

Functions can be created in two different ways. We can create an instance with the default values for the parameters like this:

powerlaw_instance = Powerlaw()

or we can specify on construction specific values for the parameters:

powerlaw_instance = Powerlaw(K=0.01, index=-2.2)

If you don’t remember the names of the parameters just call the .info() method as in powerlaw.info() as demonstrated above.

Getting information about an instance

Using the .display() method we get a representation of the instance which exploits the features of the environment we are using. If we are running inside a IPython notebook, a rich representation with the formula of the function will be displayed (if available). Otherwise, in a normal terminal, the latex formula will not be rendered:

powerlaw_instance.display()

• description: A simple power-law
• formula: $ K~\frac{x}{piv}^{index} $
• parameters:
  • K:
    • value: 0.01
    • desc: Normalization (differential flux at the pivot value)
    • min_value: 1e-30
    • max_value: 1000.0
    • unit:
    • is_normalization: True
    • delta: 0.1
    • free: True
  • piv:
    • value: 1.0
    • desc: Pivot value
    • min_value: None
    • max_value: None
    • unit:
    • is_normalization: False
    • delta: 0.1
    • free: False
  • index:
    • value: -2.2
    • desc: Photon index
    • min_value: -10.0
    • max_value: 10.0
    • unit:
    • is_normalization: False
    • delta: 0.20099999999999998
    • free: True

It is also possible to get the text-only representation by simply printing the object like this:

print(powerlaw_instance)

  * description: A simple power-law
  * formula: $ K~\frac{x}{piv}^{index} $
  * parameters:
    * K:
      * value: 0.01
      * desc: Normalization (differential flux at the pivot value)
      * min_value: 1.0e-30
      * max_value: 1000.0
      * unit: ''
      * is_normalization: true
      * delta: 0.1
      * free: true
    * piv:
      * value: 1.0
      * desc: Pivot value
      * min_value: null
      * max_value: null
      * unit: ''
      * is_normalization: false
      * delta: 0.1
      * free: false
    * index:
      * value: -2.2
      * desc: Photon index
      * min_value: -10.0
      * max_value: 10.0
      * unit: ''
      * is_normalization: false
      * delta: 0.20099999999999998
      * free: true

Note: the .display() method of an instance displays the current values of the parameters, while the .info() method demonstrated above (for which you don’t need an instance) displays the default values of the parameters.

Modifying parameters

Modifying a parameter of a function is easy:

# Modify current value

powerlaw_instance.K = 1.2

# Modify minimum
powerlaw_instance.K.min_value = 0.5

# Modify maximum
powerlaw_instance.K.max_value = 15

# We can also modify minimum and maximum at the same time
powerlaw_instance.K.bounds = (0.5, 15)

# Modifying the delta for the parameter
# (which can be used by downstream software for fitting, for example)
powerlaw_instance.K.delta = 0.25

# Fix the parameter
powerlaw_instance.K.fix = True

# or equivalently
powerlaw_instance.K.free = False

# Free it again
powerlaw_instance.K.fix = False

# or equivalently
powerlaw_instance.K.free = True

# We can verify what we just did by printing again the whole function as shown above,
# or simply printing the parameter:
powerlaw_instance.K.display()

21:46:40 WARNING   We have set the min_value of Powerlaw.K to 1e-99 because there was a postive    parameter.py:704
                  transform                                                                                        

Parameter K = 1.2 [] (min_value = 0.5, max_value = 15.0, delta = 0.25, free = True)

Using physical units

Astromodels uses the facility defined in astropy.units to make easier to convert between units during interactive analysis, when assigning to parameters. In order for functions to be aware of their units, they must be part of a `Source. Let’s create one:

powerlaw_instance = Powerlaw()

point_source = PointSource("my_point_source",ra=0,dec=0, spectral_shape=powerlaw_instance)

Now we can see the units

powerlaw_instance.display()

• description: A simple power-law
• formula: $ K~\frac{x}{piv}^{index} $
• parameters:
  • K:
    • value: 1.0
    • desc: Normalization (differential flux at the pivot value)
    • min_value: 1e-30
    • max_value: 1000.0
    • unit: keV-1 s-1 cm-2
    • is_normalization: True
    • delta: 0.1
    • free: True
  • piv:
    • value: 1.0
    • desc: Pivot value
    • min_value: None
    • max_value: None
    • unit: keV
    • is_normalization: False
    • delta: 0.1
    • free: False
  • index:
    • value: -2.01
    • desc: Photon index
    • min_value: -10.0
    • max_value: 10.0
    • unit:
    • is_normalization: False
    • delta: 0.20099999999999998
    • free: True

powerlaw_instance.x_unit

$\mathrm{keV}$

powerlaw_instance.y_unit

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

import astropy.units as u

# Express the differential flux at the pivot energy in 1 / (MeV cm2 s)

powerlaw_instance.K = (122.3 / (u.MeV * u.cm * u.cm * u.s))

print(powerlaw_instance.K)

# Express the differential flux at the pivot energy in 1 / (GeV m2 s)
powerlaw_instance.K = (122.3 / (u.GeV * u.m * u.m * u.s))


print(powerlaw_instance.K)

Parameter K = 0.1223 [1 / (keV s cm2)]
(min_value = 1e-30, max_value = 1000.0, delta = 0.1, free = True)
Parameter K = 1.2230000000000013e-08 [1 / (keV s cm2)]
(min_value = 1e-30, max_value = 1000.0, delta = 0.1, free = True)

We see that astromodels does the unit conversion for us in the background!

However, astropy units are very slow and we would not want to deal with them or set parameters with units during a fit. Thus, if you do not specify units when setting a parameter, the value is assumed to have the units specified in the construction of the function.

%timeit powerlaw_instance.K = 1

5.24 µs ± 9.92 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

%timeit powerlaw_instance.K = (122.3 / (u.MeV * u.cm * u.cm * u.s))

90.1 µs ± 1.63 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

As you can see using astropy.units requires about 10x more than using a plain assignment. In an interactive analysis you are unlikely to notice the difference, but if you use units in a loop or during a fit this slow-down will add up an become very noticeable. Note that this is a feature of astropy.units, not of astromodels.

Composing functions

We can create arbitrary complex functions by combining “primitive” functions using the normal math operators:

composite = Gaussian() + Powerlaw()

# Instead of the usual .display(), which would print all the many parameters,
# let's print just the description of the new composite functions:
print(composite.description)

(Gaussian{1} + Powerlaw{2})

These expressions can be as complex as needed. For example:

crazy_function = 3 * Sin() + Powerlaw()**2 * (5+Gaussian()) / 3.0

print(crazy_function.description)

((Sin{1} * 3) + (((Powerlaw{2} ** 2) * (Gaussian{3} + 5)) / 3.0))

The numbers between {} enumerate the unique functions which constitute a composite function. This is useful because composite functions can be created starting from pre-existing instances of functions, in which case the same instance can be used more than once. For example:

a_powerlaw = Powerlaw()
a_sin = Sin()

another_composite = 2 * a_powerlaw + (3 + a_powerlaw) * a_sin

print(another_composite.description)

((Powerlaw{1} * 2) + ((Powerlaw{1} + 3) * Sin{2}))

In this case the same instance of a power law has been used twice. Changing the value of the parameters for “a_powerlaw” will affect also the second part of the expression. Instead, by doing this:

another_composite2 = 2 * Powerlaw() + (3 + Powerlaw()) * Sin()

print(another_composite2.description)

((Powerlaw{1} * 2) + ((Powerlaw{2} + 3) * Sin{3}))

we will end up with two independent sets of parameters for the two power laws. The difference can be seen immediately from the number of parameters of the two composite functions:

print(len(another_composite.parameters)) # 6 parameters
print(len(another_composite2.parameters)) # 9 parameters

6
9

Creating custom functions

One of the most powerful aspects of astromodels is the ability to quickly build custom functions on the fly. The source code for a function 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 1D function in astromodels, 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 a python meta type 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, metaclass=FunctionMeta):
        r"""
        description :
            A  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


        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 function. 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.

Optionally, there can be an additional properties category that we will cover later.

Keep in mind that this is in YAML format.

Set units

astromodels keeps 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 argument 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. For 2D and 3D functions, the functions have y and z for their first arguments as well. ** x, y, and z are reserved names in functions**.

A functions 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 Thus, remember to import any local models you used for an analysis before trying to reload that analysis.

Custom functions 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, pybind11, f2py, etc, you can wrap these models and call them easily.

def cpp_function_wrapper(a):
    # we could wrap a c++ function here
    # with cython, pybind11, etc

    return a

cpp_function_wrapper(2.)

2.0

Now we will define a astromodels function that will handle both the unit and non-unit call.

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.

            # 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

cpp_spectrum = CppModel()

from astromodels import PointSource

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

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

1.0

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

Advanced functions

We are not limited to functions that can only take numerical parameters as arguments. We can also link other astromodels function into a function to expand its abilities. ### Properties

Let’s create a function that uses text based switches to alter its functionality

class SwitchFunction(Function1D,metaclass=FunctionMeta):
        r"""
        description :
            A demo function that can alter its state
        latex : $$

        parameters :
            a :
                desc : Normalization (differential flux)
                initial value : 1.0
                is_normalization : True
                min : 1e-30
                max : 1e3
                delta : 0.1
        properties:
            switch:
                desc: a switch for functions
                initial value: powerlaw
                allowed values:
                    - powerlaw
                    - cosine
                function: _say_hello
        """
        def _say_hello(self):
        # called when we set the value of switch

            print(self.switch.value)


        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):

            if self.switch.value == "powerlaw":

                return a * np.powerlaw(x,-2)

            elif self.switch.value == "cosine":

                return a * np.cos(x)

We have added a text parameter called switch and specified the allowed values that it can take on. We can of course allow it to take on any value. Additionally, we have specified a function to call whenever we change the value. This allows use to do things like read in a table from a dictionary or load a file, etc. Properties can be set in the constructor of a function. They behave just like parameters except that they do not participate in the function call. Thus, their state is saved whenever you serialize the model to disk.

f = SwitchFunction()

f.switch = 'cosine'

powerlaw
cosine

In the docstring, one can also specify defer: True which allows you to not set a value until instancing an object. This is useful if you have a model that reads in file at runtime, but the file name is not known until then. Check out the source code of astromodels to see how properties can be used to expand the functionality of your custom models. For example, the absorption models such as TbAbs take advantage of this to set their abundance tables.

Linking functions

What if you need to call another astromodels function from inside your custom function? This is achieved by linking functions. For example, let’s create a convolutional model that redshifts the energies of and model linked to it:

class Redshifter(Function1D, metaclass=FunctionMeta):
    r"""
    description :
        a function that can redshift the energies of any 1D function

    latex: not available

    parameters :
        redshift :
            desc : the redshift
            initial value : 0
            min : 0

    """


    def _set_units(self, x_unit, y_unit):

        self.redshift.unit = astropy_units.dimensionless_unscaled


    def set_linked_function(self, function):
         # this is an optional helper to
             # ease in the setting of the function



        if "func" in self._external_functions:

            # since we only want to link one function
                    # we unlink if we have linked before

            self.unlink_external_function("func")

        # this allows use to link in a function
            # with an internal name 'func'

        self.link_external_function(function, "func")

        self._linked_function = function

    def get_linked_function(self):

        # linked functions are stored in a dictionary
            # but you

        return self._external_functions["func"]


    linked_function = property(
        get_linked_function,
        set_linked_function,
        doc="""Get/set linked function""",
    )


    def evaluate(self, x, redshift):


        # we can call the function here
        return self._linked_function(x * (1 + redshift))

With this function, whenever we set the linked_function property to another astromodels function, its call will return that function redshifted.

p = Powerlaw()
rs = Redshifter(redshift=1)

rs.linked_function = p

print(p(10.))

print(rs(10.))

0.009772372209558112
0.0024262173759824015

We have added a lot of syntax sugar to make it easier for users to handle the function, but every function in astromodels has the members f.link_external_function(func, 'internal_name') and f.unlink_external_function('internal_name'). You can link as many functions as needed and they are accessed via an internal dictionary self._extranal_functions. As long as all functions used are part of the model, all the linking is saved when a model is saved to disk allowing you to restore all the complexity you built.