.. _dep_pars: 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 --------------------------------- .. code:: ipython3 import warnings warnings.filterwarnings('ignore') .. code:: ipython3 import numpy as np .. code:: ipython3 import jetset print('tested on jetset',jetset.__version__) .. parsed-literal:: tested on jetset 1.3.0rc7 .. code:: ipython3 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) .. parsed-literal:: ===> setting C threads to 12 .. raw:: html Table length=14
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_flaringRregion_sizecm5.000000e+151.000000e+031.000000e+30FalseFalse
jet_flaringR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_flaringBmagnetic_fieldgauss1.000000e-010.000000e+00--FalseFalse
jet_flaringNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_flaringbeam_objbeaming1.000000e+011.000000e-04--FalseFalse
jet_flaringz_cosmredshift1.000000e-020.000000e+00--FalseFalse
jet_flaringgminlow-energy-cut-offlorentz-factor*2.000000e+001.000000e+001.000000e+09FalseFalse
jet_flaringgmaxhigh-energy-cut-offlorentz-factor*1.000000e+061.000000e+001.000000e+15FalseFalse
jet_flaringNemitters_density1 / cm31.000000e+020.000000e+00--FalseFalse
jet_flaringgamma0_log_parabturn-over-energylorentz-factor*1.000000e+041.000000e+001.000000e+09FalseFalse
jet_flaringsLE_spectral_slope2.000000e+00-1.000000e+011.000000e+01FalseFalse
jet_flaringrspectral_curvature4.000000e-01-1.500000e+011.500000e+01FalseFalse
Franceschini_2008scale_factorscale_factor1.000000e+000.000000e+00--FalseTrue
Franceschini_2008z_cosmredshift1.000000e+000.000000e+00--FalseTrue
.. parsed-literal:: adding par: z_cosm to z_cosm ===> setting C threads to 12 adding par: z_cosm to z_cosm .. raw:: html Table length=14
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_flaringgminlow-energy-cut-offlorentz-factor*2.000000e+001.000000e+001.000000e+09FalseFalse
jet_flaringgmaxhigh-energy-cut-offlorentz-factor*1.000000e+061.000000e+001.000000e+15FalseFalse
jet_flaringNemitters_density1 / cm31.000000e+020.000000e+00--FalseFalse
jet_flaringgamma0_log_parabturn-over-energylorentz-factor*1.000000e+041.000000e+001.000000e+09FalseFalse
jet_flaringsLE_spectral_slope2.000000e+00-1.000000e+011.000000e+01FalseFalse
jet_flaringrspectral_curvature4.000000e-01-1.500000e+011.500000e+01FalseFalse
jet_flaringRregion_sizecm5.000000e+151.000000e+031.000000e+30FalseFalse
jet_flaringR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_flaringBmagnetic_fieldgauss1.000000e-010.000000e+00--FalseFalse
jet_flaringNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_flaringbeam_objbeaming1.000000e+011.000000e-04--FalseFalse
jet_flaringz_cosm(M)redshift3.001000e-020.000000e+00--FalseFalse
Franceschini_2008scale_factorscale_factor1.000000e+000.000000e+00--FalseTrue
Franceschini_2008z_cosm(L,jet_flaring)redshift------FalseTrue
.. parsed-literal:: 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. .. code:: ipython3 from jetset.jet_emitters import EmittersDistribution import numpy as np from jetset.jet_model import Jet j = Jet(emitters_distribution='bkn') j.parameters .. parsed-literal:: ===> setting C threads to 12 .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm5.000000e+151.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss1.000000e-010.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming1.000000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift1.000000e-010.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*2.000000e+001.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.000000e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm31.000000e+020.000000e+00--FalseFalse
jet_leptonicgamma_breakturn-over-energylorentz-factor*1.000000e+041.000000e+001.000000e+09FalseFalse
jet_leptonicpLE_spectral_slope2.500000e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicp_1HE_spectral_slope3.500000e+00-1.000000e+011.000000e+01FalseFalse
.. parsed-literal:: 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: .. code:: ipython3 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 .. parsed-literal:: 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 .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm5.000000e+151.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss1.000000e-010.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming1.000000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift1.000000e-010.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*2.000000e+001.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.000000e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm31.000000e+020.000000e+00--FalseFalse
jet_leptonicgamma_breakturn-over-energylorentz-factor*1.000000e+041.000000e+001.000000e+09FalseFalse
jet_leptonicp(M)LE_spectral_slope2.000000e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonic*p_1(D,p)HE_spectral_slope3.000000e+00-1.000000e+011.000000e+01FalseTrue
.. parsed-literal:: 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 .. code:: ipython3 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 .. parsed-literal:: ==> par p_1 is depending on ['p'] according to expr: p_1 = p+1 .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm5.000000e+151.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss1.000000e-010.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming1.000000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift1.000000e-010.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*2.000000e+001.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.000000e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm31.000000e+020.000000e+00--FalseFalse
jet_leptonicgamma_breakturn-over-energylorentz-factor*1.000000e+041.000000e+001.000000e+09FalseFalse
jet_leptonicp(M)LE_spectral_slope2.000000e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonic*p_1(D,p)HE_spectral_slope3.000000e+00-1.000000e+011.000000e+01FalseTrue
.. parsed-literal:: 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: .. code:: ipython3 #j.parameters.p_1.par_expression_source_code .. code:: ipython3 j.save_model('jet.pkl') new_jet=Jet.load_model('jet.pkl') .. parsed-literal:: ===> setting C threads to 12 adding par: p to p_1 ==> par p_1 is depending on ['p'] according to expr: p_1 = p+1 .. code:: ipython3 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 .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*2.000000e+001.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.000000e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm31.000000e+020.000000e+00--FalseFalse
jet_leptonicgamma_breakturn-over-energylorentz-factor*1.000000e+041.000000e+001.000000e+09FalseFalse
jet_leptonicp(M)LE_spectral_slope2.500000e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonic*p_1(D,p)HE_spectral_slope3.500000e+00-1.000000e+011.000000e+01FalseTrue
jet_leptonicRregion_sizecm5.000000e+151.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss1.000000e-010.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming1.000000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift1.000000e-010.000000e+00--FalseFalse
.. parsed-literal:: 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: :math:`B=B_0(R_0/R_H)^{1.1}` .. code:: ipython3 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)) .. parsed-literal:: ===> setting C threads to 12 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) .. code:: ipython3 fit_model_lsb.jet_leptonic.parameters .. raw:: html Table length=14
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm5.000000e+151.000000e+031.000000e+30FalseFalse
jet_leptonicR_H(M)region_positioncm5.000000e+170.000000e+00--FalseFalse
jet_leptonic*B(D,m_B)magnetic_fieldgauss1.000000e-010.000000e+00--FalseTrue
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming1.000000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift1.000000e-010.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*2.000000e+001.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.000000e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm31.000000e+020.000000e+00--FalseFalse
jet_leptonicgamma_cutturn-over-energylorentz-factor*1.000000e+041.000000e+001.000000e+09FalseFalse
jet_leptonicpLE_spectral_slope2.000000e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicB0(M)user_definedG1.000000e+030.000000e+00--FalseTrue
jet_leptonicR0(M)user_definedcm5.000000e+130.000000e+00--FalseTrue
jet_leptonicm_B(M)user_defined1.000000e+001.000000e+002.000000e+00FalseFalse
.. parsed-literal:: None .. code:: ipython3 fit_model_lsb.jet_leptonic.parameters .. raw:: html Table length=14
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm5.000000e+151.000000e+031.000000e+30FalseFalse
jet_leptonicR_H(M)region_positioncm5.000000e+170.000000e+00--FalseFalse
jet_leptonic*B(D,m_B)magnetic_fieldgauss1.000000e-010.000000e+00--FalseTrue
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming1.000000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift1.000000e-010.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*2.000000e+001.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.000000e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm31.000000e+020.000000e+00--FalseFalse
jet_leptonicgamma_cutturn-over-energylorentz-factor*1.000000e+041.000000e+001.000000e+09FalseFalse
jet_leptonicpLE_spectral_slope2.000000e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicB0(M)user_definedG1.000000e+030.000000e+00--FalseTrue
jet_leptonicR0(M)user_definedcm5.000000e+130.000000e+00--FalseTrue
jet_leptonicm_B(M)user_defined1.000000e+001.000000e+002.000000e+00FalseFalse
.. parsed-literal:: None .. code:: ipython3 fit_model_lsb.save_model('test.pkl') .. code:: ipython3 fit_model_lsb=FitModel.load_model('test.pkl') .. parsed-literal:: ===> setting C threads to 12 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) .. code:: ipython3 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)) 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 .. code:: ipython3 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 .. code:: ipython3 data=Data.from_file(test_SEDs[1]) .. code:: ipython3 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() .. parsed-literal:: ================================================================================ *** binning data *** ---> N bins= 89 ---> bin_widht= 0.2 ================================================================================ .. image:: depending_pars_files/depending_pars_33_1.png .. code:: ipython3 from jetset.sed_shaper import SEDShape my_shape=SEDShape(sed_data) my_shape.eval_indices(minimizer='lsb',silent=True) p=my_shape.plot_indices() .. parsed-literal:: ================================================================================ *** evaluating spectral indices for data *** ================================================================================ .. image:: depending_pars_files/depending_pars_34_1.png .. code:: ipython3 mm,best_fit=my_shape.sync_fit(check_host_gal_template=False, Ep_start=None, minimizer='lsb', silent=True, fit_range=[10.,21.]) .. parsed-literal:: ================================================================================ *** Log-Polynomial fitting of the synchrotron component *** ---> first blind fit run, fit range: [10.0, 21.0] ---> class: HSP .. raw:: html Table length=4
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
LogCubicb-1.585748e-01-1.585748e-016.470535e-03---1.000000e+00-1.000000e+010.000000e+00False
LogCubicc-1.089513e-02-1.089513e-029.764985e-04---1.000000e+00-1.000000e+011.000000e+01False
LogCubicEp1.673177e+011.673177e+012.478677e-02--1.667298e+010.000000e+003.000000e+01False
LogCubicSp-9.489417e+00-9.489417e+001.853260e-02---1.000000e+01-3.000000e+010.000000e+00False
.. parsed-literal:: ---> sync nu_p=+1.673177e+01 (err=+2.478677e-02) nuFnu_p=-9.489417e+00 (err=+1.853260e-02) curv.=-1.585748e-01 (err=+6.470535e-03) ================================================================================ .. code:: ipython3 my_shape.IC_fit(fit_range=[23.,29.],minimizer='minuit',silent=True) p=my_shape.plot_shape_fit() p.setlim(y_min=1E-15) .. parsed-literal:: ================================================================================ *** Log-Polynomial fitting of the IC component *** ---> fit range: [23.0, 29.0] ---> LogCubic fit ====> simplex ====> migrad ====> simplex ====> migrad ====> simplex ====> migrad .. raw:: html Table length=4
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
LogCubicb-1.971111e-01-1.971111e-012.679732e-02---1.000000e+00-1.000000e+010.000000e+00False
LogCubicc-4.037544e-02-4.037544e-022.119803e-02---1.000000e+00-1.000000e+011.000000e+01False
LogCubicEp2.521789e+012.521789e+011.198160e-01--2.529262e+010.000000e+003.000000e+01False
LogCubicSp-1.012535e+01-1.012535e+012.996508e-02---1.000000e+01-3.000000e+010.000000e+00False
.. parsed-literal:: ---> IC nu_p=+2.521789e+01 (err=+1.198160e-01) nuFnu_p=-1.012535e+01 (err=+2.996508e-02) curv.=-1.971111e-01 (err=+2.679732e-02) ================================================================================ .. image:: depending_pars_files/depending_pars_36_3.png .. code:: ipython3 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') .. parsed-literal:: ================================================================================ *** constrains parameters from observable *** ===> setting C threads to 12 .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicRregion_sizecm3.460321e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss5.050000e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming2.500000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.080000e-020.000000e+00--FalseFalse
jet_leptonicgminlow-energy-cut-offlorentz-factor*4.697542e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.373160e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm36.545152e-010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*3.333017e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.183468e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature7.928739e-01-1.500000e+011.500000e+01FalseFalse
.. parsed-literal:: ================================================================================ .. code:: ipython3 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') .. parsed-literal:: ===> setting C threads to 12 .. code:: ipython3 fit_model=FitModel( jet=prefit_jet, name='SSC-best-fit-lsb',template=None) fit_model.parameters .. raw:: html Table length=12
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*4.697542e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.373160e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm36.545152e-010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*3.333017e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.183468e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature7.928739e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm3.460321e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm1.000000e+170.000000e+00--FalseTrue
jet_leptonicBmagnetic_fieldgauss5.050000e-020.000000e+00--FalseFalse
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming2.500000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.080000e-020.000000e+00--FalseFalse
.. parsed-literal:: None .. code:: ipython3 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 .. parsed-literal:: 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) .. raw:: html Table length=15
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*4.697542e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.373160e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm36.545152e-010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*3.333017e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.183468e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature7.928739e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm3.460321e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_H(M)region_positioncm5.000000e+170.000000e+00--FalseFalse
jet_leptonic*B(D,m_B)magnetic_fieldgauss1.000000e-010.000000e+00--FalseTrue
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming2.500000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.080000e-020.000000e+00--FalseFalse
jet_leptonicB0(M)user_definedG1.000000e+030.000000e+00--FalseTrue
jet_leptonicR0(M)user_definedcm5.000000e+130.000000e+00--FalseTrue
jet_leptonicm_B(M)user_defined1.000000e+001.000000e+002.000000e+00FalseFalse
.. parsed-literal:: None .. code:: ipython3 %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)') .. parsed-literal:: Text(0, 0.5, 'B (G)') .. image:: depending_pars_files/depending_pars_41_1.png .. code:: ipython3 fit_model.jet_leptonic.parameters.R_H.val=5E17 .. code:: ipython3 fit_model.parameters .. raw:: html Table length=15
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*4.697542e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*1.373160e+061.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm36.545152e-010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*3.333017e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.183468e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature7.928739e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm3.460321e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_H(M)region_positioncm5.000000e+170.000000e+00--FalseFalse
jet_leptonic*B(D,m_B)magnetic_fieldgauss1.000000e-010.000000e+00--FalseTrue
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming2.500000e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift3.080000e-020.000000e+00--FalseFalse
jet_leptonicB0(M)user_definedG1.000000e+030.000000e+00--FalseTrue
jet_leptonicR0(M)user_definedcm5.000000e+130.000000e+00--FalseTrue
jet_leptonicm_B(M)user_defined1.000000e+001.000000e+002.000000e+00FalseFalse
.. parsed-literal:: 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 .. code:: ipython3 model_minimizer_minuit=ModelMinimizer('minuit') .. code:: ipython3 model_minimizer_minuit.minimizer.add_simplex=False .. code:: ipython3 best_fit_minuit=model_minimizer_minuit.fit(fit_model, sed_data, 1E11, 1E29, fitname='SSC-best-fit-minuit', repeat=3) .. parsed-literal:: filtering data in fit range = [1.000000e+11,1.000000e+29] data length 35 ================================================================================ *** start fit process *** ----- fit run: 0 .. parsed-literal:: 0it [00:00, ?it/s] .. parsed-literal:: ====> migrad - best chisq=3.30494e+01 fit run: 1 - old chisq=3.30494e+01 .. parsed-literal:: 0it [00:00, ?it/s] .. parsed-literal:: ====> migrad - best chisq=2.21161e+01 fit run: 2 - old chisq=2.21161e+01 .. parsed-literal:: 0it [00:00, ?it/s] .. parsed-literal:: ====> migrad - best chisq=2.01970e+01 ------------------------------------------------------------------------- Fit report Model: SSC-best-fit-minuit .. raw:: html Table length=15
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*5.955908e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*7.336491e+051.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm37.078650e-010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*3.225336e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.232280e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature6.488585e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm2.940961e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_H(M)region_positioncm4.786191e+170.000000e+00--FalseFalse
jet_leptonic*B(D,m_B)magnetic_fieldgauss4.766679e-020.000000e+00--FalseTrue
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming3.726640e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift5.799018e-020.000000e+00--FalseFalse
jet_leptonicB0(M)user_definedG1.000000e+030.000000e+00--FalseTrue
jet_leptonicR0(M)user_definedcm5.000000e+130.000000e+00--FalseTrue
jet_leptonicm_B(M)user_defined1.085597e+001.000000e+002.000000e+00FalseFalse
.. parsed-literal:: converged=True calls=4803 mesg= .. raw:: html
Migrad
FCN = 20.2 Nfcn = 4803
EDM = 3.32e-07 (Goal: 0.0002) time = 109.8 sec
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 par_0 595.6 1.1 10 1E+03
1 par_1 733.6e3 0.5e3 1E+05 1E+08
2 par_2 707.865e-3 0.010e-3 0
3 par_3 32.2534e3 0.0022e3 1E+03 1E+06
4 par_4 2.232280 0.000023 -10 10
5 par_5 648.86e-3 0.26e-3 -15 15
6 par_6 29.4096e15 0.0012e15 3.16E+15 3.16E+17
7 par_7 478.619e15 0.027e15 1E+15 1E+19
8 par_8 37.26640 0.00024 5 50
9 par_9 57.9902e-3 0.0034e-3 0
10 par_10 1.0855972 0.0000033 1 2
.. parsed-literal:: dof=24 chisq=20.196967, chisq/red=0.841540 null hypothesis sig=0.685522 best fit pars .. raw:: html Table length=15
model namenamevalbestfit valerr +err -start valfit range minfit range maxfrozen
jet_leptonicgmin5.955908e+025.955908e+021.099191e+00--4.697542e+021.000000e+011.000000e+03False
jet_leptonicgmax7.336491e+057.336491e+055.239237e+02--1.373160e+061.000000e+051.000000e+08False
jet_leptonicN7.078650e-017.078650e-011.030217e-05--6.545152e-010.000000e+00--False
jet_leptonicgamma0_log_parab3.225336e+043.225336e+042.156816e+00--3.333017e+041.000000e+031.000000e+06False
jet_leptonics2.232280e+002.232280e+002.282556e-05--2.183468e+00-1.000000e+011.000000e+01False
jet_leptonicr6.488585e-016.488585e-012.582949e-04--7.928739e-01-1.500000e+011.500000e+01False
jet_leptonicR2.940961e+162.940961e+161.250211e+12--3.460321e+163.162278e+153.162278e+17False
jet_leptonicR_H(M)4.786191e+174.786191e+172.668558e+13--5.000000e+171.000000e+151.000000e+19False
jet_leptonic*B(D,m_B)4.766679e-02------1.000000e-010.000000e+00--True
jet_leptonicNH_cold_to_rel_e1.000000e+00------1.000000e+000.000000e+00--True
jet_leptonicbeam_obj3.726640e+013.726640e+012.449689e-04--2.500000e+015.000000e+005.000000e+01False
jet_leptonicz_cosm5.799018e-025.799018e-023.432852e-06--3.080000e-020.000000e+00--False
jet_leptonicB0(M)1.000000e+03------1.000000e+030.000000e+00--True
jet_leptonicR0(M)5.000000e+13------5.000000e+130.000000e+00--True
jet_leptonicm_B(M)1.085597e+001.085597e+003.331357e-06--1.000000e+001.000000e+002.000000e+00False
.. parsed-literal:: ------------------------------------------------------------------------- ================================================================================ .. code:: ipython3 fit_model.plot_model(sed_data=sed_data) .. parsed-literal:: .. image:: depending_pars_files/depending_pars_48_1.png .. code:: ipython3 %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)') .. parsed-literal:: Text(0, 0.5, 'B (G)') .. image:: depending_pars_files/depending_pars_49_1.png .. code:: ipython3 fit_model.save_model('test.pkl') .. code:: ipython3 from jetset.model_manager import FitModel new_fit_model=FitModel.load_model('test.pkl') .. parsed-literal:: ===> setting C threads to 12 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) .. code:: ipython3 new_fit_model.parameters .. raw:: html Table length=15
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*5.955908e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*7.336491e+051.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm37.078650e-010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*3.225336e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.232280e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature6.488585e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm2.940961e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_H(M)region_positioncm4.786191e+170.000000e+00--FalseFalse
jet_leptonic*B(D,m_B)magnetic_fieldgauss4.766679e-020.000000e+00--FalseTrue
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming3.726640e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift5.799018e-020.000000e+00--FalseFalse
jet_leptonicB0(M)user_definedG1.000000e+030.000000e+00--FalseTrue
jet_leptonicR0(M)user_definedcm5.000000e+130.000000e+00--FalseTrue
jet_leptonicm_B(M)user_defined1.085597e+001.000000e+002.000000e+00FalseFalse
.. parsed-literal:: None .. code:: ipython3 new_fit_model.jet_leptonic.parameters.reset_dependencies() .. code:: ipython3 new_fit_model.parameters .. raw:: html Table length=15
model namenamepar typeunitsvalphys. bound. minphys. bound. maxlogfrozen
jet_leptonicgminlow-energy-cut-offlorentz-factor*5.955908e+021.000000e+001.000000e+09FalseFalse
jet_leptonicgmaxhigh-energy-cut-offlorentz-factor*7.336491e+051.000000e+001.000000e+15FalseFalse
jet_leptonicNemitters_density1 / cm37.078650e-010.000000e+00--FalseFalse
jet_leptonicgamma0_log_parabturn-over-energylorentz-factor*3.225336e+041.000000e+001.000000e+09FalseFalse
jet_leptonicsLE_spectral_slope2.232280e+00-1.000000e+011.000000e+01FalseFalse
jet_leptonicrspectral_curvature6.488585e-01-1.500000e+011.500000e+01FalseFalse
jet_leptonicRregion_sizecm2.940961e+161.000000e+031.000000e+30FalseFalse
jet_leptonicR_Hregion_positioncm4.786191e+170.000000e+00--FalseFalse
jet_leptonicBmagnetic_fieldgauss4.766679e-020.000000e+00--FalseTrue
jet_leptonicNH_cold_to_rel_ecold_p_to_rel_e_ratio1.000000e+000.000000e+00--FalseTrue
jet_leptonicbeam_objbeaming3.726640e+011.000000e-04--FalseFalse
jet_leptonicz_cosmredshift5.799018e-020.000000e+00--FalseFalse
jet_leptonicB0user_definedG1.000000e+030.000000e+00--FalseTrue
jet_leptonicR0user_definedcm5.000000e+130.000000e+00--FalseTrue
jet_leptonicm_Buser_defined1.085597e+001.000000e+002.000000e+00FalseFalse
.. parsed-literal:: None