In [1]:
from astromodels import *
list_functions()
Out[1]:
If you need more info about a function, you can obtain it by using:
In [2]:
powerlaw.info()
Note that you don't need to create an instance in order to call the info() method.
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.
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()
It is also possible to get the text-only representation by simply printing the object like this:
In [6]:
print(powerlaw_instance)
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 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()
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)
Now if we display the function we can see that other parameters got units as well:
In [9]:
powerlaw_instance.display()
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()
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)
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)
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.
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)
In [17]:
a_source = PointSource("a_source",l=24.3, b=44.3, spectral_shape=composite)
composite.display()
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)
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)
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)
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
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)))
In [ ]:
In [ ]: