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:


In [1]:
from astromodels import *

list_functions()


/home/giacomov/software/canopy-env/lib/python2.7/site-packages/astromodels-0.1-py2.7.egg/astromodels/functions/functions.py:39: NaimaNotAvailable: The naima package is not available. Models that depend on it will not be available
Out[1]:
nameDescription
bandThe 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)
biasReturn x plus a bias
broken_powerlawA broken power law function
cutoff_powerlawA power law multiplied by an exponential cutoff
cutoff_powerlaw_fluxA cutoff power law having the flux as normalization, which should reduce the correlation among parameters.
gaussianA Gaussian function
identityReturn x
lineA linear function
log_parabolaA log-parabolic function
log_uniform_priorA function which is 1/x on the interval lower_bound - upper_bound and 0 outside the interval. The extremes of the interval are NOT counted as part of the interval. Lower_bound must be >= 0.
powerlawA simple power-law
powerlaw_fluxA simple power-law with the photon flux in a band used as normalization. This will reduce the correlation between the index and the normalization.
sinA sinusodial function
uniform_priorA 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.

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


In [2]:
powerlaw.info()


  • description: A simple power-law
  • formula: $ K~\frac{x}{piv}^{index} $
  • default parameters:
    • K:
      • value: 1.0
      • desc: Normalization (differential flux at the pivot value)
      • min_value: None
      • max_value: None
      • unit:
      • delta: 0.3
      • free: True
    • piv:
      • value: 1.0
      • desc: Pivot value
      • min_value: None
      • max_value: None
      • unit:
      • delta: 0.3
      • free: False
    • index:
      • value: -2.0
      • desc: Photon index
      • min_value: -10.0
      • max_value: 10.0
      • unit:
      • delta: 0.6
      • 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:


In [3]:
powerlaw_instance = powerlaw()

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


In [4]:
powerlaw_instance = powerlaw(K=-2.0, 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:


In [5]:
powerlaw_instance.display()


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

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


In [6]:
print(powerlaw_instance)


  * description: A simple power-law
  * formula: $ K~\frac{x}{piv}^{index} $
  * parameters: 
    * K: 
      * value: -2.0
      * desc: Normalization (differential flux at the pivot value)
      * min_value: None
      * max_value: None
      * unit: 
      * delta: 0.3
      * free: True
    * piv: 
      * value: 1.0
      * desc: Pivot value
      * min_value: None
      * max_value: None
      * unit: 
      * delta: 0.3
      * free: False
    * index: 
      * value: -2.2
      * desc: Photon index
      * min_value: -10.0
      * max_value: 10.0
      * unit: 
      * delta: 0.6
      * 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:


In [7]:
# Modify current value

powerlaw_instance.K = 1.2

# Modify minimum 
powerlaw_instance.K.min_value = -10

# Modify maximum
powerlaw_instance.K.max_value = 15

# We can also modify minimum and maximum at the same time
powerlaw_instance.K.set_bounds(-10, 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()


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

Physical units

Astromodels uses the facility defined in astropy.units to make easier to convert between units during interactive analysis, when assigning to parameters. However, when functions are initialized their parameters do not have units, as it is evident from the .display calls above. They however get assigned units when they are used for something specific, like to represent a spectrum. For example, let's create a point source (see the "Point source tutorial" for more on this):


In [8]:
# Create a powerlaw instance with default values
powerlaw_instance = powerlaw()

# Right now the parameters of the power law don't have any unit
print("Unit of K is [%s]" % powerlaw_instance.K.unit)

# Let's use it as a spectrum for a point source
test_source = PointSource('test_source', ra=0.0, dec=0.0, spectral_shape=powerlaw_instance)

# Now the parameter K has units
print("Unit of K is [%s]" % powerlaw_instance.K.unit)


Unit of K is []
Unit of K is [1 / (cm2 keV s)]

Now if we display the function we can see that other parameters got units as well:


In [9]:
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: None
      • max_value: None
      • unit: 1 / (cm2 keV s)
      • delta: 0.3
      • free: True
    • piv:
      • value: 1.0
      • desc: Pivot value
      • min_value: None
      • max_value: None
      • unit: keV
      • delta: 0.3
      • free: False
    • index:
      • value: -2.0
      • desc: Photon index
      • min_value: -10.0
      • max_value: 10.0
      • unit:
      • delta: 0.6
      • free: True

Note that the index has still no units, as it is intrinsically a dimensionless quantity.

We can now change the values of the parameters using units, or pure floating point numbers. In the latter case, the current unit for the parameter will be assumed:


In [10]:
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)

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

# Express the differential flux at the pivot energy in its default unit
# (currently 1/(keV cm2 s))
powerlaw_instance.K = 122.3

powerlaw_instance.display()


  • description: A simple power-law
  • formula: $ K~\frac{x}{piv}^{index} $
  • parameters:
    • K:
      • value: 122.3
      • desc: Normalization (differential flux at the pivot value)
      • min_value: None
      • max_value: None
      • unit: 1 / (cm2 keV s)
      • delta: 0.3
      • free: True
    • piv:
      • value: 1.0
      • desc: Pivot value
      • min_value: None
      • max_value: None
      • unit: keV
      • delta: 0.3
      • free: False
    • index:
      • value: -2.0
      • desc: Photon index
      • min_value: -10.0
      • max_value: 10.0
      • unit:
      • delta: 0.6
      • free: True

NOTE : using astropy.units in an assigment makes the operation pretty slow. This is hardly noticeable in an interactive settings, but if you put an assigment with units in a for loop or in any other context where it is repeated many times, you might start to notice. For this reason, astromodels allow you to assign directly the value of the parameter in an alternative way, by using the .scaled_value property. This assume that you are providing a simple floating point number, which implicitly uses a specific set of units, which you can retrieve with .scaled_units like this:


In [14]:
print(powerlaw_instance.K.scaled_unit)


1 / (cm2 keV s)

In [15]:
# NOTE: These requires IPython

%timeit powerlaw_instance.K.scaled_value = 122.3 # 1 / (cm2 keV s)
%timeit powerlaw_instance.K = 122.3 / (u.keV * u.cm**2 * u.s)


The slowest run took 4.60 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 3.94 µs per loop
1000 loops, best of 3: 573 µs per loop

As you can see using an assignment with units is more than 100x slower than using .scaled_value. Note that this is a feature of astropy.units, not of astromodels. Thus, do not use assignment with units in computing intensive situations.

Composing functions

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


In [16]:
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})

In [17]:
a_source = PointSource("a_source",l=24.3, b=44.3, spectral_shape=composite)

composite.display()


  • description: (gaussian{1} + powerlaw{2})
  • formula: (no latex formula available)
  • parameters:
    • F_1:
      • 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: 1 / (cm2 s)
      • delta: 0.3
      • free: True
    • mu_1:
      • value: 0.0
      • desc: Central value
      • min_value: None
      • max_value: None
      • unit: keV
      • delta: 0.1
      • free: True
    • sigma_1:
      • value: 1.0
      • desc: standard deviation
      • min_value: 1e-12 keV
      • max_value: None
      • unit: keV
      • delta: 0.3
      • free: True
    • K_2:
      • value: 1.0
      • desc: Normalization (differential flux at the pivot value)
      • min_value: None
      • max_value: None
      • unit: 1 / (cm2 keV s)
      • delta: 0.3
      • free: True
    • piv_2:
      • value: 1.0
      • desc: Pivot value
      • min_value: None
      • max_value: None
      • unit: keV
      • delta: 0.3
      • free: False
    • index_2:
      • value: -2.0
      • desc: Photon index
      • min_value: -10.0
      • max_value: 10.0
      • unit:
      • delta: 0.6
      • free: True

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


In [18]:
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:


In [19]:
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:


In [20]:
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:


In [21]:
print(len(another_composite.parameters)) # 6 parameters
print(len(another_composite2.parameters)) # 9 parameters


6
9

Composing functions as in f(g(x))

Suppose you have two functions (f and g) and you want to compose them in a new function h(x) = f(g(x)). You can accomplish this by using the .of() method:


In [22]:
# Let's get two functions (for example a gaussian and a sin function)
f = gaussian()
g = sin()

# Let's compose them in a composite function h = f(g(x))

h = f.of(g)

# Verify that indeed we have composed the function

# Get a random number between 1 and 10
x = np.random.uniform(1,10)

print (h(x) == f(g(x)))


True

In [ ]:


In [ ]: