4754829731

Build Type

push

github

Committed by Johannes Buchner

Commit Message

update docs with new code

Run Details

3948 of 7059 relevant lines covered (55.93%)

0.63 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

92.94

/examples/xspec/example_simplest.py

"""
Example of doing BXA in X-spec
"""
import bxa.xspec as bxa
from xspec import *

Fit.statMethod = 'cstat'
Plot.xAxis = 'keV'

s = Spectrum('example-file.fak')
s.ignore("**"); s.notice("0.2-8.0")

# set model and parameters ranges
#                         val, delta, min, bottom, top, max
m = Model("pow")
m.powerlaw.norm.values = ",,1e-10,1e-10,1e1,1e1" # 10^-10 .. 10
m.powerlaw.PhoIndex.values = ",,1,1,3,3"         #     1 .. 3

# define prior
transformations = [
        # uniform prior for Photon Index (see other example for 
        # something more advanced)
        bxa.create_uniform_prior_for( m, m.powerlaw.PhoIndex),
        # jeffreys prior for scale variable
        bxa.create_jeffreys_prior_for(m, m.powerlaw.norm),
        # and possibly many more parameters here
]
print('running analysis ...')
# where to store intermediate and final results? this is the prefix used
outputfiles_basename = 'simplest/'
solver = bxa.BXASolver(transformations=transformations, outputfiles_basename=outputfiles_basename)
results = solver.run(resume=True)
print('running analysis ... done!')


import matplotlib.pyplot as plt

print('creating plot of posterior predictions against data ...')
plt.figure()

from xspec import Plot
from bxa.xspec.solver import set_parameters, XSilence
from ultranest.plot import PredictionBand

band = None
Plot.background = True

with XSilence():
        Plot.device = '/null'
        # plot models
        for row in solver.posterior[:400]:
                set_parameters(values=row, transformations=solver.transformations)
                Plot('counts')
                if band is None:
                        band = PredictionBand(Plot.x())
                band.add(Plot.model())

band.shade(alpha=0.5, color='blue')
band.shade(q=0.495, alpha=0.1, color='blue')
band.line(color='blue', label='convolved model')

plt.scatter(Plot.x(), Plot.y(), label='data')
#plt.errorbar(
#        x=Plot.x(), xerr=Plot.xErr(), y=Plot.y(), yerr=Plot.yErr(), 
#        marker='o', label='data')

if Plot.xAxis == 'keV':
        plt.xlabel('Energy [keV]')
elif Plot.xAxis == 'channel':
        plt.xlabel('Channel')
plt.ylabel('Counts/s/cm$^2$')
plt.savefig(outputfiles_basename + 'convolved_posterior_direct.pdf', bbox_inches='tight')
plt.close()



print('creating a rebinned plot to check the model ...')
plt.figure()
data = solver.posterior_predictions_convolved(nsamples=100)
print('binning for plot...')
binned = bxa.binning(outputfiles_basename=outputfiles_basename, 
        bins = data['bins'], widths = data['width'], 
        data = data['data'], models = data['models'])
for point in binned['marked_binned']:
        plt.errorbar(marker='o', zorder=-1, **point)
plt.xlim(binned['xlim'])
plt.ylim(binned['ylim'][0], binned['ylim'][1]*2)
plt.gca().set_yscale('log')
if Plot.xAxis == 'keV':
        plt.xlabel('Energy [keV]')
elif Plot.xAxis == 'channel':
        plt.xlabel('Channel')
plt.ylabel('Counts/s/cm$^2$')
print('saving plot...')
plt.legend()
plt.savefig(outputfiles_basename + 'convolved_posterior.pdf', bbox_inches='tight')
plt.close()




print('creating plot of posterior predictions ...')
plt.figure()
solver.posterior_predictions_unconvolved(nsamples=100)
ylim = plt.ylim()
# 3 orders of magnitude at most
plt.ylim(max(ylim[0], ylim[1] / 1000), ylim[1])
plt.yscale('log')
plt.xscale('log')
if Plot.xAxis == 'keV':
        plt.xlabel('Energy [keV]')
elif Plot.xAxis == 'channel':
        plt.xlabel('Channel')
plt.ylabel('Energy flux density [erg/s/cm$^2$/keV]')
print('saving plot...')
plt.legend()
plt.savefig(outputfiles_basename + 'unconvolved_posterior.pdf', bbox_inches='tight')
plt.close()



print('creating quantile-quantile plot ...')
plt.figure(figsize=(7,7))
with bxa.XSilence():
        solver.set_best_fit()
        bxa.qq.qq(prefix=outputfiles_basename, markers=5, annotate=True)
print('saving plot...')
plt.savefig(outputfiles_basename + 'qq_model_deviations.pdf', bbox_inches='tight')
plt.close()

1	"""
2	Example of doing BXA in X-spec
3	"""
4	import bxa.xspec as bxa	1✔
5	from xspec import *	1✔
6
7	Fit.statMethod = 'cstat'	1✔
8	Plot.xAxis = 'keV'	1✔
9
10	s = Spectrum('example-file.fak')	1✔
11	s.ignore("**"); s.notice("0.2-8.0")	1✔
12
13	# set model and parameters ranges
14	# val, delta, min, bottom, top, max
15	m = Model("pow")	1✔
16	m.powerlaw.norm.values = ",,1e-10,1e-10,1e1,1e1" # 10^-10 .. 10	1✔
17	m.powerlaw.PhoIndex.values = ",,1,1,3,3" # 1 .. 3	1✔
18
19	# define prior
20	transformations = [	1✔
21	# uniform prior for Photon Index (see other example for
22	# something more advanced)
23	bxa.create_uniform_prior_for( m, m.powerlaw.PhoIndex),
24	# jeffreys prior for scale variable
25	bxa.create_jeffreys_prior_for(m, m.powerlaw.norm),
26	# and possibly many more parameters here
27	]
28	print('running analysis ...')	1✔
29	# where to store intermediate and final results? this is the prefix used
30	outputfiles_basename = 'simplest/'	1✔
31	solver = bxa.BXASolver(transformations=transformations, outputfiles_basename=outputfiles_basename)	1✔
32	results = solver.run(resume=True)	1✔
33	print('running analysis ... done!')	1✔
34
35
36	import matplotlib.pyplot as plt	1✔
37
38	print('creating plot of posterior predictions against data ...')	1✔
39	plt.figure()	1✔
40
41	from xspec import Plot	1✔
42	from bxa.xspec.solver import set_parameters, XSilence	1✔
43	from ultranest.plot import PredictionBand	1✔
44
45	band = None	1✔
46	Plot.background = True	1✔
47
48	with XSilence():	1✔
49	Plot.device = '/null'	1✔
50	# plot models
51	for row in solver.posterior[:400]:	1✔
52	set_parameters(values=row, transformations=solver.transformations)	1✔
53	Plot('counts')	1✔
54	if band is None:	1✔
55	band = PredictionBand(Plot.x())	1✔
56	band.add(Plot.model())	1✔
57
58	band.shade(alpha=0.5, color='blue')	1✔
59	band.shade(q=0.495, alpha=0.1, color='blue')	1✔
60	band.line(color='blue', label='convolved model')	1✔
61
62	plt.scatter(Plot.x(), Plot.y(), label='data')	1✔
63	#plt.errorbar(
64	# x=Plot.x(), xerr=Plot.xErr(), y=Plot.y(), yerr=Plot.yErr(),
65	# marker='o', label='data')
66
67	if Plot.xAxis == 'keV':	1✔
68	plt.xlabel('Energy [keV]')	1✔
69	elif Plot.xAxis == 'channel':	×
70	plt.xlabel('Channel')	×
71	plt.ylabel('Counts/s/cm$^2$')	1✔
72	plt.savefig(outputfiles_basename + 'convolved_posterior_direct.pdf', bbox_inches='tight')	1✔
73	plt.close()	1✔
74
75
76
77	print('creating a rebinned plot to check the model ...')	1✔
78	plt.figure()	1✔
79	data = solver.posterior_predictions_convolved(nsamples=100)	1✔
80	print('binning for plot...')	1✔
81	binned = bxa.binning(outputfiles_basename=outputfiles_basename,	1✔
82	bins = data['bins'], widths = data['width'],
83	data = data['data'], models = data['models'])
84	for point in binned['marked_binned']:	1✔
85	plt.errorbar(marker='o', zorder=-1, **point)	1✔
86	plt.xlim(binned['xlim'])	1✔
87	plt.ylim(binned['ylim'][0], binned['ylim'][1]*2)	1✔
88	plt.gca().set_yscale('log')	1✔
89	if Plot.xAxis == 'keV':	1✔
90	plt.xlabel('Energy [keV]')	1✔
91	elif Plot.xAxis == 'channel':	×
92	plt.xlabel('Channel')	×
93	plt.ylabel('Counts/s/cm$^2$')	1✔
94	print('saving plot...')	1✔
95	plt.legend()	1✔
96	plt.savefig(outputfiles_basename + 'convolved_posterior.pdf', bbox_inches='tight')	1✔
97	plt.close()	1✔
98
99
100
101
102	print('creating plot of posterior predictions ...')	1✔
103	plt.figure()	1✔
104	solver.posterior_predictions_unconvolved(nsamples=100)	1✔
105	ylim = plt.ylim()	1✔
106	# 3 orders of magnitude at most
107	plt.ylim(max(ylim[0], ylim[1] / 1000), ylim[1])	1✔
108	plt.yscale('log')	1✔
109	plt.xscale('log')	1✔
110	if Plot.xAxis == 'keV':	1✔
111	plt.xlabel('Energy [keV]')	1✔
112	elif Plot.xAxis == 'channel':	×
113	plt.xlabel('Channel')	×
114	plt.ylabel('Energy flux density [erg/s/cm$^2$/keV]')	1✔
115	print('saving plot...')	1✔
116	plt.legend()	1✔
117	plt.savefig(outputfiles_basename + 'unconvolved_posterior.pdf', bbox_inches='tight')	1✔
118	plt.close()	1✔
119
120
121
122	print('creating quantile-quantile plot ...')	1✔
123	plt.figure(figsize=(7,7))	1✔
124	with bxa.XSilence():	1✔
125	solver.set_best_fit()	1✔
126	bxa.qq.qq(prefix=outputfiles_basename, markers=5, annotate=True)	1✔
127	print('saving plot...')	1✔
128	plt.savefig(outputfiles_basename + 'qq_model_deviations.pdf', bbox_inches='tight')	1✔
129	plt.close()	1✔

JohannesBuchner / BXA / 4754829731

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous