Depending parameters#
In the following we show how to link parameters in the same model or among different models, and how to make a paramter depending on other parameters according to a mathematical expression.
Example: linked paramters for EBL#
import warnings
warnings.filterwarnings('ignore')
import jetset
print('tested with',jetset.__version__)
tested on with 1.4.0rc0
from jetset.jet_model import Jet
from jetset.template_2Dmodel import EBLAbsorptionTemplate
from jetset.model_manager import FitModel
my_jet = Jet(electron_distribution='lppl', name='jet_flaring')
my_jet.parameters.z_cosm.val = 0.01
ebl_franceschini = EBLAbsorptionTemplate.from_name('Franceschini_2008')
composite_model = FitModel(nu_size=500, name='EBL corrected')
composite_model.add_component(my_jet)
composite_model.add_component(ebl_franceschini)
composite_model.show_pars()
composite_model.link_par(par_name='z_cosm', from_model='Franceschini_2008', to_model='jet_flaring')
v=0.03001
my_jet.parameters.z_cosm.val = v
assert (composite_model.Franceschini_2008.parameters.z_cosm.val==v)
assert (composite_model.Franceschini_2008.parameters.z_cosm.linked==True)
composite_model.composite_expr = '%s*%s'%(my_jet.name,ebl_franceschini.name)
composite_model.eval()
#if plot is True:
# composite_model.plot_model()
composite_model.save_model('ebl_jet.pkl')
new_composite_model=FitModel.load_model('ebl_jet.pkl')
new_composite_model.show_pars()
v=2.0
new_composite_model.jet_flaring.parameters.z_cosm.val=v
print('new_composite_model.Franceschini_2008.parameters.z_cosm.val',new_composite_model.Franceschini_2008.parameters.z_cosm.val,'v',v)
assert (new_composite_model.Franceschini_2008.parameters.z_cosm.val == v)
assert (new_composite_model.Franceschini_2008.parameters.z_cosm.linked == True)
| model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
|---|---|---|---|---|---|---|---|---|
| jet_flaring | R | region_size | cm | 5.000000e+15 | 1.000000e+03 | 1.000000e+30 | False | False |
| jet_flaring | R_H | region_position | cm | 1.000000e+17 | 0.000000e+00 | -- | False | True |
| jet_flaring | B | magnetic_field | gauss | 1.000000e-01 | 0.000000e+00 | -- | False | False |
| jet_flaring | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_flaring | beam_obj | beaming | 1.000000e+01 | 1.000000e-04 | -- | False | False | |
| jet_flaring | z_cosm | redshift | 1.000000e-02 | 0.000000e+00 | -- | False | False | |
| jet_flaring | gmin | low-energy-cut-off | lorentz-factor* | 2.000000e+00 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_flaring | gmax | high-energy-cut-off | lorentz-factor* | 1.000000e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_flaring | N | emitters_density | 1 / cm3 | 1.000000e+02 | 0.000000e+00 | -- | False | False |
| jet_flaring | gamma0_log_parab | turn-over-energy | lorentz-factor* | 1.000000e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_flaring | s | LE_spectral_slope | 2.000000e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_flaring | r | spectral_curvature | 4.000000e-01 | -1.500000e+01 | 1.500000e+01 | False | False | |
| Franceschini_2008 | scale_factor | scale_factor | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| Franceschini_2008 | z_cosm | redshift | 1.000000e+00 | 0.000000e+00 | -- | False | True |
| model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
|---|---|---|---|---|---|---|---|---|
| jet_flaring | gmin | low-energy-cut-off | lorentz-factor* | 2.000000e+00 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_flaring | gmax | high-energy-cut-off | lorentz-factor* | 1.000000e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_flaring | N | emitters_density | 1 / cm3 | 1.000000e+02 | 0.000000e+00 | -- | False | False |
| jet_flaring | gamma0_log_parab | turn-over-energy | lorentz-factor* | 1.000000e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_flaring | s | LE_spectral_slope | 2.000000e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_flaring | r | spectral_curvature | 4.000000e-01 | -1.500000e+01 | 1.500000e+01 | False | False | |
| jet_flaring | R | region_size | cm | 5.000000e+15 | 1.000000e+03 | 1.000000e+30 | False | False |
| jet_flaring | R_H | region_position | cm | 1.000000e+17 | 0.000000e+00 | -- | False | True |
| jet_flaring | B | magnetic_field | gauss | 1.000000e-01 | 0.000000e+00 | -- | False | False |
| jet_flaring | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_flaring | beam_obj | beaming | 1.000000e+01 | 1.000000e-04 | -- | False | False | |
| jet_flaring | z_cosm(M) | redshift | 3.001000e-02 | 0.000000e+00 | -- | False | False | |
| Franceschini_2008 | scale_factor | scale_factor | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| Franceschini_2008 | z_cosm(L,jet_flaring) | redshift | -- | -- | -- | False | True |
new_composite_model.Franceschini_2008.parameters.z_cosm.val 2.0 v 2.0
Example: depending pars for bkn power-law emitters#
here we create a custom bkn distribution where we impose a
functional dependence among the low and high-energy spectral index.
from jetset.jet_emitters import EmittersDistribution
import numpy as np
from jetset.jet_model import Jet
j = Jet(emitters_distribution='bkn')
j.parameters
| model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
|---|---|---|---|---|---|---|---|---|
| jet_leptonic | R | region_size | cm | 5.000000e+15 | 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 | 1.000000e-01 | 0.000000e+00 | -- | False | False |
| jet_leptonic | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_leptonic | beam_obj | beaming | 1.000000e+01 | 1.000000e-04 | -- | False | False | |
| jet_leptonic | z_cosm | redshift | 1.000000e-01 | 0.000000e+00 | -- | False | False | |
| jet_leptonic | gmin | low-energy-cut-off | lorentz-factor* | 2.000000e+00 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 1.000000e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_leptonic | N | emitters_density | 1 / cm3 | 1.000000e+02 | 0.000000e+00 | -- | False | False |
| jet_leptonic | gamma_break | turn-over-energy | lorentz-factor* | 1.000000e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | p | LE_spectral_slope | 2.500000e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_leptonic | p_1 | HE_spectral_slope | 3.500000e+00 | -1.000000e+01 | 1.000000e+01 | False | False |
None
the functional dependence can be provided by a python function, where
the argument (p in this case) is the same name as the parameter:
def f_p(p):
return p+1
j.make_dependent_par(par='p_1',depends_on=['p'],par_expr=f_p)
j.parameters.p.val=2
np.testing.assert_allclose(j.parameters.p_1.val, j.parameters.p.val + 1)
j.parameters
adding par: p to p_1
==> par p_1 is depending on ['p'] according to expr: p_1 =
def f_p(p):
return p+1
| model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
|---|---|---|---|---|---|---|---|---|
| jet_leptonic | R | region_size | cm | 5.000000e+15 | 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 | 1.000000e-01 | 0.000000e+00 | -- | False | False |
| jet_leptonic | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_leptonic | beam_obj | beaming | 1.000000e+01 | 1.000000e-04 | -- | False | False | |
| jet_leptonic | z_cosm | redshift | 1.000000e-01 | 0.000000e+00 | -- | False | False | |
| jet_leptonic | gmin | low-energy-cut-off | lorentz-factor* | 2.000000e+00 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 1.000000e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_leptonic | N | emitters_density | 1 / cm3 | 1.000000e+02 | 0.000000e+00 | -- | False | False |
| jet_leptonic | gamma_break | turn-over-energy | lorentz-factor* | 1.000000e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | p(M) | LE_spectral_slope | 2.000000e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_leptonic | *p_1(D,p) | HE_spectral_slope | 3.000000e+00 | -1.000000e+01 | 1.000000e+01 | False | True |
None
as you can notice, now a message is shown describing the dependence of the parameters
It is also possible to set the dependence function as a string that can be evaluated
j.make_dependent_par(par='p_1',depends_on=['p'],par_expr='p+1')
j.parameters.p.val=2
np.testing.assert_allclose(j.parameters.p_1.val, j.parameters.p.val + 1)
j.parameters
==> par p_1 is depending on ['p'] according to expr: p_1 =
p+1
| model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
|---|---|---|---|---|---|---|---|---|
| jet_leptonic | R | region_size | cm | 5.000000e+15 | 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 | 1.000000e-01 | 0.000000e+00 | -- | False | False |
| jet_leptonic | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_leptonic | beam_obj | beaming | 1.000000e+01 | 1.000000e-04 | -- | False | False | |
| jet_leptonic | z_cosm | redshift | 1.000000e-01 | 0.000000e+00 | -- | False | False | |
| jet_leptonic | gmin | low-energy-cut-off | lorentz-factor* | 2.000000e+00 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 1.000000e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_leptonic | N | emitters_density | 1 / cm3 | 1.000000e+02 | 0.000000e+00 | -- | False | False |
| jet_leptonic | gamma_break | turn-over-energy | lorentz-factor* | 1.000000e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | p(M) | LE_spectral_slope | 2.000000e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_leptonic | *p_1(D,p) | HE_spectral_slope | 3.000000e+00 | -1.000000e+01 | 1.000000e+01 | False | True |
None
In principle, you can use strings for short expressions, and functions for more complicated formulas.
You can print the actual expression/function for the depending parameter
using the print_par_expr method:
j.parameters.p_1.par_expression_source_code
==> par p_1 is depending on ['p'] according to expr: p_1 =
p+1
j.save_model('jet.pkl')
new_jet=Jet.load_model('jet.pkl')
new_jet.parameters.p.val=2.5
np.testing.assert_allclose(new_jet.parameters.p_1.val, new_jet.parameters.p.val + 1)
new_jet.parameters
| model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
|---|---|---|---|---|---|---|---|---|
| jet_leptonic | gmin | low-energy-cut-off | lorentz-factor* | 2.000000e+00 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 1.000000e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_leptonic | N | emitters_density | 1 / cm3 | 1.000000e+02 | 0.000000e+00 | -- | False | False |
| jet_leptonic | gamma_break | turn-over-energy | lorentz-factor* | 1.000000e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | p(M) | LE_spectral_slope | 2.500000e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_leptonic | *p_1(D,p) | HE_spectral_slope | 3.500000e+00 | -1.000000e+01 | 1.000000e+01 | False | True | |
| jet_leptonic | R | region_size | cm | 5.000000e+15 | 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 | 1.000000e-01 | 0.000000e+00 | -- | False | False |
| jet_leptonic | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_leptonic | beam_obj | beaming | 1.000000e+01 | 1.000000e-04 | -- | False | False | |
| jet_leptonic | z_cosm | redshift | 1.000000e-01 | 0.000000e+00 | -- | False | False |
None
Example depending par: Building a Jet model with B function of R_H and R_0#
In this example we create a fuctional dependence among the paramters
B, R_H introducing user custom pararameters. Wewant that the
value of the mangentic field in the jet is a function or R_H, and of
the initial value of B=B0 at R=R_H0, according to the
expression:
\(B=B_0(R_0/R_H)^{1.1}\)
jet=Jet(emitters_distribution='plc')
fit_model_lsb=FitModel( jet=jet, name='SSC-best-fit-lsb',template=None)
fit_model_lsb.jet_leptonic.parameters.beam_obj.fit_range = [5, 50]
fit_model_lsb.jet_leptonic.parameters.R_H.val=5E17
fit_model_lsb.jet_leptonic.parameters.R_H.frozen=False
fit_model_lsb.jet_leptonic.parameters.R_H.fit_range = [1E15, 1E19]
fit_model_lsb.jet_leptonic.parameters.R.fit_range = [10 ** 15.5, 10 ** 17.5]
fit_model_lsb.jet_leptonic.add_user_par(name='B0',units='G',val=1E3,val_min=0,val_max=None)
fit_model_lsb.jet_leptonic.add_user_par(name='R0', units='cm', val=5E13, val_min=0, val_max=None)
fit_model_lsb.jet_leptonic.add_user_par(name='m_B', val=1, val_min=1, val_max=2)
fit_model_lsb.jet_leptonic.parameters.R0.frozen=True
fit_model_lsb.jet_leptonic.parameters.B0.frozen=True
def par_func(R0,B0,R_H,m_B):
return B0*np.power((R0/R_H),m_B)
fit_model_lsb.jet_leptonic.make_dependent_par(par='B', depends_on=['B0', 'R0', 'R_H','m_B'], par_expr=par_func)
B0=fit_model_lsb.jet_leptonic.parameters.B0.val
R0 = fit_model_lsb.jet_leptonic.parameters.R0.val
R_H = fit_model_lsb.jet_leptonic.parameters.R_H.val
m_B= fit_model_lsb.jet_leptonic.parameters.m_B.val
np.testing.assert_allclose(fit_model_lsb.jet_leptonic.parameters.B.val, par_func(R0,B0,R_H,m_B))
adding par: B0 to B
adding par: R0 to B
adding par: R_H to B
adding par: m_B to B
==> par B is depending on ['B0', 'R0', 'R_H', 'm_B'] according to expr: B =
def par_func(R0,B0,R_H,m_B):
return B0*np.power((R0/R_H),m_B)
fit_model_lsb.jet_leptonic.parameters
| model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
|---|---|---|---|---|---|---|---|---|
| jet_leptonic | R | region_size | cm | 5.000000e+15 | 1.000000e+03 | 1.000000e+30 | False | False |
| jet_leptonic | R_H(M) | region_position | cm | 5.000000e+17 | 0.000000e+00 | -- | False | False |
| jet_leptonic | *B(D,m_B) | magnetic_field | gauss | 1.000000e-01 | 0.000000e+00 | -- | False | True |
| jet_leptonic | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_leptonic | beam_obj | beaming | 1.000000e+01 | 1.000000e-04 | -- | False | False | |
| jet_leptonic | z_cosm | redshift | 1.000000e-01 | 0.000000e+00 | -- | False | False | |
| jet_leptonic | gmin | low-energy-cut-off | lorentz-factor* | 2.000000e+00 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 1.000000e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_leptonic | N | emitters_density | 1 / cm3 | 1.000000e+02 | 0.000000e+00 | -- | False | False |
| jet_leptonic | gamma_cut | turn-over-energy | lorentz-factor* | 1.000000e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | p | LE_spectral_slope | 2.000000e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_leptonic | B0(M) | user_defined | G | 1.000000e+03 | 0.000000e+00 | -- | False | True |
| jet_leptonic | R0(M) | user_defined | cm | 5.000000e+13 | 0.000000e+00 | -- | False | True |
| jet_leptonic | m_B(M) | user_defined | 1.000000e+00 | 1.000000e+00 | 2.000000e+00 | False | False |
None
fit_model_lsb.save_model('test.pkl')
fit_model_lsb=FitModel.load_model('test.pkl')
B0=fit_model_lsb.jet_leptonic.parameters.B0.val
R0 = fit_model_lsb.jet_leptonic.parameters.R0.val
R_H = fit_model_lsb.jet_leptonic.parameters.R_H.val
m_B= fit_model_lsb.jet_leptonic.parameters.m_B.val
np.testing.assert_allclose(fit_model_lsb.jet_leptonic.parameters.B.val, par_func(R0,B0,R_H,m_B))
%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(dpi=150)
R_H_array=np.logspace(13,18,100)
B_array=np.zeros(R_H_array.shape)
for ID,R_H in enumerate(R_H_array):
fit_model_lsb.jet_leptonic.parameters.R_H.val=R_H
B_array[ID]=fit_model_lsb.jet_leptonic.parameters.B.val
plt.loglog(R_H_array,B_array)
plt.xlabel('R_H (cm)')
plt.ylabel('B (G)')
Text(0, 0.5, 'B (G)')
Removing the dependenencies#
for the entire fit_model
fit_model_lsb.parameters.reset_dependencies()
or, for a specific component
fit_model_lsb.jet_leptonic.parameters.reset_dependencies()
fit_model_lsb.parameters
| model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
|---|---|---|---|---|---|---|---|---|
| jet_leptonic | gmin | low-energy-cut-off | lorentz-factor* | 2.000000e+00 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 1.000000e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_leptonic | N | emitters_density | 1 / cm3 | 1.000000e+02 | 0.000000e+00 | -- | False | False |
| jet_leptonic | gamma_cut | turn-over-energy | lorentz-factor* | 1.000000e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | p | LE_spectral_slope | 2.000000e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_leptonic | R | region_size | cm | 5.000000e+15 | 1.000000e+03 | 1.000000e+30 | False | False |
| jet_leptonic | R_H | region_position | cm | 1.000000e+18 | 0.000000e+00 | -- | False | False |
| jet_leptonic | B | magnetic_field | gauss | 5.000000e-02 | 0.000000e+00 | -- | False | True |
| jet_leptonic | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_leptonic | beam_obj | beaming | 1.000000e+01 | 1.000000e-04 | -- | False | False | |
| jet_leptonic | z_cosm | redshift | 1.000000e-01 | 0.000000e+00 | -- | False | False | |
| jet_leptonic | B0 | user_defined | G | 1.000000e+03 | 0.000000e+00 | -- | False | True |
| jet_leptonic | R0 | user_defined | cm | 5.000000e+13 | 0.000000e+00 | -- | False | True |
| jet_leptonic | m_B | user_defined | 1.000000e+00 | 1.000000e+00 | 2.000000e+00 | False | False |
None
%matplotlib inline
plt.figure(dpi=150)
R_H_array=np.logspace(13,18,100)
B_array=np.zeros(R_H_array.shape)
for ID,R_H in enumerate(R_H_array):
fit_model_lsb.jet_leptonic.parameters.R_H.val=R_H
B_array[ID]=fit_model_lsb.jet_leptonic.parameters.B.val
plt.loglog(R_H_array,B_array)
plt.xlabel('R_H (cm)')
plt.ylabel('B (G)')
Text(0, 0.5, 'B (G)')
Example depending par: fitting with a Jet model with depending pars#
In this example we show how to use the previous model during a Fit
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
data=Data.from_file(test_SEDs[1])
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= 88 ---> bin_width= 0.2 ================================================================================
from jetset.sed_shaper import SEDShape
my_shape=SEDShape(sed_data)
my_shape.eval_indices(minimizer='lsb',silent=True)
p=my_shape.plot_indices()
================================================================================ * evaluating spectral indices for data * ================================================================================
mm,best_fit=my_shape.sync_fit(check_host_gal_template=False,
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: HSPTable length=4
| model name | name | val | bestfit val | err + | err - | start val | fit range min | fit range max | frozen |
|---|---|---|---|---|---|---|---|---|---|
| LogCubic | b | -1.563747e-01 | -1.563747e-01 | 5.975434e-03 | -- | -1.000000e+00 | -1.000000e+01 | 0.000000e+00 | False |
| LogCubic | c | -1.052802e-02 | -1.052802e-02 | 8.781942e-04 | -- | -1.000000e+00 | -1.000000e+01 | 1.000000e+01 | False |
| LogCubic | Ep | 1.675324e+01 | 1.675324e+01 | 2.396636e-02 | -- | 1.670206e+01 | 0.000000e+00 | 3.000000e+01 | False |
| LogCubic | Sp | -9.494365e+00 | -9.494365e+00 | 1.704982e-02 | -- | -1.000000e+01 | -3.000000e+01 | 0.000000e+00 | False |
---> sync nu_p=+1.675324e+01 (err=+2.396636e-02) nuFnu_p=-9.494365e+00 (err=+1.704982e-02) curv.=-1.563747e-01 (err=+5.975434e-03)
================================================================================
my_shape.IC_fit(fit_range=[23.,29.],minimizer='minuit',silent=True)
p=my_shape.plot_shape_fit()
p.setlim(y_min=1E-15)
================================================================================ * Log-Polynomial fitting of the IC component * ---> fit range: [23.0, 29.0] ---> LogCubic fitTable length=4
| model name | name | val | bestfit val | err + | err - | start val | fit range min | fit range max | frozen |
|---|---|---|---|---|---|---|---|---|---|
| LogCubic | b | -2.274590e-01 | -2.274590e-01 | 3.262165e-02 | -- | -1.000000e+00 | -1.000000e+01 | 0.000000e+00 | False |
| LogCubic | c | -6.259967e-02 | -6.259967e-02 | 1.629407e-02 | -- | -1.000000e+00 | -1.000000e+01 | 1.000000e+01 | False |
| LogCubic | Ep | 2.527207e+01 | 2.527207e+01 | 8.149544e-02 | -- | 2.528644e+01 | 0.000000e+00 | 3.000000e+01 | False |
| LogCubic | Sp | -1.014119e+01 | -1.014119e+01 | 2.734754e-02 | -- | -1.000000e+01 | -3.000000e+01 | 0.000000e+00 | False |
---> IC nu_p=+2.527207e+01 (err=+8.149544e-02) nuFnu_p=-1.014119e+01 (err=+2.734754e-02) curv.=-2.274590e-01 (err=+3.262165e-02)
================================================================================
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 *
/Users/orion/miniforge3/envs/jetset/lib/python3.12/site-packages/jetset/obs_constrain.py:1514: RankWarning: Polyfit may be poorly conditioned
p=polyfit(nu_p_IC_model_log,B_grid_log,2)
| model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
|---|---|---|---|---|---|---|---|---|
| jet_leptonic | R | region_size | cm | 3.452668e+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 | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_leptonic | beam_obj | beaming | 2.500000e+01 | 1.000000e-04 | -- | False | False | |
| jet_leptonic | z_cosm | redshift | 3.080000e-02 | 0.000000e+00 | -- | False | False | |
| jet_leptonic | gmin | low-energy-cut-off | lorentz-factor* | 4.697542e+02 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 1.300733e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_leptonic | N | emitters_density | 1 / cm3 | 6.119093e-01 | 0.000000e+00 | -- | False | False |
| jet_leptonic | gamma0_log_parab | turn-over-energy | lorentz-factor* | 3.290961e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | s | LE_spectral_slope | 2.169388e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_leptonic | r | spectral_curvature | 7.818737e-01 | -1.500000e+01 | 1.500000e+01 | False | False |
================================================================================
from jetset.minimizer import fit_SED,ModelMinimizer
from jetset.model_manager import FitModel
from jetset.jet_model import Jet
prefit_jet=Jet.load_model('prefit_jet.pkl')
fit_model=FitModel( jet=prefit_jet, name='SSC-best-fit-lsb',template=None)
fit_model.parameters
| 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.697542e+02 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 1.300733e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_leptonic | N | emitters_density | 1 / cm3 | 6.119093e-01 | 0.000000e+00 | -- | False | False |
| jet_leptonic | gamma0_log_parab | turn-over-energy | lorentz-factor* | 3.290961e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | s | LE_spectral_slope | 2.169388e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_leptonic | r | spectral_curvature | 7.818737e-01 | -1.500000e+01 | 1.500000e+01 | False | False | |
| jet_leptonic | R | region_size | cm | 3.452668e+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 | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_leptonic | beam_obj | beaming | 2.500000e+01 | 1.000000e-04 | -- | False | False | |
| jet_leptonic | z_cosm | redshift | 3.080000e-02 | 0.000000e+00 | -- | False | False |
None
fit_model.jet_leptonic.parameters.beam_obj.fit_range = [5, 50]
fit_model.jet_leptonic.parameters.R_H.val=5E17
fit_model.jet_leptonic.parameters.R_H.frozen=False
fit_model.jet_leptonic.parameters.R_H.fit_range = [1E15, 1E19]
fit_model.jet_leptonic.parameters.R.fit_range = [10 ** 15.5, 10 ** 17.5]
fit_model.jet_leptonic.parameters.gamma0_log_parab.fit_range = [1E3,1E6]
fit_model.jet_leptonic.parameters.gmin.fit_range = [10,1000]
fit_model.jet_leptonic.parameters.gmax.fit_range = [1E5,1E8]
fit_model.jet_leptonic.add_user_par(name='B0',units='G',val=1E3,val_min=0,val_max=None)
fit_model.jet_leptonic.add_user_par(name='R0', units='cm', val=5E13, val_min=0, val_max=None)
fit_model.jet_leptonic.add_user_par(name='m_B', val=1, val_min=1, val_max=2)
fit_model.jet_leptonic.parameters.R0.frozen=True
fit_model.jet_leptonic.parameters.B0.frozen=True
def par_func(R0,B0,R_H,m_B):
return B0*np.power((R0/R_H),m_B)
fit_model.jet_leptonic.make_dependent_par(par='B', depends_on=['B0', 'R0', 'R_H','m_B'], par_expr=par_func)
fit_model.parameters
adding par: B0 to B
adding par: R0 to B
adding par: R_H to B
adding par: m_B to B
==> par B is depending on ['B0', 'R0', 'R_H', 'm_B'] according to expr: B =
def par_func(R0,B0,R_H,m_B):
return B0*np.power((R0/R_H),m_B)
| 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.697542e+02 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 1.300733e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_leptonic | N | emitters_density | 1 / cm3 | 6.119093e-01 | 0.000000e+00 | -- | False | False |
| jet_leptonic | gamma0_log_parab | turn-over-energy | lorentz-factor* | 3.290961e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | s | LE_spectral_slope | 2.169388e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_leptonic | r | spectral_curvature | 7.818737e-01 | -1.500000e+01 | 1.500000e+01 | False | False | |
| jet_leptonic | R | region_size | cm | 3.452668e+16 | 1.000000e+03 | 1.000000e+30 | False | False |
| jet_leptonic | R_H(M) | region_position | cm | 5.000000e+17 | 0.000000e+00 | -- | False | False |
| jet_leptonic | *B(D,m_B) | magnetic_field | gauss | 1.000000e-01 | 0.000000e+00 | -- | False | True |
| jet_leptonic | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_leptonic | beam_obj | beaming | 2.500000e+01 | 1.000000e-04 | -- | False | False | |
| jet_leptonic | z_cosm | redshift | 3.080000e-02 | 0.000000e+00 | -- | False | False | |
| jet_leptonic | B0(M) | user_defined | G | 1.000000e+03 | 0.000000e+00 | -- | False | True |
| jet_leptonic | R0(M) | user_defined | cm | 5.000000e+13 | 0.000000e+00 | -- | False | True |
| jet_leptonic | m_B(M) | user_defined | 1.000000e+00 | 1.000000e+00 | 2.000000e+00 | False | False |
None
%matplotlib inline
plt.figure(dpi=150)
R_H_array=np.logspace(13,18,100)
B_array=np.zeros(R_H_array.shape)
for ID,R_H in enumerate(R_H_array):
fit_model.jet_leptonic.parameters.R_H.val=R_H
B_array[ID]=fit_model.jet_leptonic.parameters.B.val
plt.loglog(R_H_array,B_array)
plt.xlabel('R_H (cm)')
plt.ylabel('B (G)')
Text(0, 0.5, 'B (G)')
fit_model.jet_leptonic.parameters.R_H.val=5E17
fit_model.parameters
| 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.697542e+02 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 1.300733e+06 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_leptonic | N | emitters_density | 1 / cm3 | 6.119093e-01 | 0.000000e+00 | -- | False | False |
| jet_leptonic | gamma0_log_parab | turn-over-energy | lorentz-factor* | 3.290961e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | s | LE_spectral_slope | 2.169388e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_leptonic | r | spectral_curvature | 7.818737e-01 | -1.500000e+01 | 1.500000e+01 | False | False | |
| jet_leptonic | R | region_size | cm | 3.452668e+16 | 1.000000e+03 | 1.000000e+30 | False | False |
| jet_leptonic | R_H(M) | region_position | cm | 5.000000e+17 | 0.000000e+00 | -- | False | False |
| jet_leptonic | *B(D,m_B) | magnetic_field | gauss | 1.000000e-01 | 0.000000e+00 | -- | False | True |
| jet_leptonic | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_leptonic | beam_obj | beaming | 2.500000e+01 | 1.000000e-04 | -- | False | False | |
| jet_leptonic | z_cosm | redshift | 3.080000e-02 | 0.000000e+00 | -- | False | False | |
| jet_leptonic | B0(M) | user_defined | G | 1.000000e+03 | 0.000000e+00 | -- | False | True |
| jet_leptonic | R0(M) | user_defined | cm | 5.000000e+13 | 0.000000e+00 | -- | False | True |
| jet_leptonic | m_B(M) | user_defined | 1.000000e+00 | 1.000000e+00 | 2.000000e+00 | False | False |
None
As a resuslt of the best fit modeling, we are able to determine the
value of R_H. We now perform the fit with minuit to get a better
estimate of the errors
model_minimizer_minuit=ModelMinimizer('minuit')
model_minimizer_minuit.minimizer.add_simplex=False
best_fit_minuit=model_minimizer_minuit.fit(fit_model,
sed_data,
1E11,
1E29,
fitname='SSC-best-fit-minuit',
repeat=3)
filtering data in fit range = [1.000000e+11,1.000000e+29] data length 34 ================================================================================ * start fit process * ----- fit run: 0
0it [00:00, ?it/s]
- best chisq=4.32456e+01
fit run: 1
- old chisq=4.32456e+01
0it [00:00, ?it/s]
- best chisq=1.77754e+01
fit run: 2
- old chisq=1.77754e+01
0it [00:00, ?it/s]
- best chisq=1.76015e+01
-------------------------------------------------------------------------
Fit report
Model: SSC-best-fit-minuit
| model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
|---|---|---|---|---|---|---|---|---|
| jet_leptonic | gmin | low-energy-cut-off | lorentz-factor* | 6.394937e+02 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 6.862548e+05 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_leptonic | N | emitters_density | 1 / cm3 | 7.478459e-01 | 0.000000e+00 | -- | False | False |
| jet_leptonic | gamma0_log_parab | turn-over-energy | lorentz-factor* | 2.800712e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | s | LE_spectral_slope | 2.194676e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_leptonic | r | spectral_curvature | 6.353119e-01 | -1.500000e+01 | 1.500000e+01 | False | False | |
| jet_leptonic | R | region_size | cm | 2.402660e+16 | 1.000000e+03 | 1.000000e+30 | False | False |
| jet_leptonic | R_H(M) | region_position | cm | 5.096309e+17 | 0.000000e+00 | -- | False | False |
| jet_leptonic | *B(D,m_B) | magnetic_field | gauss | 6.210883e-02 | 0.000000e+00 | -- | False | True |
| jet_leptonic | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_leptonic | beam_obj | beaming | 3.117104e+01 | 1.000000e-04 | -- | False | False | |
| jet_leptonic | z_cosm | redshift | 4.355533e-02 | 0.000000e+00 | -- | False | False | |
| jet_leptonic | B0(M) | user_defined | G | 1.000000e+03 | 0.000000e+00 | -- | False | True |
| jet_leptonic | R0(M) | user_defined | cm | 5.000000e+13 | 0.000000e+00 | -- | False | True |
| jet_leptonic | m_B(M) | user_defined | 1.049538e+00 | 1.000000e+00 | 2.000000e+00 | False | False |
converged=True
calls=2953
mesg=
| Migrad | |
|---|---|
| FCN = 17.6 | Nfcn = 2953 |
| EDM = 7.09 (Goal: 0.0002) | time = 11.5 sec |
| INVALID Minimum | ABOVE EDM threshold (goal x 10) |
| No parameters at limit | Below call limit |
| Hesse ok | Covariance FORCED pos. def. |
| Name | Value | Hesse Error | Minos Error- | Minos Error+ | Limit- | Limit+ | Fixed | |
|---|---|---|---|---|---|---|---|---|
| 0 | par_0 | 640 | 60 | 10 | 1E+03 | |||
| 1 | par_1 | 0.69e6 | 0.04e6 | 1E+05 | 1E+08 | |||
| 2 | par_2 | 747.85e-3 | 0.24e-3 | 0 | ||||
| 3 | par_3 | 28.0e3 | 1.1e3 | 1E+03 | 1E+06 | |||
| 4 | par_4 | 2.195 | 0.032 | -10 | 10 | |||
| 5 | par_5 | 635.31e-3 | 0.26e-3 | -15 | 15 | |||
| 6 | par_6 | 24.0e15 | 0.4e15 | 3.16E+15 | 3.16E+17 | |||
| 7 | par_7 | 509.63e15 | 0.26e15 | 1E+15 | 1E+19 | |||
| 8 | par_8 | 31.2 | 0.4 | 5 | 50 | |||
| 9 | par_9 | 0.0436 | 0.0010 | 0 | ||||
| 10 | par_10 | 1.04954 | 0.00005 | 1 | 2 |
dof=23
chisq=17.601512, chisq/red=0.765283 null hypothesis sig=0.778731
best fit pars
| model name | name | val | bestfit val | err + | err - | start val | fit range min | fit range max | frozen |
|---|---|---|---|---|---|---|---|---|---|
| jet_leptonic | gmin | 6.394937e+02 | 6.394937e+02 | 6.397673e+01 | -- | 4.697542e+02 | 1.000000e+01 | 1.000000e+03 | False |
| jet_leptonic | gmax | 6.862548e+05 | 6.862548e+05 | 4.102357e+04 | -- | 1.300733e+06 | 1.000000e+05 | 1.000000e+08 | False |
| jet_leptonic | N | 7.478459e-01 | 7.478459e-01 | 2.416723e-04 | -- | 6.119093e-01 | 0.000000e+00 | -- | False |
| jet_leptonic | gamma0_log_parab | 2.800712e+04 | 2.800712e+04 | 1.096038e+03 | -- | 3.290961e+04 | 1.000000e+03 | 1.000000e+06 | False |
| jet_leptonic | s | 2.194676e+00 | 2.194676e+00 | 3.229843e-02 | -- | 2.169388e+00 | -1.000000e+01 | 1.000000e+01 | False |
| jet_leptonic | r | 6.353119e-01 | 6.353119e-01 | 2.577970e-04 | -- | 7.818737e-01 | -1.500000e+01 | 1.500000e+01 | False |
| jet_leptonic | R | 2.402660e+16 | 2.402660e+16 | 3.922868e+14 | -- | 3.452668e+16 | 3.162278e+15 | 3.162278e+17 | False |
| jet_leptonic | R_H(M) | 5.096309e+17 | 5.096309e+17 | 2.609058e+14 | -- | 5.000000e+17 | 1.000000e+15 | 1.000000e+19 | False |
| jet_leptonic | *B(D,m_B) | 6.210883e-02 | -- | -- | -- | 1.000000e-01 | 0.000000e+00 | -- | True |
| jet_leptonic | NH_cold_to_rel_e | 1.000000e+00 | -- | -- | -- | 1.000000e+00 | 0.000000e+00 | -- | True |
| jet_leptonic | beam_obj | 3.117104e+01 | 3.117104e+01 | 3.601592e-01 | -- | 2.500000e+01 | 5.000000e+00 | 5.000000e+01 | False |
| jet_leptonic | z_cosm | 4.355533e-02 | 4.355533e-02 | 1.047168e-03 | -- | 3.080000e-02 | 0.000000e+00 | -- | False |
| jet_leptonic | B0(M) | 1.000000e+03 | -- | -- | -- | 1.000000e+03 | 0.000000e+00 | -- | True |
| jet_leptonic | R0(M) | 5.000000e+13 | -- | -- | -- | 5.000000e+13 | 0.000000e+00 | -- | True |
| jet_leptonic | m_B(M) | 1.049538e+00 | 1.049538e+00 | 5.070712e-05 | -- | 1.000000e+00 | 1.000000e+00 | 2.000000e+00 | False |
-------------------------------------------------------------------------
================================================================================
fit_model.plot_model(sed_data=sed_data)
<jetset.plot_sedfit.PlotSED at 0x144a7e180>
%matplotlib inline
plt.figure(dpi=150)
R_H_array=np.logspace(13,18,100)
B_array=np.zeros(R_H_array.shape)
for ID,R_H in enumerate(R_H_array):
fit_model.jet_leptonic.parameters.R_H.val=R_H
B_array[ID]=fit_model.jet_leptonic.parameters.B.val
plt.loglog(R_H_array,B_array)
plt.xlabel('R_H (cm)')
plt.ylabel('B (G)')
Text(0, 0.5, 'B (G)')
fit_model.save_model('test.pkl')
from jetset.model_manager import FitModel
new_fit_model=FitModel.load_model('test.pkl')
new_fit_model.parameters
| model name | name | par type | units | val | phys. bound. min | phys. bound. max | log | frozen |
|---|---|---|---|---|---|---|---|---|
| jet_leptonic | gmin | low-energy-cut-off | lorentz-factor* | 6.394937e+02 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | gmax | high-energy-cut-off | lorentz-factor* | 6.862548e+05 | 1.000000e+00 | 1.000000e+15 | False | False |
| jet_leptonic | N | emitters_density | 1 / cm3 | 7.478459e-01 | 0.000000e+00 | -- | False | False |
| jet_leptonic | gamma0_log_parab | turn-over-energy | lorentz-factor* | 2.800712e+04 | 1.000000e+00 | 1.000000e+09 | False | False |
| jet_leptonic | s | LE_spectral_slope | 2.194676e+00 | -1.000000e+01 | 1.000000e+01 | False | False | |
| jet_leptonic | r | spectral_curvature | 6.353119e-01 | -1.500000e+01 | 1.500000e+01 | False | False | |
| jet_leptonic | R | region_size | cm | 2.402660e+16 | 1.000000e+03 | 1.000000e+30 | False | False |
| jet_leptonic | R_H(M) | region_position | cm | 1.000000e+18 | 0.000000e+00 | -- | False | False |
| jet_leptonic | *B(D,m_B) | magnetic_field | gauss | 3.061310e-02 | 0.000000e+00 | -- | False | True |
| jet_leptonic | NH_cold_to_rel_e | cold_p_to_rel_e_ratio | 1.000000e+00 | 0.000000e+00 | -- | False | True | |
| jet_leptonic | beam_obj | beaming | 3.117104e+01 | 1.000000e-04 | -- | False | False | |
| jet_leptonic | z_cosm | redshift | 4.355533e-02 | 0.000000e+00 | -- | False | False | |
| jet_leptonic | B0(M) | user_defined | G | 1.000000e+03 | 0.000000e+00 | -- | False | True |
| jet_leptonic | R0(M) | user_defined | cm | 5.000000e+13 | 0.000000e+00 | -- | False | True |
| jet_leptonic | m_B(M) | user_defined | 1.049538e+00 | 1.000000e+00 | 2.000000e+00 | False | False |
None