See the documents about creating and getting information about functions, point sources and extended sources for details about these operations. Here we only focus on the global model.
In [3]:
from astromodels import *
Now let's define a point source (see "Creating point sources" for details and alternative way to accomplish this):
In [5]:
# Let's start by defining a simple point source with a power law spectrum
pts1 = PointSource('source_1', ra=125.6, dec=-75.3,
spectral_shape=Powerlaw())
# Get some info about what we just created
pts1.display()
Now let's define another source, this time at Galactic Coordinates l = 11.25, b = -22.5, and with two spectral components:
In [7]:
# Another point source with two spectral components
spectrum1 = Powerlaw(K=1e-5, index=-0.75)
component1 = SpectralComponent('synchrotron',spectrum1)
# spectral component with polarization
spectrum2 = Powerlaw(K=2e-6, index=-1.0)
polarization = LinearPolarization(angle=10, degree=50)
component2 = SpectralComponent('IC',spectrum2 ,polarization=polarization)
point_source2 = PointSource('source_2', l=11.25, b=-22.5, components=[component1,component2])
# Have a look at what we just created
point_source2.display()
Now let's create our model, which comprises our two sources:
In [4]:
# Build a model with the two point sources
my_model = Model(pts1, point_source2)
Of course you can use as many sources as needed, like my_model = Model(pts1, pts2, pts3...)
Using the .display() method we can see all free parameters currently in the model:
In [5]:
my_model.display()
A dictionary of free parameters can be obtained like this:
In [6]:
free_parameters = my_model.free_parameters
We can use such dictionary to loop over all free parameters:
In [7]:
for parameter_name, parameter in free_parameters.iteritems():
print("Parameter %s is free" % parameter_name)
More information on a particular source can be obtained like:
In [8]:
my_model.source_1.display()
More information about a particular instance of a function can be obtained like:
In [9]:
my_model.source_1.spectrum.main.powerlaw.display()
Each source and each parameter has a precise path within the model. These paths are displayed by the .display() method of the model instance (see above), and can be used like my_model.[path]. For example:
In [10]:
my_model.display()
In [11]:
# Access the K parameters of the powerlaw spectrum of the main component for source 1:
my_model.source_1.spectrum.main.powerlaw.K = 1e-5
# Access the logK parameters of the spectrum of the IC component of source 2:
my_model.source_2.spectrum.IC.powerlaw.K = 1e-6
The structure of these paths is easy to understand. The model is a tree-like structure. The root of the tree is always the model instance itself. The second level is constituted by the various sources. The structure within a source can be understood by calling the .display method:
In [12]:
my_model.source_1.display()
Each indentation represents one level, so to access the "ra" element we can follow the levels shown by the .display() method:
In [13]:
ra_parameter = my_model.source_1.position.ra
ra_parameter.display()
# NOTE: this is a Parameter instance. To get the position of the source as a
# floating point number, use:
# my_model.source_1.position.get_ra()
# which will work for any source
while to access the K parameter of the power law function we can do:
In [14]:
K_parameter = my_model.source_1.spectrum.main.powerlaw.K
K_parameter.display()
Finally, there is an alternative way of using paths, which might be more adapt for scripts:
In [15]:
my_model['source_1.spectrum.main.powerlaw.K'].display()
You can find much more information in the document "Additional features for scripts and applications".
These fully-qualified paths are unique to each element, are very descriptive and easy to understand. They can always be used and are encouraged in general, but especially in scripts, when the effort spent writing them is paid off in terms of clarity. However, there is an alternative way which might be more convenient in certain situation, especially when models are simple and the chances of getting confused are low. This alternative method is described below.
Exploiting the feature of the python language, we can create names ("shortcuts") for objects:
In [16]:
# Create a "shortcut" for the spectrum of a source
powerlaw_1 = my_model.source_1.spectrum.main.powerlaw
# Now we can change the values of that power law as:
powerlaw_1.logK = -1.2
# GOTCHA: while it is possible to create shortcuts for parameters, it is not encouraged
# Indeed, this will not work:
# logK_1 = my_model.source_1.spectrum.main.powerlaw.logK
# logK_1 = -1.2 # WILL NOT WORK
# In order to use a shortcut for a parameter to change its value, you have to explicitly
# set its property 'value':
# logK_1.value = -1.2 # This will work
Shortcut can point at any point of the tree:
In [17]:
# Create a shortcut of a source
source_1 = my_model.source_1
# Now we can do:
source_1.spectrum.main.powerlaw.index = -2.3
# Create a shortcut for a component
main_component = my_model.source_1.spectrum.main
# Now we can do:
main_component.powerlaw.index = -1.3
If you are ever in doubt of what a particular shortcut stands for, you can always retrieve the full path of the element the shortcut is pointing to like this:
In [18]:
print(main_component.path)
An existing model can be saved to a file with:
In [19]:
# Save the model to a file, overwriting it if already existing
my_model.save('my_model.yml', overwrite=True)
The content of the file is YAML code, which is human-readable and very easy to understand. Let's have a look:
In [20]:
with open('my_model.yml') as yaml_file:
print("".join(yaml_file.readlines()))
Now suppose that you want to load back a file you created in a previous session. You can do it with:
In [21]:
my_model = load_model('my_model.yml')
# Explore the model we just loaded back
my_model.display()
In [22]:
# Now evaluate and plot our models. You need matplotlib for this
import matplotlib.pyplot as plt
%matplotlib inline
# Energies where we want to evaluate the model
e = np.logspace(0,3,100)
# Loop over the sources
for src_name, src in my_model.point_sources.iteritems():
# Loop over the components of each source
for comp_name, component in src.components.iteritems():
# Get the differential flux (in ph/cm2/s)
flux = component.shape(e * u.keV)
# this can also be accomplished with:
# flux = component.powerlaw(e)
# but this requires to know the name of the
# spectral shape which was used
# Plot this component for this source
plt.loglog(e, flux,label="%s of %s" % (component.name, src.name))
plt.legend(loc=0)
_ = plt.xlabel("Energy")
_ = plt.ylabel(r"Flux (ph cm$^{-2}$ s$^{-1}$ keV$^{-1}$")
Sometimes you want to link two parameters of a model so that they have the same value. This can be easily accomplished in astromodels:
In [23]:
# Link the photon index of the first source with the
# photon index of the IC component of the second source
my_model.link(my_model.source_2.spectrum.IC.powerlaw.index,
my_model.source_1.spectrum.main.powerlaw.index)
my_model.display()
Astromodels takes this a step further. Parameters can be linked to each other through any function. The parameters of the linking function become parameters of the model like any other, and can be left free to vary or fixed. For example, let's consider the case where we want the photon index of the IC component of the second source (p2) to be equal to the photon index of the first source (p1) plus a constant. We can link the two parameters with the 'bias' function f(x) = x + k, so that p2(p1) = p1 + k:
In [24]:
# Link the photon indexes through the 'bias' function, i.e.,
# the photon index of the IC component of the second source is fixed to be the
# photon index of the first source plus a constant k
link_function = bias()
my_model.link(my_model.source_2.spectrum.IC.powerlaw.index,
my_model.source_1.spectrum.main.powerlaw.index,
link_function)
# The parameters of the linking function become parameters
# of the model, and are put in the model tree under the parameter they are
# linking.
# In this case the only parameter of the 'bias' function ('k') becomes then
# my_model.source_2.spectrum.IC.powerlaw.logK.bias.k
my_model.display()
If we want to fix say p2 = p1 - 1.2, we can fix k to that:
In [25]:
my_model.source_2.spectrum.IC.powerlaw.index.bias.k = -1.2
my_model.source_2.spectrum.IC.powerlaw.index.bias.k.fix = True
my_model.display()
As another example, we might link the two parameters using a power law function:
In [26]:
my_model.link(my_model.source_2.spectrum.IC.powerlaw.index,
my_model.source_1.spectrum.main.powerlaw.index,
powerlaw())
my_model.display()
We can use arbitrarily complex functions as link function, if needed (see "Creating and modifying functions" for more info on how to create composite functions):
In [27]:
# A random composite function (see "Creating and modifying functions" for more info)
crazy_link = powerlaw() + gaussian()
my_model.link(my_model.source_2.spectrum.IC.powerlaw.index,
my_model.source_1.spectrum.main.powerlaw.index,
crazy_link)
my_model.display()
In astromodels parameters can become functions of independent variables such as time. This is accomplished in a way which is similar to the procedure to link parameters described above. First, let's create an independent variable. An IndependentVariable instance is created and added to the model like this:
In [28]:
# Add the time as an independent variable of the model
time = IndependentVariable("time",0.0, unit='s')
my_model.add_independent_variable(time)
The IndependentVariable instance is inserted at the root of the model tree. In this case, can be accessed as:
In [29]:
my_model.time.display()
We can now link any parameter to be a function of time, like this:
In [30]:
# First define the function. In this case, a linear function ax+b
law = line(a=-0.02,b=-2.0)
# Now link the index of the sync. component of source_2 to be law(t)
# i.e., index = law(t) = a*t + b
my_model.link(my_model.source_2.spectrum.synchrotron.powerlaw.index,
time,
law)
In [31]:
# This would show the link:
# my_model.display()
# Now changing the value of time will change the value of the parameter
# according to the law. For example, let's loop over 10 s and print
# the value of the parameter
# Reset time
my_model.time = 0.0
for i in range(10):
my_model.time = my_model.time.value + 1.0 * u.s
print("At time %s the value of the parameter is %s" % (my_model.time.value,
my_model.source_2.spectrum.synchrotron.powerlaw.index.value))
In [32]:
# Now plot the synch. spectrum of the source at different times
# (you will need matplotlib for this)
import matplotlib.pyplot as plt
# This is needed only in the ipython notebook
%matplotlib inline
# Prepare 100 logarithmically distributed energies between 1 and 100 keV
energies = np.logspace(0,2,100)
# Compute and plot the sync. spectrum every 10 s between 0 and 50 seconds
times = np.linspace(0,50,6)
my_model.time = 0.0
for tt in times:
my_model.time = tt
plt.loglog(energies, my_model.source_2.spectrum.synchrotron(energies),label='t = %s' % my_model.time.value)
plt.legend(loc=1,ncol=2)
_ = plt.xlabel("Energy (keV)")
_ = plt.ylabel(r"Differential flux (ph. cm$^{-2}$ s$^{-1}$ keV$^{-1}$)")
In [ ]: