Point source plotting basics¶
In 3ML, we distinguish between data and model plotting. Data plots contian real data points and the over-plotted model is (sometimes) folded through an instrument response. Therefore, the x-axis is not always in the same units across instruments if there is energy dispersion.
However, all instuments see the same model and a multi-wavelength fit can be viewed in model space without complication. 3ML uses one interface to plot both MLE and Bayesian fitted models. To demonstrate we will use toy data simulated from a powerlaw and two gaussians for MLE fits and an exponentially cutoff power law with one gaussian for Bayesian fits.
First we load the analysis results:
[1]:
%matplotlib inline
jtplot.style(context="talk", fscale=1, ticks=True, grid=False)
import matplotlib.pyplot as plt
plt.style.use("mike")
import numpy as np
from threeML import *
from threeML.io.package_data import get_path_of_data_file
/Users/jburgess/coding/tml/threeml/threeML/__init__.py:12: UserWarning: No DISPLAY variable set. Using backend for graphics without display (Agg)
warnings.warn("No DISPLAY variable set. Using backend for graphics without display (Agg)")
/Users/jburgess/coding/tml/astromodels/astromodels/functions/function.py:139: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
function_definition = my_yaml.load(dct['__doc__'])
/Users/jburgess/coding/tml/astromodels/astromodels/core/parameter.py:555: UserWarning: We have set the min_value of K to 1e-99 because there was a postive transform
warnings.warn('We have set the min_value of %s to 1e-99 because there was a postive transform' % self.path)
/Users/jburgess/coding/tml/astromodels/astromodels/core/parameter.py:555: UserWarning: We have set the min_value of xc to 1e-99 because there was a postive transform
warnings.warn('We have set the min_value of %s to 1e-99 because there was a postive transform' % self.path)
Configuration read from /Users/jburgess/.threeML/threeML_config.yml
WARNING RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject
WARNING: AstropyDeprecationWarning: The timefunc function is deprecated and may be removed in a future version.
Use astroquery.utils.timer.timefunc instead. [astropy.utils.timer]
WARNING:astropy:AstropyDeprecationWarning: The timefunc function is deprecated and may be removed in a future version.
Use astroquery.utils.timer.timefunc instead.
WARNING: AstropyDeprecationWarning: The timefunc function is deprecated and may be removed in a future version.
Use astroquery.utils.timer.timefunc instead. [astroquery.vo_conesearch.conesearch]
WARNING:astropy:AstropyDeprecationWarning: The timefunc function is deprecated and may be removed in a future version.
Use astroquery.utils.timer.timefunc instead.
WARNING YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
[4]:
#mle1 = load_analysis_results(get_path_of_data_file("datasets/toy_xy_mle1.fits"))
bayes1 = load_analysis_results(get_path_of_data_file("datasets/toy_xy_bayes2.fits"))
WARNING YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
WARNING UserWarning: We have set the min_value of composite.K_1 to 1e-99 because there was a postive transform
---------------------------------------------------------------------------
UnitConversionError Traceback (most recent call last)
~/coding/tml/astromodels/astromodels/core/parameter.py in _set_unit(self, input_unit)
291
--> 292 self.value = self.as_quantity.to(new_unit).value
293
~/.environs/3ml/lib/python3.7/site-packages/astropy-4.0.1.post1-py3.7-macosx-10.15-x86_64.egg/astropy/units/quantity.py in to(self, unit, equivalencies)
688 unit = Unit(unit)
--> 689 return self._new_view(self._to_value(unit, equivalencies), unit)
690
~/.environs/3ml/lib/python3.7/site-packages/astropy-4.0.1.post1-py3.7-macosx-10.15-x86_64.egg/astropy/units/quantity.py in _to_value(self, unit, equivalencies)
660 return self.unit.to(unit, self.view(np.ndarray),
--> 661 equivalencies=equivalencies)
662
~/.environs/3ml/lib/python3.7/site-packages/astropy-4.0.1.post1-py3.7-macosx-10.15-x86_64.egg/astropy/units/core.py in to(self, other, value, equivalencies)
986 else:
--> 987 return self._get_converter(other, equivalencies=equivalencies)(value)
988
~/.environs/3ml/lib/python3.7/site-packages/astropy-4.0.1.post1-py3.7-macosx-10.15-x86_64.egg/astropy/units/core.py in _get_converter(self, other, equivalencies)
917
--> 918 raise exc
919
~/.environs/3ml/lib/python3.7/site-packages/astropy-4.0.1.post1-py3.7-macosx-10.15-x86_64.egg/astropy/units/core.py in _get_converter(self, other, equivalencies)
903 return self._apply_equivalencies(
--> 904 self, other, self._normalize_equivalencies(equivalencies))
905 except UnitsError as exc:
~/.environs/3ml/lib/python3.7/site-packages/astropy-4.0.1.post1-py3.7-macosx-10.15-x86_64.egg/astropy/units/core.py in _apply_equivalencies(self, unit, other, equivalencies)
887 "{} and {} are not convertible".format(
--> 888 unit_str, other_str))
889
UnitConversionError: 'cm2 keV s' and '1 / (cm2 keV s)' are not convertible
During handling of the above exception, another exception occurred:
CannotConvertValueToNewUnits Traceback (most recent call last)
<ipython-input-4-76c4bc7cc062> in <module>
1 #mle1 = load_analysis_results(get_path_of_data_file("datasets/toy_xy_mle1.fits"))
----> 2 bayes1 = load_analysis_results(get_path_of_data_file("datasets/toy_xy_bayes2.fits"))
~/coding/tml/threeml/threeML/analysis_results.py in load_analysis_results(fits_file)
79 if n_results == 1:
80
---> 81 return _load_one_results(f['ANALYSIS_RESULTS', 1])
82
83 else:
~/coding/tml/threeml/threeML/analysis_results.py in _load_one_results(fits_extension)
94 model_dict = my_yaml.load(serialized_model)
95
---> 96 optimized_model = ModelParser(model_dict=model_dict).get_model()
97
98 # Gather statistics values
~/coding/tml/astromodels/astromodels/core/model_parser.py in __init__(self, model_file, model_dict)
90 self._model_dict = model_dict
91
---> 92 self._parse()
93
94 def _parse(self):
~/coding/tml/astromodels/astromodels/core/model_parser.py in _parse(self)
136 else:
137
--> 138 this_parser = SourceParser(source_or_var_name, source_or_var_definition)
139
140 res = this_parser.get_source()
~/coding/tml/astromodels/astromodels/core/model_parser.py in __init__(self, source_name, source_definition)
331 if source_type == POINT_SOURCE:
332
--> 333 self._parsed_source = self._parse_point_source(source_definition)
334
335 elif source_type == EXTENDED_SOURCE:
~/coding/tml/astromodels/astromodels/core/model_parser.py in _parse_point_source(self, pts_source_definition)
420
421 this_point_source = point_source.PointSource(self._source_name, sky_position=this_sky_direction,
--> 422 components=components)
423
424 except:
~/coding/tml/astromodels/astromodels/sources/point_source.py in __init__(self, source_name, ra, dec, spectral_shape, l, b, components, sky_position)
141 for component in list(self._components.values()):
142
--> 143 component.shape.set_units(x_unit, y_unit)
144
145 def __call__(self, x, tag=None):
~/coding/tml/astromodels/astromodels/functions/function.py in set_units(self, x_unit, y_unit, relaxed)
1468 else:
1469
-> 1470 function.set_units(x_unit, y_unit)
1471
1472 @property
~/coding/tml/astromodels/astromodels/functions/function.py in set_units(self, in_x_unit, in_y_unit)
850
851 # Now call the underlying method to set units, which is defined by each function
--> 852 new_units = self._set_units(in_x_unit, in_y_unit)
853
854 # Store the units.
~/coding/tml/astromodels/astromodels/functions/functions.py in _set_units(self, x_unit, y_unit)
325 # The normalization has the same units as the y
326
--> 327 self.K.unit = y_unit
328
329 # noinspection PyPep8Naming
~/coding/tml/astromodels/astromodels/core/parameter.py in _set_unit(self, input_unit)
303
304 raise CannotConvertValueToNewUnits("Cannot convert the value %s from %s to the "
--> 305 "new units %s" % (self.value, self._unit, new_unit_name))
306
307 else:
CannotConvertValueToNewUnits: Cannot convert the value 10.61718649156905 from cm2 keV s to the new units 1 / (cm2 keV s)
Plotting a single analysis result¶
The easiest way to plot is to call plot_point_source_spectra. By default, it plots in photon space with a range of 10-40000 keV evaluated at 100 logrithmic points:
[ ]:
_ = plot_point_source_spectra(mle1,ene_min=1,ene_max=1E3)
Flux and energy units¶
We use astropy units to specify both the flux and energy units. * The plotting routine understands photon, energy (\(F_{\nu}\)) and \(\nu F_{ \nu}\) flux units;
energy units can be energy, frequency, or wavelength
a custom range can be applied.
changing flux units¶
[ ]:
_ = plot_point_source_spectra(mle1,ene_min=1,ene_max=1E3,flux_unit='1/(m2 s MeV)')
_ = plot_point_source_spectra(mle1,ene_min=1,ene_max=1E3,flux_unit='erg/(cm2 day keV)')
_ = plot_point_source_spectra(mle1,ene_min=1,ene_max=1E3,flux_unit='keV2/(cm2 s keV)')
changing energy units¶
[ ]:
_ = plot_point_source_spectra(mle1,
ene_min=.001,
ene_max=1E3,
energy_unit='MeV')
# energy ranges can also be specified in units
_ = plot_point_source_spectra(mle1,
ene_min=1*astropy_units.keV,
ene_max=1*astropy_units.MeV)
_ = plot_point_source_spectra(mle1,
ene_min=1E3*astropy_units.Hz,
ene_max=1E7*astropy_units.Hz)
_ = plot_point_source_spectra(mle1,
ene_min=1E1*astropy_units.nm,
ene_max=1E3*astropy_units.nm,
xscale='linear') # plotting with a linear scale
Plotting components¶
Sometimes it is interesting to see the components in a composite model. We can specify the use_components switch. Here we will use Bayesian results. Note that all features work with MLE of Bayesian results.
[ ]:
_ = plot_point_source_spectra(bayes1,
ene_min=1,
ene_max=1E3,
use_components=True
)
_=plt.ylim(bottom=1)
Notice that the duplicated components have the subscripts n1 and n2. If we want to specify which components to plot, we must use these subscripts.
[ ]:
_ = plot_point_source_spectra(mle1,
flux_unit='erg/(cm2 s keV)',
ene_min=1,
ene_max=1E3,
use_components=True,
components_to_use=['Gaussian_n1','Gaussian_n2'])
_=plt.ylim(bottom=1E-20)
If we want to see the total model with the components, just add total to the components list.
Additionally, we can change the confidence interval for the contours from the default of 1\(\sigma\) (0.68) to 2\(\sigma\) (0.95).
[ ]:
_ = plot_point_source_spectra(bayes1,
flux_unit='erg/(cm2 s keV)',
ene_min=1,
ene_max=1E3,
use_components=True,
components_to_use=['total','Gaussian'],
confidence_level=0.95)
_=plt.ylim(bottom=1E-9)
[ ]:
_ = plot_point_source_spectra(mle1,
flux_unit='erg/(cm2 s keV)',
ene_min=1,
ene_max=1E3,
use_components=True,
fit_cmap='jet', # specify a color map
contour_colors='k', # specify a color for all contours
components_to_use=['total','Gaussian_n2','Gaussian_n1'])
_=plt.ylim(bottom=1E-16)
Additional features¶
Explore the docstring to see all the available options. Default configurations can be altered in the 3ML config file.
Use asymmetric errors and alter the default color map
[ ]:
threeML_config['model plot']['point source plot']['fit cmap'] = 'plasma'
_ = plot_point_source_spectra(mle1, equal_tailed=False)
turn of contours and the legend and increase the number of points plotted
[ ]:
_ = plot_point_source_spectra(mle1, show_legend=False, show_contours=False, num_ene=500)
colors or color maps can be specfied
[ ]:
_ = plot_point_source_spectra(mle1, fit_colors='orange', contour_colors='blue')
Further modifications to plotting style, legend style, etc. can be modified either in the 3ML configuration:
[ ]:
threeML_config['model plot']['point source plot']
or by directly passing dictionary arguments to the the plot command. Examine the docstring for more details!
Plotting multiple results¶
Any number of results can be plotted together. Simply provide them as arguments. You can mix and match MLE and Bayesian results as well as plotting their components.
[ ]:
_ = plot_point_source_spectra(mle1, bayes1,ene_min=1)
_=plt.ylim(bottom=1E-1)
Specify particular colors for each analysis and broaden the contours
[ ]:
_ = plot_point_source_spectra(mle1,
bayes1,
ene_min=1.,
confidence_level=.95,
equal_tailed=False,
fit_colors=['orange','green'],
contour_colors='blue')
_ =plt.ylim(bottom=1E-1)
As with single results, we can choose to plot the components for all the sources.
[ ]:
_ = plot_point_source_spectra(mle1,
bayes1,
ene_min=1.,
use_components=True)
_=plt.ylim(bottom=1E-4)
[ ]: