Example to use the sherpa plugin¶

import warnings
warnings.filterwarnings('ignore')

import matplotlib.pylab as plt
import jetset
from jetset.test_data_helper import  test_SEDs
from jetset.data_loader import ObsData,Data
from jetset.plot_sedfit import PlotSED
from jetset.test_data_helper import  test_SEDs

test_SEDs

['/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_3C345.ecsv',
 '/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk421_EBL_DEABS.ecsv',
 '/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_ABS.ecsv',
 '/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_DEABS.ecsv']

Loading data¶

see the data_format user guide for further information about loading data

print(test_SEDs[2])
data=Data.from_file(test_SEDs[2])

/Users/orion/anaconda3/envs/jetset/lib/python3.8/site-packages/jetset/test_data/SEDs_data/SED_MW_Mrk501_EBL_ABS.ecsv

%matplotlib inline
sed_data=ObsData(data_table=data)
sed_data.group_data(bin_width=0.2)

sed_data.add_systematics(0.1,[10.**6,10.**29])
p=sed_data.plot_sed()

================================================================================

*  binning data  *
---> N bins= 90
---> bin_widht= 0.2
================================================================================

sed_data.save('Mrk_501.pkl')

phenomenological model constraining¶

see the Phenomenological model constraining: application user guide for further information about phenomenological constraining

spectral indices¶

from jetset.sed_shaper import  SEDShape
my_shape=SEDShape(sed_data)
my_shape.eval_indices(minimizer='lsb',silent=True)
p=my_shape.plot_indices()
p.rescale(y_min=-15,y_max=-6)

================================================================================

* evaluating spectral indices for data *
================================================================================

sed shaper¶

mm,best_fit=my_shape.sync_fit(check_host_gal_template=True,
                  Ep_start=None,
                  minimizer='lsb',
                  silent=True,
                  fit_range=[10.,21.])

================================================================================

* Log-Polynomial fitting of the synchrotron component *
---> first blind fit run,  fit range: [10.0, 21.0]
---> class:  HSP

---> class:  HSP

Table length=6

model name	name	val	bestfit val	err +	err -	start val	fit range min	fit range max	frozen
LogCubic	b	-6.411144e-02	-6.411144e-02	7.838965e-03	--	-4.778764e-02	-1.000000e+01	0.000000e+00	False
LogCubic	c	-1.751721e-03	-1.751721e-03	1.127030e-03	--	3.576201e-03	-1.000000e+01	1.000000e+01	False
LogCubic	Ep	1.703747e+01	1.703747e+01	9.437354e-02	--	1.626870e+01	0.000000e+00	3.000000e+01	False
LogCubic	Sp	-1.030068e+01	-1.030068e+01	1.884114e-02	--	-1.025412e+01	-3.000000e+01	0.000000e+00	False
host_galaxy	nuFnu_p_host	-1.006557e+01	-1.006557e+01	5.462528e-02	--	-1.025412e+01	-1.225412e+01	-8.254123e+00	False
host_galaxy	nu_scale	1.730764e-02	1.730764e-02	3.694887e-03	--	0.000000e+00	-5.000000e-01	5.000000e-01	False

---> sync       nu_p=+1.703747e+01 (err=+9.437354e-02)  nuFnu_p=-1.030068e+01 (err=+1.884114e-02) curv.=-6.411144e-02 (err=+7.838965e-03)
================================================================================

my_shape.IC_fit(fit_range=[23.,29.],minimizer='minuit',silent=True)
p=my_shape.plot_shape_fit()
p.rescale(y_min=-15)

================================================================================

* Log-Polynomial fitting of the IC component *
---> fit range: [23.0, 29.0]
---> LogCubic fit

Table length=4

model name	name	val	bestfit val	err +	err -	start val	fit range min	fit range max	frozen
LogCubic	b	-1.310993e-01	-1.310993e-01	3.244188e-02	--	-1.000000e+00	-1.000000e+01	0.000000e+00	False
LogCubic	c	-3.300446e-02	-3.300446e-02	2.072521e-02	--	-1.000000e+00	-1.000000e+01	1.000000e+01	False
LogCubic	Ep	2.549603e+01	2.549603e+01	2.235473e-01	--	2.556357e+01	0.000000e+00	3.000000e+01	False
LogCubic	Sp	-1.057945e+01	-1.057945e+01	4.332978e-02	--	-1.000000e+01	-3.000000e+01	0.000000e+00	False

---> IC         nu_p=+2.549603e+01 (err=+2.235473e-01)  nuFnu_p=-1.057945e+01 (err=+4.332978e-02) curv.=-1.310993e-01 (err=+3.244188e-02)
================================================================================

Model constraining¶

In this step we are not fitting the model, we are just obtaining the phenomenological pre_fit model, that will be fitted in using minuit ore least-square bound, as shown below

from jetset.obs_constrain import ObsConstrain
from jetset.model_manager import  FitModel
sed_obspar=ObsConstrain(beaming=25,
                        B_range=[0.001,0.1],
                        distr_e='lppl',
                        t_var_sec=3*86400,
                        nu_cut_IR=1E12,
                        SEDShape=my_shape)


prefit_jet=sed_obspar.constrain_SSC_model(electron_distribution_log_values=False,silent=True)
prefit_jet.save_model('prefit_jet.pkl')

================================================================================

*  constrains parameters from observable *

Table length=11

model name	name	par type	units	val	phys. bound. min	phys. bound. max	log	frozen
jet_leptonic	R	region_size	cm	1.046299e+16	1.000000e+03	1.000000e+30	False	False
jet_leptonic	R_H	region_position	cm	1.000000e+17	0.000000e+00	--	False	True
jet_leptonic	B	magnetic_field	gauss	5.050000e-02	0.000000e+00	--	False	False
jet_leptonic	beam_obj	beaming	lorentz-factor*	2.500000e+01	1.000000e-04	--	False	False
jet_leptonic	z_cosm	redshift		3.360000e-02	0.000000e+00	--	False	False
jet_leptonic	gmin	low-energy-cut-off	lorentz-factor*	4.703917e+02	1.000000e+00	1.000000e+09	False	False
jet_leptonic	gmax	high-energy-cut-off	lorentz-factor*	2.310708e+06	1.000000e+00	1.000000e+15	False	False
jet_leptonic	N	emitters_density	1 / cm3	7.305900e+00	0.000000e+00	--	False	False
jet_leptonic	gamma0_log_parab	turn-over-energy	lorentz-factor*	1.045843e+04	1.000000e+00	1.000000e+09	False	False
jet_leptonic	s	LE_spectral_slope		2.248787e+00	-1.000000e+01	1.000000e+01	False	False
jet_leptonic	r	spectral_curvature		3.205572e-01	-1.500000e+01	1.500000e+01	False	False

================================================================================

pl=prefit_jet.plot_model(sed_data=sed_data)
pl.add_model_residual_plot(prefit_jet,sed_data)
pl.rescale(y_min=-15,x_min=7,x_max=29)

Model fitting with Sherpa¶

from jetset.sherpa_plugin import JetsetSherpaModel

from jetset.template_2Dmodel import EBLAbsorptionTemplate
ebl_franceschini=EBLAbsorptionTemplate.from_name('Franceschini_2008')

from jetset.jet_model import Jet
prefit_jet=Jet.load_model('prefit_jet.pkl')

Table length=11

model name	name	par type	units	val	phys. bound. min	phys. bound. max	log	frozen
jet_leptonic	gmin	low-energy-cut-off	lorentz-factor*	4.703917e+02	1.000000e+00	1.000000e+09	False	False
jet_leptonic	gmax	high-energy-cut-off	lorentz-factor*	2.310708e+06	1.000000e+00	1.000000e+15	False	False
jet_leptonic	N	emitters_density	1 / cm3	7.305900e+00	0.000000e+00	--	False	False
jet_leptonic	gamma0_log_parab	turn-over-energy	lorentz-factor*	1.045843e+04	1.000000e+00	1.000000e+09	False	False
jet_leptonic	s	LE_spectral_slope		2.248787e+00	-1.000000e+01	1.000000e+01	False	False
jet_leptonic	r	spectral_curvature		3.205572e-01	-1.500000e+01	1.500000e+01	False	False
jet_leptonic	R	region_size	cm	1.046299e+16	1.000000e+03	1.000000e+30	False	False
jet_leptonic	R_H	region_position	cm	1.000000e+17	0.000000e+00	--	False	True
jet_leptonic	B	magnetic_field	gauss	5.050000e-02	0.000000e+00	--	False	False
jet_leptonic	beam_obj	beaming	lorentz-factor*	2.500000e+01	1.000000e-04	--	False	False
jet_leptonic	z_cosm	redshift		3.360000e-02	0.000000e+00	--	False	False

sherpa_model_jet=JetsetSherpaModel(prefit_jet)
sherpa_model_gal=JetsetSherpaModel(my_shape.host_gal)
sherpa_model_ebl=JetsetSherpaModel(ebl_franceschini)

jetset model name R renamed to  R_sh due to sherpa internal naming convention

sherpa_model=(sherpa_model_jet+sherpa_model_gal)*sherpa_model_ebl

sherpa_model

Model

Expression: (jet_leptonic + host_galaxy) * Franceschini_2008
Component	Parameter	Value	Min	Max	Units
jet_leptonic	gmin	470.39174855643597	1.0	1000000000.0	lorentz-factor*
	gmax	2310708.197406515	1.0	1000000000000000.0	lorentz-factor*
	N	7.305900046565794	0.0	MAX	1 / cm3
	gamma0_log_parab	10458.433237258416	1.0	1000000000.0	lorentz-factor*
	s	2.2487867709713574	-10.0	10.0
	r	0.32055721323507314	-15.0	15.0
	R_sh	1.046298751763524e+16	1000.0	1e+30	cm
	R_H	1e+17	0.0	MAX	cm
	B	0.0505	0.0	MAX	gauss
	beam_obj	25.0	0.0001	MAX	lorentz-factor*
	z_cosm	0.0336	0.0	MAX
host_galaxy	nuFnu_p_host	-10.065573840795157	-20.0	20.0	erg / (cm2 s)
host_galaxy	nu_scale	0.01730764400081941	-2.0	2.0	Hz
Franceschini_2008	z_cosm	1.0	0.0	MAX

sherpa_model_ebl.z_cosm  = sherpa_model_jet.z_cosm

sherpa_model

Model

Expression: (jet_leptonic + host_galaxy) * Franceschini_2008
Component	Parameter	Thawed	Value	Min	Max	Units
jet_leptonic	gmin		470.39174855643597	1.0	1000000000.0	lorentz-factor*
	gmax		2310708.197406515	1.0	1000000000000000.0	lorentz-factor*
	N		7.305900046565794	0.0	MAX	1 / cm3
	gamma0_log_parab		10458.433237258416	1.0	1000000000.0	lorentz-factor*
	s		2.2487867709713574	-10.0	10.0
	r		0.32055721323507314	-15.0	15.0
	R_sh		1.046298751763524e+16	1000.0	1e+30	cm
	R_H		1e+17	0.0	MAX	cm
	B		0.0505	0.0	MAX	gauss
	beam_obj		25.0	0.0001	MAX	lorentz-factor*
	z_cosm		0.0336	0.0	MAX
host_galaxy	nuFnu_p_host		-10.065573840795157	-20.0	20.0	erg / (cm2 s)
host_galaxy	nu_scale		0.01730764400081941	-2.0	2.0	Hz
Franceschini_2008	z_cosm	linked	0.0336	⇐ jet_leptonic.z_cosm

sherpa_model_jet.R_H.freeze()
sherpa_model_jet.z_cosm.freeze()

sherpa_model

Model

Expression: (jet_leptonic + host_galaxy) * Franceschini_2008
Component	Parameter	Thawed	Value	Min	Max	Units
jet_leptonic	gmin		470.39174855643597	1.0	1000000000.0	lorentz-factor*
	gmax		2310708.197406515	1.0	1000000000000000.0	lorentz-factor*
	N		7.305900046565794	0.0	MAX	1 / cm3
	gamma0_log_parab		10458.433237258416	1.0	1000000000.0	lorentz-factor*
	s		2.2487867709713574	-10.0	10.0
	r		0.32055721323507314	-15.0	15.0
	R_sh		1.046298751763524e+16	1000.0	1e+30	cm
	R_H		1e+17	0.0	MAX	cm
	B		0.0505	0.0	MAX	gauss
	beam_obj		25.0	0.0001	MAX	lorentz-factor*
	z_cosm		0.0336	0.0	MAX
host_galaxy	nuFnu_p_host		-10.065573840795157	-20.0	20.0	erg / (cm2 s)
host_galaxy	nu_scale		0.01730764400081941	-2.0	2.0	Hz
Franceschini_2008	z_cosm	linked	0.0336	⇐ jet_leptonic.z_cosm

from sherpa import data
from sherpa.fit import Fit
from sherpa.stats import Chi2
from sherpa.optmethods import LevMar, NelderMead

sherpa_data=data.Data1D("sed", sed_data.table['nu_data'], sed_data.table['nuFnu_data'], staterror=sed_data.table['dnuFnu_data'])

fitter = Fit(sherpa_data, sherpa_model, stat=Chi2(), method=LevMar())
fit_range=[1e11,1e29]

sherpa_data.notice(fit_range[0], fit_range[1])

results = fitter.fit()
print("-- fit succesful?", results.succeeded)
print(results.format())

-- fit succesful? True
Method                = levmar
Statistic             = chi2
Initial fit statistic = 159.009
Final fit statistic   = 12.8563 at function evaluation 379
Data points           = 31
Degrees of freedom    = 20
Probability [Q-value] = 0.883462
Reduced statistic     = 0.642817
Change in statistic   = 146.152
   jet_leptonic.gmin   161.587      +/- 537.608
   jet_leptonic.gmax   2.6044e+06   +/- 0
   jet_leptonic.N   12.4908      +/- 49.9373
   jet_leptonic.gamma0_log_parab   12476.6      +/- 0
   jet_leptonic.s   2.22409      +/- 0.0334469
   jet_leptonic.r   0.246937     +/- 0.0213915
   jet_leptonic.R_sh   1.69738e+16  +/- 0
   jet_leptonic.B   0.00724202   +/- 0.00193696
   jet_leptonic.beam_obj   48.2305      +/- 6.91246
   host_galaxy.nuFnu_p_host   -10.0495     +/- 0.0846434
   host_galaxy.nu_scale   0.0235454    +/- 0.00763641

sherpa_model

Model

Expression: (jet_leptonic + host_galaxy) * Franceschini_2008
Component	Parameter	Thawed	Value	Min	Max	Units
jet_leptonic	gmin		161.58718772367948	1.0	1000000000.0	lorentz-factor*
	gmax		2604402.6487550396	1.0	1000000000000000.0	lorentz-factor*
	N		12.490793758971192	0.0	MAX	1 / cm3
	gamma0_log_parab		12476.58647443473	1.0	1000000000.0	lorentz-factor*
	s		2.2240928430873117	-10.0	10.0
	r		0.24693749268885934	-15.0	15.0
	R_sh		1.6973782230807432e+16	1000.0	1e+30	cm
	R_H		1e+17	0.0	MAX	cm
	B		0.00724201775285936	0.0	MAX	gauss
	beam_obj		48.23048000534593	0.0001	MAX	lorentz-factor*
	z_cosm		0.0336	0.0	MAX
host_galaxy	nuFnu_p_host		-10.04953502514069	-20.0	20.0	erg / (cm2 s)
host_galaxy	nu_scale		0.02354535123290668	-2.0	2.0	Hz
Franceschini_2008	z_cosm	linked	0.0336	⇐ jet_leptonic.z_cosm

from jetset.sherpa_plugin import plot_sherpa_model

def plot_sherpa_model(sherpa_model, fit_range=None, model_range=[1E10, 1E30], nu_grid_size=200, sed_data=None,add_res=False,plot_obj=None,label=None,line_style=None):

    x = np.logspace(np.log10(model_range[0]), np.log10(model_range[1]), nu_grid_size)
    y = sherpa_model(x)
    if plot_obj is None:
        plot_obj = PlotSED(sed_data=sed_data, frame='obs', density=False)

    plot_obj.add_xy_plot(np.log10(x), np.log10(y),label=label,line_style=line_style)

    if add_res is True and fit_range is not None:
        nufnu_res = sherpa_model(sed_data.data['nu_data'])
        y_res = (sed_data.data['nuFnu_data'] - nufnu_res) / sed_data.data['dnuFnu_data']
        x_res = sed_data.data['nu_data']
        plot_obj.add_xy_residual_plot(x=np.log10(x_res), y=y_res, fit_range=np.log10([fit_range[0], fit_range[1]]))

    return  plot_obj

p=plot_sherpa_model(sherpa_model=sherpa_model,sed_data=sed_data,fit_range=fit_range,add_res=True,label='(SSC+host gal)*ebl')
p=plot_sherpa_model(sherpa_model_jet,plot_obj=p,label='SSC',line_style='--')
p=plot_sherpa_model(sherpa_model_gal,plot_obj=p,label='host gal',line_style='--')

from sherpa.plot import IntervalProjection
iproj = IntervalProjection()
iproj.prepare(fac=5, nloop=15)
iproj.calc(fitter, sherpa_model.s)
iproj.plot()

from sherpa.sim import MCMC

mcmc = MCMC()

mcmc.get_sampler_name()

'MetropolisMH'

fitter.estmethod

<Covariance error-estimation method instance 'covariance'>

eres = fitter.est_errors()

WARNING: hard minimum hit for parameter jet_leptonic.gmax
WARNING: hard maximum hit for parameter jet_leptonic.gmax
WARNING: hard minimum hit for parameter jet_leptonic.beam_obj
WARNING: hard maximum hit for parameter jet_leptonic.beam_obj
WARNING: hard minimum hit for parameter host_galaxy.nuFnu_p_host
WARNING: hard maximum hit for parameter host_galaxy.nuFnu_p_host
WARNING: hard minimum hit for parameter host_galaxy.nu_scale
WARNING: hard maximum hit for parameter host_galaxy.nu_scale

print(eres.format())

Confidence Method     = covariance
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = chi2
covariance 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   jet_leptonic.gmin      161.587     -17.2539      17.2539
   jet_leptonic.gmax   2.6044e+06        -----        -----
   jet_leptonic.N      12.4908     -2.43741      2.43741
   jet_leptonic.gamma0_log_parab      12476.6     -1835.53      1835.53
   jet_leptonic.s      2.22409    -0.018924     0.018924
   jet_leptonic.r     0.246937   -0.0231805    0.0231805
   jet_leptonic.R_sh  1.69738e+16 -1.54339e+15  1.54339e+15
   jet_leptonic.B   0.00724202 -0.000565301  0.000565301
   jet_leptonic.beam_obj      48.2305        -----        -----
   host_galaxy.nuFnu_p_host     -10.0495        -----        -----
   host_galaxy.nu_scale    0.0235454        -----        -----

cmatrix = eres.extra_output
pnames = [p.split('.')[1] for p in eres.parnames]
plt.imshow(cmatrix, interpolation='nearest', cmap='viridis')
plt.xticks(np.arange(len(pnames)), pnames)
plt.yticks(np.arange(len(pnames)), pnames)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7f8d86c26700>

draws = mcmc.get_draws(fitter, cmatrix, niter=200)

---------------------------------------------------------------------------

ValueError                                Traceback (most recent call last)

<ipython-input-40-e4874d746f87> in <module>
----> 1 draws = mcmc.get_draws(fitter, cmatrix, niter=200)


~/anaconda3/envs/jetset/lib/python3.8/site-packages/sherpa/sim/__init__.py in get_draws(self, fit, sigma, niter, cache)
    685         """
    686         if not isinstance(fit.stat, (Cash, CStat, WStat)):
--> 687             raise ValueError("Fit statistic must be cash, cstat or " +
    688                              "wstat, not %s" % fit.stat.name)
    689


ValueError: Fit statistic must be cash, cstat or wstat, not chi2