diff --git a/.github/workflows/adapt_reqirements_to_git_action.py b/.github/workflows/adapt_reqirements_to_git_action.py index b7d3999d..8f615324 100644 --- a/.github/workflows/adapt_reqirements_to_git_action.py +++ b/.github/workflows/adapt_reqirements_to_git_action.py @@ -1,5 +1,6 @@ _skip_list=['pyqt','swig'] #_skip_list=['pyqt'] +_skip_list=[] f = open("./requirements.txt",'r') req=f.readlines() f.close() diff --git a/.github/workflows/test-conda-workflow.yml b/.github/workflows/test-conda-workflow.yml index 2820d4d5..c9297b6f 100644 --- a/.github/workflows/test-conda-workflow.yml +++ b/.github/workflows/test-conda-workflow.yml @@ -46,58 +46,54 @@ jobs: fail-fast: false matrix: include: - - os: macOS-12 - name: macOS-x86_64 - cpython-version: 39 - python-version: "3.9" - platform_id: macosx_x86_64 - cibw_manylinux: manylinux2014 - cibw_arch: x86_64 - - os: ubuntu-latest - name: x86_64 - cpython-version: 39 - python-version: "3.9" - platform_id: manylinux_x86_64 - cibw_manylinux: manylinux2014 - cibw_arch: x86_64 - - - os: macOS-12 - name: macOS-x86_64 - cpython-version: 310 - python-version: "3.10" - platform_id: macosx_x86_64 - cibw_manylinux: manylinux2014 - cibw_arch: x86_64 - - os: ubuntu-latest - name: x86_64 - cpython-version: 310 - python-version: "3.10" - platform_id: manylinux_x86_64 - cibw_manylinux: manylinux2014 - cibw_arch: x86_64 - - - - os: macOS-12 - name: macOS-x86_64 - cpython-version: 311 - python-version: "3.11" - platform_id: macosx_x86_64 - cibw_manylinux: manylinux2014 - cibw_arch: x86_64 - - os: ubuntu-latest - name: x86_64 - cpython-version: 311 - python-version: "3.11" - platform_id: manylinux_x86_64 - cibw_manylinux: manylinux2014 - cibw_arch: x86_64 - - os: macOS-14 - name: MAC_ARM - cpython-version: 311 - python-version: "3.11" - platform_id: macosx_arm64 - cibw_manylinux: macosx_arm64 - + - os: macOS-13 + name: macOS-x86_64 + cpython-version: 39 + python-version: "3.9" + platform_id: macosx_x86_64 + cibw_arch: x86_64 + + - os: ubuntu-latest + name: x86_64 + cpython-version: 39 + python-version: "3.9" + platform_id: manylinux_x86_64 + cibw_manylinux: manylinux2014 + cibw_arch: x86_64 + + - os: macOS-13 + name: macOS-x86_64 + cpython-version: 310 + python-version: "3.10" + platform_id: macosx_x86_64 + cibw_arch: x86_64 + - os: ubuntu-latest + name: x86_64 + cpython-version: 310 + python-version: "3.10" + platform_id: manylinux_x86_64 + cibw_manylinux: manylinux2014 + cibw_arch: x86_64 + + - os: macOS-13 + name: macOS-x86_64 + cpython-version: 311 + python-version: "3.11" + platform_id: macosx_x86_64 + cibw_arch: x86_64 + - os: ubuntu-latest + name: x86_64 + cpython-version: 311 + python-version: "3.11" + platform_id: manylinux_x86_64 + cibw_manylinux: manylinux2014 + cibw_arch: x86_64 + - os: macOS-14 + name: MAC_ARM + cpython-version: 311 + python-version: "3.11" + platform_id: macosx_arm64 + cibw_manylinux: macosx_arm64 steps: - name: 'Clone Repository (Latest)' uses: actions/checkout@v4 diff --git a/.github/workflows/test-pip-workflow.yml b/.github/workflows/test-pip-workflow.yml index 666abb9f..ee1f9cbc 100644 --- a/.github/workflows/test-pip-workflow.yml +++ b/.github/workflows/test-pip-workflow.yml @@ -45,13 +45,13 @@ jobs: fail-fast: false matrix: include: - - os: macOS-12 + - os: macOS-13 name: macOS-x86_64 cpython-version: 39 python-version: "3.9" platform_id: macosx_x86_64 - cibw_manylinux: manylinux2014 cibw_arch: x86_64 + - os: ubuntu-latest name: x86_64 cpython-version: 39 @@ -60,12 +60,11 @@ jobs: cibw_manylinux: manylinux2014 cibw_arch: x86_64 - - os: macOS-12 + - os: macOS-13 name: macOS-x86_64 cpython-version: 310 python-version: "3.10" platform_id: macosx_x86_64 - cibw_manylinux: manylinux2014 cibw_arch: x86_64 - os: ubuntu-latest name: x86_64 @@ -75,13 +74,11 @@ jobs: cibw_manylinux: manylinux2014 cibw_arch: x86_64 - - - os: macOS-12 + - os: macOS-13 name: macOS-x86_64 cpython-version: 311 python-version: "3.11" platform_id: macosx_x86_64 - cibw_manylinux: manylinux2014 cibw_arch: x86_64 - os: ubuntu-latest name: x86_64 @@ -98,7 +95,6 @@ jobs: cibw_manylinux: macosx_arm64 - steps: - name: 'Clone Repository (Latest)' uses: actions/checkout@v4 diff --git a/jetkernel_src/src/Energetic.c b/jetkernel_src/src/Energetic.c index 4e424c4d..bef3b6a6 100644 --- a/jetkernel_src/src/Energetic.c +++ b/jetkernel_src/src/Energetic.c @@ -391,6 +391,7 @@ struct jet_energetic EnergeticOutput(struct blob * pt) { //char f_Energetic[static_file_name_max_legth]; //FILE *fp_Energetic; + //Lum factor consistent with Eq. 3 and 4, Ghisellini 2010, doi:10.1111/j.1365-2966.2009.15898.x lum_factor_rad =0.25 *eval_beta_gamma(pt->BulkFactor) * pt->BulkFactor * pt->BulkFactor ; lum_factor = pi * pt->R * pt->R * vluce_cm * eval_beta_gamma(pt->BulkFactor) * pt->BulkFactor * pt->BulkFactor ; energetic.U_B= pt->UB; diff --git a/jetkernel_src/src/PyInterface.c b/jetkernel_src/src/PyInterface.c index 2fae0397..bfdcfea9 100644 --- a/jetkernel_src/src/PyInterface.c +++ b/jetkernel_src/src/PyInterface.c @@ -357,24 +357,30 @@ struct blob MakeBlob() { void set_seed_freq_start(struct blob *pt_base){ //pt_base->nu_start_Sync = min(1e6, pt_base->nu_start_grid); pt_base->nu_start_Sync = 1E6; + pt_base->NU_INT_STOP_Sync_SSC=0; pt_base->nu_stop_Sync = 1E20; pt_base->nu_start_SSC = 1E14; - + pt_base->NU_INT_STOP_COMPTON_SSC=0; + //pt_base->nu_stop_SSC = max(1e30, pt_base->nu_stop_grid); pt_base->nu_stop_SSC=1E30; pt_base->nu_start_EC_Disk = 1E13; //pt_base->nu_stop_EC_Disk = max(1e30, pt_base->nu_stop_grid); pt_base->nu_stop_EC_Disk =1E30; + pt_base->NU_INT_MAX_Disk=0; pt_base->nu_start_EC_BLR = 1E13; //pt_base->nu_stop_EC_BLR = max(1e30, pt_base->nu_stop_grid); pt_base->nu_stop_EC_BLR = 1E30; + pt_base->NU_INT_MAX_BLR=0; pt_base->nu_start_EC_DT = 1E13; + pt_base->NU_INT_MAX_DT=0; pt_base->nu_start_EC_CMB = 1E13; - + pt_base->NU_INT_MAX_CMB=0; + pt_base->NU_INT_MAX_Star=0; //pt_base->nu_stop_EC_CMB = max(1e30, pt_base->nu_stop_grid); pt_base->nu_stop_EC_CMB = 1E30; } @@ -410,8 +416,8 @@ void InitRadiative(struct blob *pt_base,unsigned int update_EC){ //COSTANTI PER ALFA=FIXED E KERNEL DELTA O KERNEL 2 pt_base->C1_Sync_K53 = pow(3, 0.5) * pow(q_esu, 3.0) * pt_base->sin_psi; pt_base->C1_Sync_K53 *= pt_base->B / (MEC2) * one_by_four_pi; - pt_base->C2_Sync_K53 = 2.0/(3*pt_base->nu_B); - + pt_base->C2_Sync_K53 = 2.0/(3*pt_base->nu_B* pt_base->sin_psi); + pt_base->C1_Sync_K_AVE= 4*pi*pow(3, 0.5) * pow(q_esu, 2.0)*pt_base->nu_B/(vluce_cm) * one_by_four_pi; pt_base->C2_Sync_K_AVE=1.0/(3*pt_base->nu_B); diff --git a/jetset/base_model.py b/jetset/base_model.py index 21a90eed..f94f14f9 100644 --- a/jetset/base_model.py +++ b/jetset/base_model.py @@ -185,7 +185,7 @@ def plot_model(self,plot_obj=None,clean=False,sed_data=None,frame='obs',skip_com label = self.name if hasattr(self,'SED'): - plot_obj.add_model_plot(self.SED, line_style=line_style,label =label,flim=self.flux_plot_lim,density=density, frame=frame) + plot_obj.add_model_plot(self.SED, line_style=line_style,label =label,flim=self.flux_plot_lim, frame=frame) @@ -197,13 +197,11 @@ def plot_model(self,plot_obj=None,clean=False,sed_data=None,frame='obs',skip_com line_style = '--' if comp_label!='Sum': if hasattr(c, 'SED'): - plot_obj.add_model_plot(c.SED, line_style=line_style, label=' -%s'%comp_label, flim=self.flux_plot_lim, density=density, frame=frame) + plot_obj.add_model_plot(c.SED, line_style=line_style, label=' -%s'%comp_label, flim=self.flux_plot_lim, frame=frame) line_style = '-' - #plot_obj.add_model_residual_plot(data=sed_data, model=self,fit_range=np.log10([self.nu_min_fit,self.nu_max_fit]) ) - return plot_obj def set_nu_grid(self,nu_min=None,nu_max=None,nu_size=None): @@ -247,14 +245,16 @@ def get_residuals(self, data, log_log=False,filter_UL=True): return np.log10(nu_residuals[msk]), residuals[msk] def save_model(self, file_name): + pickle.dump(self, open(file_name, 'wb'), protocol=pickle.HIGHEST_PROTOCOL) + def save_model(self, file_name): pickle.dump(self, open(file_name, 'wb'), protocol=pickle.HIGHEST_PROTOCOL) @classmethod - def load_model(cls, file_name): + def load_model(cls, file_name_or_obj,from_string=False): try: - c = pickle.load(open(file_name, "rb")) + c=cls._load_pickle(file_name_or_obj,from_string=from_string) c._fix_par_dep_on_load() if isinstance(c, Model): c.eval() @@ -276,18 +276,16 @@ def _fix_par_dep_on_load(self,): p.reset_dependencies() self.make_dependent_par(p.name, _master_par_list, _depending_par_expr,set_par_expr_source_code=True) - - #def _set_pars_dep(self): - # for p in self.parameters.par_array: - # if - - - #def _set_pars_dep(self): - # for p in self.parameters.par_array: - # if - + @staticmethod + def _load_pickle(file_name_or_obj,from_string=False): + if from_string: + c = pickle.loads(file_name_or_obj) + else: + c = pickle.load(open(file_name_or_obj, "rb")) + return c + def clone(self): - return pickle.loads(pickle.dumps(self)) + return self.load_model(pickle.dumps(self, protocol=pickle.HIGHEST_PROTOCOL),from_string=True) def show_model(self): print("") @@ -450,9 +448,11 @@ def build_table(self, restframe='obs'): _names.append('nuLnu_src') _cols.append(self.SED.nuLnu_src) - _meta=dict(model_name=self.name) - _meta['restframe']= restframe - self._SED_table = Table(_cols, names=_names,meta=_meta) + _meta=dict(model_name=self.name) + _meta['restframe']= restframe + self._SED_table = Table(_cols, names=_names,meta=_meta) + else: + self._SED_table = None def sed_table(self, restframe='obs'): diff --git a/jetset/data_loader.py b/jetset/data_loader.py index 41622f88..9cc66041 100644 --- a/jetset/data_loader.py +++ b/jetset/data_loader.py @@ -6,6 +6,7 @@ from astropy.table import Table,Column +from astropy.table import vstack from astropy import units as u from astropy.units import cds import pickle @@ -534,8 +535,8 @@ def _build_data(self,dupl_filter=False): self.data['dnuFnu_data_log']= np.ones(self.data['nu_data_log'].size) * self.fake_error - print ("Warning: error were not provided ") - print (" assigning %f ") + print ("Warning: errors were not provided ") + print (" assigning ") print (" set error with .set_error"%self.fake_error) @@ -881,7 +882,7 @@ def remove_dupl_entries(self,data): - def group_data(self,N_bin=None,bin_width=None,correct_dispersions=True): + def group_data(self,N_bin=None,bin_width=None,correct_dispersions=True,nu_min=None,nu_max=None): """ function to perform a spectral group of the data @@ -904,52 +905,48 @@ def group_data(self,N_bin=None,bin_width=None,correct_dispersions=True): print ("you must provide either N_bin or bin_width") raise ValueError - #if nu_min is None: - xmin= self.data['nu_data_log'].min() * 0.99 - #else: - # xmin= np.log10(nu_min) + if nu_min is None: + xmin= self.data['nu_data'].min() * 0.99 + else: + xmin= nu_min - #if nu_max is None: - xmax= self.data['nu_data_log'].max() * 1.01 - #else: - #xmin = np.log10(nu_max) + if nu_max is None: + xmax= self.data['nu_data'].max() * 1.01 + else: + xmax= nu_max if N_bin is None: - N_bin=int((xmax-xmin)/bin_width) + N_bin=int((np.log10(xmax)-np.log10(xmin))/bin_width) if bin_width is None: - bin_width=(xmax-xmin)/N_bin - - - bin_grid=np.linspace(xmin,xmax,N_bin) - - + bin_width=(np.log10(xmax)-np.log10(xmin))/N_bin + + + bin_grid=np.logspace(np.log10(xmin),np.log10(xmax),N_bin+1) + dx_bin=(bin_grid[1:]-bin_grid[:-1])*0.5 + x_bin=(bin_grid[1:]+bin_grid[:-1])*0.5 + y_bin=np.zeros(x_bin.size) + dy_bin=np.zeros(x_bin.size) + x_UL=np.zeros(0) + dx_UL=np.zeros(0) + dy_UL=np.zeros(0) + y_UL=np.zeros(0) print ("*** binning data ***") print ("---> N bins=",N_bin) - print ("---> bin_widht=",bin_width) - - self.data_reb=self._build_empty_table(n_rows=bin_grid.size) - - - x_bin=np.zeros(bin_grid.size) - y_bin=np.zeros(bin_grid.size) - dx_bin=np.zeros(bin_grid.size) - dy_bin=np.zeros(bin_grid.size) - - #gives the id of element falling in each bin - bin_elements_id=np.digitize(self.data['nu_data_log'], bin_grid) - - - for id in range(bin_grid.size): - msk1=[bin_elements_id==id][0] - #Remove UL - msk2=np.invert(self.data['UL']) - msk=msk1*msk2 + print ("---> bin_width=",bin_width) + + for id in range(bin_grid.size-1): + if id=bin_grid[id]) + else: + msk_bin=np.logical_and(self.data['nu_data']<=bin_grid[id+1],self.data['nu_data']>=bin_grid[id]) - if msk.any()==True>0: - sample_size=len(self.data['dnuFnu_data'][msk]) - if sample_size>1: + if msk_bin.sum()>0: + msk_UL=np.invert(self.data['UL']) + msk=np.logical_and(msk_bin,msk_UL) + if msk.sum()>1: + sample_size=len(self.data['dnuFnu_data'][msk]) sigma_2_i=self.data['dnuFnu_data'][msk]**2 w=1.0/sigma_2_i @@ -963,45 +960,57 @@ def group_data(self,N_bin=None,bin_width=None,correct_dispersions=True): corr_term=np.sum(w * (self.data['nuFnu_data'][msk] - y_bin[id]) * (self.data['nuFnu_data'][msk] - y_bin[id]))/(sample_size-1) else: corr_term=1 - - + dy_bin[id]=np.sqrt(sigma_y_2_bar*corr_term) - dx_bin[id]=bin_width/2 - x_bin[id]=bin_grid[id]-dx_bin[id] - #print"x", x_bin[id],len(w) - else: - #print "xxx" - x_bin[id]=self.data['nu_data_log'][msk][0] + + elif msk.sum()==1 : y_bin[id]=self.data['nuFnu_data'][msk][0] - dx_bin[id]=bin_width/2 dy_bin[id]=self.data['dnuFnu_data'][msk][0] + elif np.invert(msk_UL).sum()>0: + x_UL=np.append(x_UL,self.data['nu_data'][ np.invert(msk_UL)]) + dx_UL=np.append(dx_UL,self.data['dnu_data'][ np.invert(msk_UL)]) + y_UL=np.append(y_UL,self.data['nuFnu_data'][ np.invert(msk_UL)]) + dy_UL=np.append(dy_UL,self.data['dnuFnu_data'][ np.invert(msk_UL)]) + + else: + y_bin[id]=-1. - self.data_reb['nu_data_log']=x_bin + + + self.data_reb=self._build_empty_table(n_rows=x_bin.size) + self.data_reb['nu_data']=x_bin + self.data_reb['dnu_data']=dx_bin self.data_reb['nuFnu_data']=y_bin - self.data_reb['dnu_data_log']=dx_bin self.data_reb['dnuFnu_data']=dy_bin + if x_UL.size>0: + self.data_reb_UL=self._build_empty_table(n_rows=x_UL.size) + self.data_reb_UL['nu_data']=x_UL + self.data_reb_UL['dnu_data']=dx_UL + self.data_reb_UL['nuFnu_data']=y_UL + self.data_reb_UL['dnuFnu_data']=dy_UL + self.data_reb_UL['UL']=True + self.data_reb=vstack([self.data_reb,self.data_reb_UL]) + + #remove empty bins - msk=self.data_reb['nu_data_log']!=0 - #print('msk',msk) + msk=self.data_reb['nuFnu_data']>0 self.data_reb=self.data_reb[msk] - - - #self.data_reb['nuFnu_data'],self.data_reb['dnuFnu_data']=self.log_to_lin(log_val=self.data_reb['nuFnu_data_log'], log_err=self.data_reb['dnuFnu_data_log']) - - self.data_reb['nu_data'],self.data_reb['dnu_data']=self.log_to_lin(log_val=self.data_reb['nu_data_log'], log_err=self.data_reb['dnu_data_log']) - self.data_reb['nuFnu_data_log'],self.data_reb['dnuFnu_data_log']=self.lin_to_log(val=self.data_reb['nuFnu_data'], err=self.data_reb['dnuFnu_data']) - - #self.data_reb['nu_data_log'],self.data_reb['dnu_data']=self.lin_to_log(val=self.data_reb['nu_data'], err=self.data_reb['dnu_data']) + self.data_reb['nu_data_log'],self.data_reb['dnu_data_log']=self.lin_to_log(val=self.data_reb['nu_data'], err=self.data_reb['dnu_data']) #set original units for c in self.data.colnames: if self.data[c].unit is not None: self.data_reb[c] = self.data_reb[c]*self.data[c].unit - self.data=self.data_reb + if nu_min is None and nu_max is None: + self.data=self.data_reb + else: + msk=np.logical_and(self.data['nu_data']>=xmin,self.data['nu_data']<=xmax) + self.data.remove_rows(msk) + self.data=vstack([self.data,self.data_reb]) self.set_fake_error(self.fake_error) @@ -1171,11 +1180,11 @@ def plot_sed(self,plot_obj=None,frame='obs',color=None,fmt='o',ms=4,mew=0.5,figs if show_dataset is False: - plot_obj.add_data_plot(self, color=color, fmt=fmt, ms=ms, mew=mew, density=density) + plot_obj.add_data_plot(self, color=color, fmt=fmt, ms=ms, mew=mew) else: for ds in self.get_data_sets(): self.filter_data_set(filters=ds,silent=True,exclude=False) - plot_obj.add_data_plot(self, color=color, fmt=fmt, ms=ms, mew=mew ,label='dataset %s'%ds, density=density) + plot_obj.add_data_plot(self, color=color, fmt=fmt, ms=ms, mew=mew ,label='dataset %s'%ds) self.reset_data() self.reset_data() diff --git a/jetset/jet_emitters.py b/jetset/jet_emitters.py index f2c18d12..d6db10aa 100644 --- a/jetset/jet_emitters.py +++ b/jetset/jet_emitters.py @@ -761,7 +761,8 @@ def _set_L_inj(self, L_inj_target_erg, volume): pass #print('==>',self.eval_U()*volume*delta_t) self.update() - + + #NOTE: not used, to be removed def set_temp_ev(self): self.e_gamma_ptr = getattr(self._temp_ev, self._gammae_name) self._Q_inj_e_second_ptr = getattr(self._temp_ev._blob, self._Q_inj_e_second_name) diff --git a/jetset/jet_emitters_factory.py b/jetset/jet_emitters_factory.py index c137f178..6b7cf153 100644 --- a/jetset/jet_emitters_factory.py +++ b/jetset/jet_emitters_factory.py @@ -9,6 +9,7 @@ 'pl': 'powerlaw', 'lppl': 'log-parabola with low-energy powerlaw branch', 'lpep': 'log-parabola defined by peak energy', + 'lppl_pileup': 'log-parabola with low-energy powerlaw branch and pileup', 'plc': 'powerlaw with cut-off', 'bkn': 'broken powerlaw', 'superexp': 'powerlaw with super-exp cut-off'} @@ -48,6 +49,28 @@ def distr_func_lppl(gamma, gamma0_log_parab, r, s): return f +def distr_func_lppl_pileup(gamma, gamma0_log_parab, gamma_inj, r, s, gamma_eq, ratio_pile_up , alpha, gamma_cut_acc): + b = np.zeros(gamma.shape) + a = np.zeros(gamma.shape) + m = gamma < gamma_inj + s1=s+0.5 + b[m] = np.power(gamma[m]/gamma0_log_parab, s1) + + + f1 = np.zeros(gamma.shape) + m1 = np.logical_and(gamma < gamma0_log_parab,gamma>gamma_inj) + f1[m1] = np.power(gamma[m1]/gamma0_log_parab, -s) + f1[~m1] = np.power(gamma[~m1]/gamma0_log_parab, (-s -r*np.log10(gamma[~m1] / gamma0_log_parab ))) + b[~m] = np.power(gamma_inj/gamma0_log_parab,s1+s)*f1[~m]*np.exp(-np.power((gamma[~m]/gamma_cut_acc),alpha)/alpha) + + pile_up=gamma*gamma*np.exp(-np.power((gamma/gamma_eq),alpha)/alpha) + pile_up=pile_up/(gamma_eq*gamma_eq*np.exp(-np.power((gamma_eq/gamma_eq),alpha)/alpha)) + g_ID=np.argmin(np.fabs(gamma-gamma_eq)) + c=b[g_ID]*ratio_pile_up + a=pile_up*c + return a+b + + class EmittersFactory: def __repr__(self): return str(pprint.pprint(self._available_dict)) @@ -60,6 +83,7 @@ def __init__(self): 'plc': self._create_plc, 'lp': self._create_lp, 'lppl': self._create_lppl, + 'lppl_pileup': self._create_lppl_pileup, 'lpep': self._create_lpep, 'superexp': self._create_super_exp} @@ -237,6 +261,33 @@ def _create_lppl(self, gamma_grid_size, log_values, normalize, skip_build, emitt return n_lppl + def _create_lppl_pileup(self, gamma_grid_size, log_values, normalize, skip_build, emitters_type): + n_lppl_pileup = self._emitters_class(name='lppl_pileup', + spectral_type='lp', + normalize=normalize, + emitters_type=emitters_type, + log_values=log_values, + skip_build=skip_build, + gamma_grid_size=gamma_grid_size) + + + a_t, b_t = n_lppl_pileup.set_bounds(1, 1E9, log_val=n_lppl_pileup._log_values) + gamma0_log_parab_val = n_lppl_pileup._set_log_val(1E4,log_val=log_values) + + n_lppl_pileup.add_par('gamma0_log_parab', par_type='turn-over-energy', val=gamma0_log_parab_val, vmin=a_t, vmax=b_t, + unit='lorentz-factor',log=log_values) + n_lppl_pileup.add_par('s', par_type='LE_spectral_slope', val=2.0, vmin=-10., vmax=10, unit='') + n_lppl_pileup.add_par('r', par_type='spectral_curvature', val=0.4, vmin=-15., vmax=15., unit='') + n_lppl_pileup.add_par('alpha', par_type='spectral_curvature', val=0.4, vmin=1E-3, vmax=2., unit='') + n_lppl_pileup.add_par('gamma_inj', par_type='low-energy-cut-off', val=1E2, vmin=1, vmax=1E6, unit='') + n_lppl_pileup.add_par('gamma_eq', par_type='turn-over-energy', val=1E5, vmin=1E1, vmax=1E9, unit='') + n_lppl_pileup.add_par('gamma_cut_acc', par_type='high-energy-cut-off', val=1E6, vmin=1E1, vmax=1E9, unit='') + n_lppl_pileup.add_par('ratio_pile_up', par_type='scaling_factor', val=1E-3, vmin=1E-10, vmax=1E3, unit='') + + n_lppl_pileup.set_distr_func(distr_func_lppl_pileup) + + return n_lppl_pileup + class InjEmittersFactory(EmittersFactory): def __init__(self): diff --git a/jetset/jet_model.py b/jetset/jet_model.py index dc93a6b8..c5447c7e 100644 --- a/jetset/jet_model.py +++ b/jetset/jet_model.py @@ -267,13 +267,13 @@ def _serialize_model(self): _model['internal_pars']['Norm_distr'] = self.Norm_distr return _model - def save_model(self,file_name): - pickle.dump(self, open(file_name, 'wb'), protocol=pickle.HIGHEST_PROTOCOL) + #def save_model(self,file_name,to_string=False): + # pickle.dump(self, open(file_name, 'wb'), protocol=pickle.HIGHEST_PROTOCOL) @classmethod @safe_run - def load_model(cls, file_name): + def load_model(cls, file_name_or_obj, from_string=False): """Load a save model Parameters @@ -294,10 +294,10 @@ def load_model(cls, file_name): _description_ """ try: - jet=pickle.load(open(file_name, "rb")) - jet.set_blob() - jet._update_spectral_components() - return jet + jet=cls._load_pickle(file_name_or_obj,from_string=from_string) + jet.set_blob() + jet._update_spectral_components() + return jet except Exception as e: raise RuntimeError('The model you loaded is not valid please check the file name', e) @@ -1437,7 +1437,7 @@ def plot_model(self,plot_obj=None,clean=False,label=None,comp=None,sed_data=None else: comp_label=None if c.state!='off': - plot_obj.add_model_plot(c.SED, line_style=line_style, label=comp_label,flim=self.flux_plot_lim,color=color,auto_label=auto_label, density=density,update=False, frame=frame) + plot_obj.add_model_plot(c.SED, line_style=line_style, label=comp_label,flim=self.flux_plot_lim,color=color,auto_label=auto_label,update=False, frame=frame) else: for c in self.spectral_components_list: @@ -1446,7 +1446,7 @@ def plot_model(self,plot_obj=None,clean=False,label=None,comp=None,sed_data=None comp_label=label #print('comp label',comp_label) if c.state != 'off' and c.name!='Sum': - plot_obj.add_model_plot(c.SED, line_style=line_style, label=comp_label,flim=self.flux_plot_lim,auto_label=auto_label,color=color, density=density,update=False, frame=frame) + plot_obj.add_model_plot(c.SED, line_style=line_style, label=comp_label,flim=self.flux_plot_lim,auto_label=auto_label,color=color,update=False, frame=frame) c=self.get_spectral_component_by_name('Sum') if label is not None: @@ -1454,7 +1454,7 @@ def plot_model(self,plot_obj=None,clean=False,label=None,comp=None,sed_data=None else: comp_label='Sum' - plot_obj.add_model_plot(c.SED, line_style='--', label=comp_label, flim=self.flux_plot_lim,color=color, density=density,update=False) + plot_obj.add_model_plot(c.SED, line_style='--', label=comp_label, flim=self.flux_plot_lim,color=color,update=False) plot_obj.update_plot() @@ -2113,7 +2113,7 @@ def set_B_eq(self, nuFnu_obs, nu_obs, B_min=1E-9,B_max=1.0,N_pts=20,plot=False): print ('setting N to ',N[ID_min]) return b_grid[ID_min],b_grid,U_B,U_e - def eval_synch_pol(self,nu_range): + def eval_synch_pol(self,nu_range_obs): """_summary_ evluates the synchrotron polarization for Parameters @@ -2126,10 +2126,12 @@ def eval_synch_pol(self,nu_range): ------- arrrays of polarization and nuF_nu """ - nuF_nu=self.eval(get_model=True,nu=nu_range) + nuF_nu=self.eval(get_model=True,nu=nu_range_obs) - pol_nu=np.zeros(nu_range.size) - for ID,nu in enumerate(nu_range): + pol_nu=np.zeros(nu_range_obs.size) + #TODO: this will be removed when eval_Sync_polarization will follow the same pattern of synch flux + nu_range_pol_blob=nu_range_obs/self.get_beaming()*(1+self.parameters.z_cosm.val) + for ID,nu in enumerate(nu_range_pol_blob): pol_nu[ID]=BlazarSED.eval_Sync_polarization(self._blob,nu) m=np.logical_or(nuF_nu<=0,np.isnan(pol_nu)) diff --git a/jetset/mcmc.py b/jetset/mcmc.py index 3ffb76d3..4f940643 100644 --- a/jetset/mcmc.py +++ b/jetset/mcmc.py @@ -16,7 +16,7 @@ import warnings import time import copy -import threading +from collections import OrderedDict from .plot_sedfit import plt, PlotSED, set_mpl __all__=['McmcSampler'] @@ -61,15 +61,16 @@ class McmcSampler(object): def __init__(self,model_minimizer): if emcee.__version__ < "3": raise RuntimeError('Please update to emcee v>=3.0.0') - - self.model = model_minimizer.fit_model + #self.model_minimizer + self.model = model_minimizer.fit_model.clone() self.data = model_minimizer.data - self.fit_par_free = model_minimizer.fit_par_free - - self._progress_iter = cycle(['|', '/', '-', '\\']) - self.labels=None - self.par_array=None + #self._fit_par_free = self.model_minimizer.fit_par_free + self._par_array=None self._bounds=None + + + self._progress_iter = cycle(['|', '/', '-', '\\']) + def set_labels(self,use_labels_dict=None): @@ -80,33 +81,132 @@ def set_labels(self,use_labels_dict=None): use_labels_dict : _type_, optional _description_, by default None """ - if use_labels_dict is None: - self.par_array = self.fit_par_free - self.labels = [par.name for par in self.par_array] - self.labels_units = [par.units for par in self.par_array] - self.labels_start_val = [p.best_fit_val for p in self.par_array] + self._par_array=[] + self._bounds=None + if use_labels_dict is None: + for comp in self.model.components._components_list: + mod = getattr(self.model,comp.name) + for p in mod.parameters.par_array: + if p.frozen == False: + self._append_label_dict(comp.name,p) else: - self.labels=[] - self.par_array=[] - self.labels_units=[] - self.labels_start_val=[] for model_name in use_labels_dict.keys(): for par_name in use_labels_dict[model_name]: p= self.model.parameters.get_par_by_name(model_name,par_name) if p is not None: - self.par_array.append(p) - self.labels.append( p.name) - self.labels_units.append(p.units) - self.labels_start_val.append(p.best_fit_val) - else: - warnings.warn('par %s'%par_name+' not present in model, will be skipped') + if p.frozen == False: + self._append_label_dict(model_name,p) + + def _append_label_dict(self,comp_name,p): + self._par_array.append(OrderedDict(comp_name=comp_name, + name=p.name, + full_name=f"{comp_name}.{p.name}", + val=p.val, + val_mix=p.val_min, + val_max=p.val_max, + labels_start_val=p.val, + minimizer_best_fit_val=p.best_fit_val, + minimizer_best_fit_err=p.best_fit_err, + minimizer_fit_range_min=p.fit_range_min, + minimizer_fit_range_max=p.fit_range_max, + units=p.units, + plot_label=p.name, + bounds=[None,None])) + + + def set_bounds(self,bound=0.2,bound_rel=False,preserve_fit_range=True): + self._set_bounds(bound=bound,bound_rel=bound_rel,preserve_fit_range=preserve_fit_range) + + def _set_bounds(self, bound=0.2,bound_rel=True,preserve_fit_range=True): + + self._bounds=[] + + if np.shape(bound) == (): + bound=[bound,bound] + elif np.shape(bound) == (2,): + pass + else: + raise RuntimeError('bound shape', np.shape(bound), 'it is wrong, has to be a scalar or (2,)') + for par in self._par_array: + if not bound_rel and par['best_fit_err'] is not None: + delta_p = par['minimizer_best_fit_err'] * bound[1] + delta_m = par['minimizer_best_fit_err'] * bound[0] - def set_bounds(self,bound=0.2,bound_rel=False,): + else: + delta_p = np.fabs(par['minimizer_best_fit_val'])*bound[1] + delta_m = np.fabs(par['minimizer_best_fit_val'])*bound[0] + + _min = par['minimizer_best_fit_val'] - delta_m + _max = par['minimizer_best_fit_val'] + delta_p + + + if par['minimizer_fit_range_min'] is not None and preserve_fit_range is True: + _min= max(_min, par['minimizer_fit_range_min'] ) + elif par['val_min'] is not None: + _min= max(_min, par['val_min']) + + if par['minimizer_fit_range_max'] is not None: + _max= min(_max, par['minimizer_fit_range_max']) + elif par['val_max'] is not None: + _max= min(_max, par['val_max']) + + print('par:',par['name'],' best fit value: ',par['minimizer_best_fit_val'],' mcmc bounds:',[_min, _max]) + par['bounds']=[_min, _max] + + def _build_sampler_bounds(self): + self._bounds=[] + for par in self._par_array: + self._bounds.append(par['bounds']) + + def get_par(self, par_name_or_idx, comp_name=None, get_index=False): + if type(par_name_or_idx) == int: + par_idx=par_name_or_idx + + else: + par_name=par_name_or_idx + try: + if comp_name is None: + par_idx = [par['name'] for par in self._par_array].index(par_name) + else: + par_idx = [par['name']+par['comp_name'] for par in self._par_array].index(par_name+comp_name) + except: + raise RuntimeError('parameter p', par_name, 'not found') - self._build_bounds(bound=bound,bound_rel=bound_rel) - + if par_idx > len(self._par_array): + raise RuntimeError('label id larger then labels size') + + if get_index is True: + return self._par_array[par_idx], par_idx + else: + return self._par_array[par_idx] + + + + def set_plot_label(self,par_name,plot_label,comp_name=None): + p=self.get_par(par_name,comp_name=comp_name) + p['plot_label']=plot_label + + def reset_to_minimizer_best_fit(self): + for par in self._par_array: + par['val'] = par['minimizer_best_fit_val'] + + #def reset_to_mcmc(self,quantile=0.5): + # for ID,par in enumerate(self.par_array): + # q_vals=self.get_par_quantiles(ID,quantiles=quantile) + # par.val = q_vals + + def reset_to_mcmc_best_fit(self): + _prob_max = np.argmax(self.samples_log_prob) + _id_prob_max = np.unravel_index(_prob_max, self.samples_log_prob.shape) + + print("----------------------------") + print("MCMC best fit solution") + for ID,par in enumerate(self._par_array): + par['val'] = self.get_sample(ID)[_id_prob_max] + print(f"{par['name']}: {par['val']}") + print("----------------------------") def run_sampler(self, @@ -118,37 +218,30 @@ def run_sampler(self, threads=None, walker_start_bound=0.005, loglog = False, - progress='notebook', - use_labels_dict=None, - bound=None, - bound_rel=False): + progress='notebook'): + + + if self._par_array is None: + raise RuntimeError('please set the labels, using .set_labels, before running the sampler') + + self._build_sampler_bounds() - self.calls = 0 self.calls_OK = 0 self.use_UL = use_UL - - - if use_labels_dict is not None or bound is not None: - warnings.warn('use_labels_dict, and bounds in run_sampler are deprecated and will result in an error in the next version, please use the .set_labels and .set_bounds methods as explained in the documentation') - self.set_labels(use_labels_dict=use_labels_dict) - self.set_bounds(bound=bound,bound_rel=bound_rel) - - - self.plot_labels=copy.deepcopy(self.labels) - self.par_array_best_fit=copy.deepcopy(self.par_array) - self.ndim = len(self.labels) + + self.ndim = len(self._par_array) self.pos = pos self.burnin=burnin if nwalkers is None: - self.nwalkers = 4*len(self.labels) + self.nwalkers = 4*len(self._par_array) print('setting nwalkers to:', self.nwalkers) else: self.nwalkers = nwalkers - if self.nwalkers < 2*len(self.labels): + if self.nwalkers < 2*len(self._par_array): raise RuntimeError("numbers of walkers has to be at least two times the number of sampling pars") counter=Counter( self.nwalkers*steps) @@ -156,24 +249,22 @@ def run_sampler(self, self.calls_tot = self.nwalkers * steps if pos is None: - pos = sample_ball(np.array([p.best_fit_val for p in self.par_array]), - np.array([p.best_fit_val * walker_start_bound for p in self.par_array]), + pos = sample_ball(np.array([p['minimizer_best_fit_val'] for p in self._par_array]), + np.array([p['minimizer_best_fit_val'] * walker_start_bound for p in self._par_array]), self.nwalkers) print('mcmc run starting') print('') start = time.time() - if threads is not None and threads>1: - warnings.warn('python multithreading is not effective, JetSeT uses C threads to speedup computation') threads=1 - self.sampler = emcee.EnsembleSampler(self.nwalkers, self.ndim, log_prob, args=(self.model, self.data, use_UL, counter, self._bounds, self.par_array, loglog)) + self.sampler = emcee.EnsembleSampler(self.nwalkers, self.ndim, log_prob, args=(self.model, self.data, use_UL, counter, self._bounds, self._par_array, loglog)) self.sampler.run_mcmc(pos, steps, progress=progress) else: threads=1 - self.sampler = emcee.EnsembleSampler(self.nwalkers, self.ndim, log_prob, args=(self.model, self.data, use_UL, counter, self._bounds, self.par_array, loglog)) + self.sampler = emcee.EnsembleSampler(self.nwalkers, self.ndim, log_prob, args=(self.model, self.data, use_UL, counter, self._bounds, self._par_array, loglog)) self.sampler.run_mcmc(pos, steps, progress=progress) @@ -182,145 +273,125 @@ def run_sampler(self, comp_time = end - start print("mcmc run done, with %d threads took %2.2f seconds"%(threads,comp_time)) - #self.samples = self.sampler.chain[:, burnin:, :].reshape((-1, self.ndim)) + self.samples = self.sampler.get_chain(flat=True,discard=burnin) self.samples_log_prob = self.sampler.get_log_prob(flat=True,discard=burnin) self.acceptance_fraction=np.mean(self.sampler.acceptance_fraction) + self.reset_to_minimizer_best_fit() - self.reset_to_best_fit() - - - def set_plot_label(self,label,plot_label): - self.plot_labels[self.labels.index(label)]=plot_label - - - def get_par_quantiles(self,p,quantiles=(0.16,0.5,0.84)): - _d, idx=self.get_par(p) - return np.array(np.quantile(_d,quantiles)) + def get_par_quantiles(self,par_name,comp_name=None,quantiles=(0.16,0.5,0.84)): + return np.array(np.quantile(self.get_sample(par_name,comp_name=comp_name),quantiles)) - def reset_to_best_fit(self): - for ID,par in enumerate(self.par_array): - par.val = self.par_array_best_fit[ID].val - - - def reset_to_mcmc(self,quantile=0.5): - for ID,par in enumerate(self.par_array): - q_vals=self.get_par_quantiles(ID,quantiles=quantile) - par.val = q_vals - - def _build_bounds(self, bound=0.2,bound_rel=True): + - self._bounds=[] + def corner_plot(self, labels = None, comp_name=None,quantiles = (0.16, 0.5, 0.84), levels = None, title_kwargs = {}, **kwargs): + """_summary_ - if np.shape(bound) == (): - bound=[bound,bound] - elif np.shape(bound) == (2,): - pass + Parameters + ---------- + labels : _type_, optional + _description_, by default None + quantiles : tuple, optional + _description_, by default (0.16, 0.5, 0.84) + levels :levels=(1 - np.exp(-0.5),) + check for proper levels definition for 2d: https://corner.readthedocs.io/en/latest/pages/sigmas/ + title_kwargs : dict, optional + _description_, by default {} + + Returns + ------- + _type_ + _description_ + """ + + if comp_name is None: + components=np.unique([p['comp_name'] for p in self._par_array]) else: - raise RuntimeError('bound shape', np.shape(bound), 'it is wrong, has to be a scalar or (2,)') - - for par in self.par_array: - if bound_rel is False and par.best_fit_err is not None: - - #_min = par.best_fit_val - par.best_fit_err * bound[0] - #_max = par.best_fit_val + par.best_fit_err * bound[1] - delta_p = par.best_fit_err * bound[1] - delta_m = par.best_fit_err * bound[0] - - else: - delta_p = np.fabs(par.best_fit_val)*bound[1] - delta_m = np.fabs(par.best_fit_val)*bound[0] - #_min = par.best_fit_val * (1.0 - bound[0]) - #_max = par.best_fit_val * (1.0 + bound[1]) - - _min = par.best_fit_val - delta_m - _max = par.best_fit_val + delta_p - - - if par.fit_range_min is not None: - _min= max(_min, par.fit_range_min) - - if par.fit_range_max is not None: - _max= min(_max, par.fit_range_max) + components=[comp_name] + f_list=[] + for c in components: + _idxs = [] + truths = [] - print('par:',par.name,' best fit value: ',par.best_fit_val,' mcmc bounds:',[_min, _max]) - self._bounds.append([_min, _max]) - - - - def corner_plot(self, labels = None, quantiles = (0.16, 0.5, 0.84), levels = None, title_kwargs = {}, **kwargs): - _id = [] - - if labels is None: - labels=self.labels - - if type(labels) == list: - pass - else: - labels = [labels] - - for l in labels: - _id.append(self.labels.index(l)) - - f = corner.corner(self.samples[:, _id], - quantiles=quantiles, labels=self.plot_labels, - truths=[self.labels_start_val[i] for i in _id], - title_kwargs=title_kwargs,show_titles = True, - levels = levels,**kwargs) - - title = 'quantiles ='+str(quantiles) - f.suptitle(title,y=1.0) - return f - - def get_par(self, p,): - if type(p) == int: - pass - - else: - try: - p = self.labels.index(p) - except: - raise RuntimeError('parameter p', p, 'not found') + + msk=np.array([p['comp_name']==c for p in self._par_array]) - if p > len(self.labels): - raise RuntimeError('label id larger then labels size') + if labels is None: + plot_labels= [p['name'] for p in np.array(self._par_array)[msk]] + + elif type(labels) == list: + plot_labels=labels + else: + plot_labels = [labels] - return self.samples[:, p].flatten(), p + if len(plot_labels)>0: + for l in plot_labels: + + _idxs.append(self.get_par(l,comp_name=c,get_index=True)[1]) + + for _idx in _idxs: + truths.append(self.get_par(_idx)['minimizer_best_fit_val']) + plot_labels.append(self.get_par(_idx)['plot_label']) + + f = corner.corner(self.samples[:, _idxs], + quantiles=quantiles, + labels=plot_labels, + truths=truths, + title_kwargs=title_kwargs, + show_titles = True, + levels = levels,**kwargs) + + + #print(c,str(quantiles)) + title = c + ' quantiles ='+str(quantiles) + + f.suptitle(title,y=1.0) + f_list.append(f) + return f_list - def plot_chain(self,p=None,log_plot=False): - if p is None: - p=self.labels + + def plot_chain(self,par_name=None, comp_name=None,log_plot=False): + f_list=[] + if comp_name is None: + components=np.unique([p['comp_name'] for p in self._par_array]) else: - p=np.atleast_1d(p) - - f, axes = plt.subplots(len(p), sharex=True) - axes=np.atleast_1d(axes) - for ID,_p in enumerate(p): - self._plot_chain(_p,axes[ID],log_plot=log_plot) - - axes[-1].set_xlabel('steps') - return f + components=[comp_name] + for c in components: + msk=np.array([p['comp_name']==c for p in self._par_array]) + if par_name is None: + + par_names=[par['name'] for par in np.array(self._par_array)[msk]] + else: + par_names=np.atleast_1d(par_name) - def _plot_chain(self, p,ax, log_plot=False): - _d, idx = self.get_par(p) + f, axes = plt.subplots(len(par_names), sharex=True) + axes=np.atleast_1d(axes) + for ID,_p_name in enumerate(par_names): + self._plot_chain(_p_name,axes[ID],comp_name=comp_name,log_plot=log_plot) + + axes[-1].set_xlabel('steps') + f_list.append(f) + + return f_list - n = self.plot_labels[idx] + def _plot_chain(self, par_name,ax,comp_name=None, log_plot=False): + par = self.get_par(par_name,comp_name=comp_name) - traces=self.sampler.chain[:, :, idx] + n = par['plot_label'] - if self.labels_units is not None: - if self.labels_units[idx] is not None: - n += ' (%s)' % self.labels_units[idx] + traces=self.get_trace(par_name,comp_name=comp_name) - alpha_true = np.median(_d) + if par['units'] is not None: + n += ' (%s)' % par['units'] - + _s=self.get_sample(par_name,comp_name=comp_name) + alpha_true = np.median(_s) if log_plot == True: n = 'log10(%s)'%n - if np.any(_d <= 0): + if np.any(_s <= 0): raise RuntimeWarning('negative values in p') else: traces = np.log10(traces) @@ -336,23 +407,30 @@ def _plot_chain(self, p,ax, log_plot=False): ax.set_ylabel(n) - + + + def get_trace(self, par_name,comp_name=None): + _p,p_idx=self.get_par(par_name,comp_name=comp_name,get_index=True) + return self.sampler.chain[:, :, p_idx] - def plot_par(self, p, nbins=20, log_plot=False,quantiles=(0.16,0.5,0.84),figsize=None): + + def plot_par(self, par_name, comp_name=None,nbins=20, log_plot=False,quantiles=(0.16,0.5,0.84),figsize=None): set_mpl() - _d, idx = self.get_par(p) + par = self.get_par(par_name,comp_name=comp_name) - par_name = self.plot_labels[idx] + par_name = par['name'] x_name = par_name - if self.labels_units is not None: - if self.labels_units[idx] is not None and str(self.labels_units[idx]).strip() != '' : - x_name += ' (%s)' % self.labels_units[idx] + if par['units'] is not None: + x_name += ' (%s)' % par['units'] + + _d=self.get_sample(par_name,comp_name=comp_name) + f = plt.figure(figsize=figsize) ax = f.add_subplot(111) - + if log_plot == True: x_name = 'log10(%s)' % x_name if np.any(_d <= 0): @@ -360,7 +438,7 @@ def plot_par(self, p, nbins=20, log_plot=False,quantiles=(0.16,0.5,0.84),figsize else: _d = np.log10(_d) - q_vals=self.get_par_quantiles(p,quantiles=quantiles) + q_vals=self.get_par_quantiles(par_name,comp_name=comp_name,quantiles=quantiles) q_diff = np.diff(q_vals) @@ -379,9 +457,40 @@ def plot_par(self, p, nbins=20, log_plot=False,quantiles=(0.16,0.5,0.84),figsize f.suptitle('quantiles = %s'%str(quantiles)) ax.legend(loc='center left', bbox_to_anchor=(1.0, 0.5), ncol=1) return f + + def get_sample(self, par_name,comp_name=None): + _p,p_idx=self.get_par(par_name,comp_name=comp_name,get_index=True) + return self.samples[:,p_idx] - def plot_model(self, sed_data=None, fit_range=None, size=100, frame='obs', density=False,quantiles=None, get_model=False, plot_mcmc_best_fit_model=False): + def plot_model(self, sed_data=None, fit_range=None, size=100, frame='obs', density=False,quantiles=None, get_model=False, plot_mcmc_best_fit_model=False,rnd_seed=0): + """_summary_ + Parameters + ---------- + sed_data : _type_, optional + _description_, by default None + fit_range : _type_, optional + _description_, by default None + size : int, optional + _description_, by default 100 + frame : str, optional + _description_, by default 'obs' + density : bool, optional + _description_, by default False + quantiles : _type_, optional + _description_, by default None + get_model : bool, optional + _description_, by default False + plot_mcmc_best_fit_model : bool, optional + _description_, by default False + rnd_seed : int, optional + _description_, by default 0 + + Returns + ------- + _type_ + _description_ + """ if sed_data is None: sed_data=self.sed_data @@ -389,55 +498,36 @@ def plot_model(self, sed_data=None, fit_range=None, size=100, frame='obs', densi fit_range = [self.model.nu_min_fit, self.model.nu_max_fit] p = self.model._set_up_plot(None, sed_data, frame, density) - - self.reset_to_best_fit() - self.model.eval(fill_SED=True) - - x, y = self.model.SED.get_model_points(log_log=False, frame=frame) - #if density is True: - # y=y-x - if size is None: - size = len(self.samples) - ID_mcmc = np.arange(size) - else: - size = min(len(self.samples), int(size)) - ID_mcmc = np.random.randint(len(self.samples), size=size) - - y = np.zeros((size,x.size)) - - for ID,ID_rand in enumerate(ID_mcmc): - - for ID_par, pi in enumerate(self.par_array): - pi.set(val=self.get_par(ID_par)[0][ID_rand]) - - self.model.eval(fill_SED=True) - x, y[ID] = self.model.SED.get_model_points(log_log=False, frame=frame) - + x,y=self._get_model_samples(size=size,rnd_seed=rnd_seed,frame=frame) + if density is True: y=y/x if quantiles is None: y_min=np.amin(y, axis=0) y_max=np.amax(y, axis=0) - msk = y_min > self.model.flux_plot_lim - l=p.sedplot.fill_between(x[msk],y_max[msk],y_min[msk],color='gray',alpha=0.3,label='mcmc model range') + _l='mcmc model range' + #msk = y_min > self.model.flux_plot_lim + #l=p.sedplot.fill_between(x[msk],y_max[msk],y_min[msk],color='gray',alpha=0.3,label='mcmc model range') else: + _l='mcmc model conf. %s'%quantiles y_min,y_max=np.quantile(y, quantiles, axis=0) - msk = y_min > self.model.flux_plot_lim - l=p.sedplot.fill_between(x[msk],y_max[msk],y_min[msk],color='gray',alpha=0.3,label='mcmc model conf %s'%quantiles) - p.lines_model_list.append(l) - self.reset_to_best_fit() - self.model.eval(fill_SED=True) + + msk = y_min > self.model.flux_plot_lim + l=p.sedplot.fill_between(x[msk],y_max[msk],y_min[msk],color='gray',alpha=0.3,label=_l) + p.lines_model_list.append(l) + msk = y_min > self.model.flux_plot_lim if plot_mcmc_best_fit_model is False: - p.add_model_plot(self.model, color='red',fit_range = fit_range,density=density,flim=self.model.flux_plot_lim) - p.add_model_residual_plot(model = self.model, data = sed_data, fit_range = fit_range, color='red') + self.reset_to_minimizer_best_fit() + label=None else: - self.reset_to_mcmc(quantile=0.5) - self.model.eval() - p.add_model_plot(self.model, color='red',fit_range = fit_range,density=density,flim=self.model.flux_plot_lim,label='mcmc 0.5 quantile') - p.add_model_residual_plot(model = self.model, data = sed_data, fit_range = fit_range, color='red') - #self.reset_to_best_fit() + label='mcmc best fit' + self.reset_to_mcmc_best_fit() + + self.model.eval(fill_SED=True) + p.add_model_plot(self.model, color='red',fit_range = fit_range,flim=self.model.flux_plot_lim,label=label) + p.add_model_residual_plot(model = self.model, data = sed_data, fit_range = fit_range, color='red') if get_model is True: return p, [x[msk],y_min[msk],y_max[msk]] @@ -445,6 +535,30 @@ def plot_model(self, sed_data=None, fit_range=None, size=100, frame='obs', densi return p + def _get_model_samples(self,size,rnd_seed,frame): + self.model.eval() + x, _y = self.model.SED.get_model_points(log_log=False, frame=frame) + + if size is None: + size = len(self.samples) + else: + size = min(len(self.samples), int(size)) + + rng = np.random.default_rng(rnd_seed) + ID_mcmc = rng.integers(0,len(self.samples),size=size) + + y = np.zeros((size,x.size)) + + for ID,ID_rand in enumerate(ID_mcmc): + for id_p,par in enumerate(self._par_array): + par['val']=self.get_sample(id_p)[ID_rand] + self.model.parameters.set_par(model_name= par['comp_name'],par_name=par['name'],val=par['val']) + + self.model.eval(fill_SED=True) + x, y[ID] = self.model.SED.get_model_points(log_log=False, frame=frame) + + return x,y + def _progess_bar(self,): if np.mod(self.calls, 10) == 0 and self.calls != 0: print("\r%s progress=%3.3f%% calls=%d accepted=%d" % (next(self._progress_iter),float(100*self.calls)/(self.calls_tot),self.calls,self.calls_OK), end="") @@ -454,27 +568,43 @@ def save(self, name): pickle.dump(self, output, pickle.HIGHEST_PROTOCOL) @classmethod - def load(self, name): - with open(name, 'rb') as input: - return pickle.load(input) + #def load(self, name): + # with open(name, 'rb') as input: + # return pickle.load(input) + def load(cls, file_name): + + try: + c = pickle.load(open(file_name, "rb")) + if isinstance(c, McmcSampler): + #c.__init__(c.minimizer) + if hasattr(c,'model'): + c.model=c.model._build_model(c.model) + return c + else: + raise RuntimeError('The model you loaded is not valid please check the file name') + + except Exception as e: + raise RuntimeError(e) def emcee_log_like(theta,fit_model,data,use_UL,par_array,loglog): _warn = False for pi in range(len(theta)): - par_array[pi].set(val=theta[pi]) + + + par_array[pi]['val']=theta[pi] + fit_model.parameters.set_par(model_name= par_array[pi]['comp_name'],par_name=par_array[pi]['name'],val=par_array[pi]['val']) if np.isnan(theta[pi]): _warn=True - _m = fit_model.eval(nu=data['x'], fill_SED=False, get_model=True, loglog=loglog) + _model = fit_model.eval(nu=data['x'], fill_SED=False, get_model=True, loglog=loglog) _res_sum, _res, _res_UL = _eval_res(data['y'], - _m, + _model, data['dy'], data['UL'], use_UL=use_UL) - #_progess_bar(counter) return _res_sum *-0.5 @@ -489,19 +619,19 @@ def log_prob(theta,fit_model,data,use_UL,counter,bounds,par_array,loglog): lp = log_prior(theta,bounds) counter.count += 1 if not np.isfinite(lp): - lp = -np.inf + res = -np.inf else: - lp += emcee_log_like(theta,fit_model,data,use_UL,par_array,loglog) + ll= emcee_log_like(theta,fit_model,data,use_UL,par_array,loglog) + res= lp+ll counter.count_OK += 1 - return lp + return res def log_prior(theta,bounds): _r=0. - #bounds = [(par.fit_range_min, par.fit_range_max) for par in model_minimizer.fit_par_free] - #skip=False + for pi in range(len(theta)): if bounds[pi][1] is not None: @@ -512,10 +642,6 @@ def log_prior(theta,bounds): if theta[pi] reval',mc.name) mc.eval() _c=mc.spectral_components.get_spectral_component_by_name(c.name) - plot_obj.add_model_plot(_c.SED, line_style=line_style, label=' -%s'%comp_label, flim=self.flux_plot_lim, density=density, frame=frame) + plot_obj.add_model_plot(_c.SED, line_style=line_style, label=' -%s'%comp_label, flim=self.flux_plot_lim, frame=frame) except Exception as e: raise RuntimeError('for model', mc.name, "spectral component",c.name, "problem with plotting SED", e) line_style = '-' if label is None: label=self.name - plot_obj.add_model_plot(self.SED, line_style=line_style, label=label, flim=self.flux_plot_lim,fit_range=[self.nu_min_fit,self.nu_max_fit], density=density, frame=frame ) + plot_obj.add_model_plot(self.SED, line_style=line_style, label=label, flim=self.flux_plot_lim,fit_range=[self.nu_min_fit,self.nu_max_fit], frame=frame ) plot_obj.add_model_residual_plot(data=sed_data, model=self, fit_range=[self.nu_min_fit, self.nu_max_fit]) #if frame == 'src' and sed_data is not None: @@ -351,8 +351,8 @@ def eval(self,nu=None,fill_SED=True,get_model=False,loglog=False,label=None,phys @classmethod - def load_model(cls, file_name): - c = pickle.load(open(file_name, "rb")) + def load_model(cls, file_name_or_obj, from_string=False): + c = cls._load_pickle(file_name_or_obj,from_string=from_string) return cls._build_model(c) @staticmethod @@ -392,8 +392,8 @@ def set_fit_range(self,down_tol=0.1,up_tol=100): m.set_fit_range(down_tol=down_tol,up_tol=up_tol) - def clone(self): - return self._build_model(pickle.loads(pickle.dumps(self, protocol=pickle.HIGHEST_PROTOCOL))) + #def clone(self): + # return self.load_model(pickle.loads(pickle.dumps(self, protocol=pickle.HIGHEST_PROTOCOL))) def show_model_components(self): print("") @@ -435,6 +435,7 @@ def sed_tables_dict(self, restframe='obs'): self._sed_tables_dict[self.name]=self.sed_table(restframe=restframe) for comp in self.components.components_list: if hasattr(comp,'sed_table'): - self._sed_tables_dict[comp.name]=comp.sed_table(restframe=restframe) + if comp.sed_table(restframe=restframe) is not None: + self._sed_tables_dict[comp.name]=comp.sed_table(restframe=restframe) return self._sed_tables_dict diff --git a/jetset/model_parameters.py b/jetset/model_parameters.py index 45e9d3e9..15d3ffd0 100644 --- a/jetset/model_parameters.py +++ b/jetset/model_parameters.py @@ -28,7 +28,8 @@ def _show_table(t): if is_notebook(): try: from IPython.display import display - display(t.show_in_notebook(show_row_index=False, display_length=100)) + display(t.show_in_notebook(show_row_index=False, display_length=100, backend='classic')) + #display(t.show_in_notebook(show_row_index=False, display_length=100, auto_fit_columns=True )) except Exception as e: try: from IPython.display import display @@ -321,8 +322,9 @@ def fit_range(self, fit_range=None): def reset_dependencies(self): for mp in self._master_pars: - if self in mp._depending_pars: - mp._depending_pars.remove(self) + if hasattr(mp, '_depending_pars'): + if self in mp._depending_pars: + mp._depending_pars.remove(self) self._linked = False self._linked_root_model = None self._func=None diff --git a/jetset/pkg_info.json b/jetset/pkg_info.json index 9668b1d1..ff64f562 100644 --- a/jetset/pkg_info.json +++ b/jetset/pkg_info.json @@ -1 +1 @@ -{"version": "1.3.0", "label": ""} \ No newline at end of file +{"version": "1.3.1rc2", "label": ""} \ No newline at end of file diff --git a/jetset/plot_sedfit.py b/jetset/plot_sedfit.py index b03b5e97..79150bc8 100644 --- a/jetset/plot_sedfit.py +++ b/jetset/plot_sedfit.py @@ -78,7 +78,8 @@ def __init__(self, check_frame(frame) self.frame=frame - + self._sed_data=None + self.density=density self.axis_kw=['x_min','x_max','y_min','y_max'] self.interactive=interactive @@ -113,7 +114,7 @@ def __init__(self, self.sedplot= self.fig.add_subplot(self.gs[0]) self._add_res_plot() - self.set_plot_axis_labels(density=density) + self.set_plot_axis_labels(density=self.density) #if autoscale==True: self.sedplot.set_autoscalex_on(True) @@ -146,9 +147,10 @@ def __init__(self, self.fig.canvas.manager.toolbar.update() except: pass - + + if sed_data is not None : - self.add_data_plot(sed_data,density=density) + self.add_data_plot(sed_data) if model is not None: self.add_model_plot(model) @@ -296,14 +298,22 @@ def setlim_res(self,x_min=None,x_max=None,y_min=None,y_max=None): self.update_plot() def update_plot(self): - self.fig.canvas.draw() - - y_s = [] - x_min = [] - x_max = [] - y_min = None - y_max = None - if len(self.sedplot.lines)>0: + + if self._sed_data is not None: + x,y,dx,dy,=self._sed_data.get_data_points(log_log=False,frame=self.frame, density=self.density) + + self.fig.canvas.draw() + self.sedplot.relim() + self.sedplot.set_xlim(np.min(x)/10,np.max(x)*10) + self.sedplot.set_ylim(np.min(y)/10,np.max(y)*10) + + + elif len(self.sedplot.lines)>0: + y_s = [] + x_min = [] + x_max = [] + y_min = None + y_max = None for l in self.sedplot.lines: if len(l.get_ydata())>0: @@ -353,7 +363,7 @@ def update_legend(self,label=None): - def add_model_plot(self, model, label=None, color=None, line_style=None, flim=None,auto_label=True,fit_range=None,density=False, update=True, lw=1.0 ,frame=None): + def add_model_plot(self, model, label=None, color=None, line_style=None, flim=None,auto_label=True,fit_range=None, update=True, lw=1.0 ,frame=None): frame=self._check_frame(frame=frame) @@ -368,7 +378,7 @@ def add_model_plot(self, model, label=None, color=None, line_style=None, flim=No except Exception as e: raise RuntimeError('for model',model.name, "problem with SED.get_model_points()",e) - if density is True: + if self.density is True: y=y/x if line_style is None: line_style = '-' @@ -415,7 +425,6 @@ def plot_tempev_model(self, time_slice_bin=None, time=None, time_bin=None, - density=False, use_cached=False, sed_data=None, average=False): @@ -493,7 +502,7 @@ def plot_tempev_model(self, label = 'stop, t=%2.2e (s)' % t self.add_model_plot(model=s, label=label, line_style=ls, color=color, update=False, lw=lw, - auto_label=False,density=density) + auto_label=False) if sed_data is not None: self.add_data_plot(sed_data) @@ -502,11 +511,11 @@ def plot_tempev_model(self, return - def add_data_plot(self,sed_data,label=None,color=None,frame=None,fmt='o',ms=4,mew=0.5,fit_range=None, density = False): - + def add_data_plot(self,sed_data,label=None,color=None,frame=None,fmt='o',ms=4,mew=0.5,fit_range=None): + self._sed_data=sed_data frame = self._check_frame(frame) try: - x,y,dx,dy,=sed_data.get_data_points(log_log=False,frame=self.frame, density=density) + x,y,dx,dy,=sed_data.get_data_points(log_log=False,frame=self.frame) except Exception as e: raise RuntimeError("!!! ERROR failed to get data points from", sed_data,e) diff --git a/jetset/sed_shaper.py b/jetset/sed_shaper.py index 323e44c1..2fac6383 100644 --- a/jetset/sed_shaper.py +++ b/jetset/sed_shaper.py @@ -71,7 +71,7 @@ def __init__(self,name=None,data_type=None,val=None,err=None,idx_range=[]): if name in index_names: self.name=name else: - raise RuntimeError("index name=%s not allowed, allowed names=%"%(name,index_names)) + raise RuntimeError("index name=%s not allowed, allowed names=s%"%(name,index_names)) if data_type in data_type_allowed: self.data_type=data_type @@ -89,8 +89,10 @@ def __init__(self,name=None,data_type=None,val=None,err=None,idx_range=[]): self.err=None if idx_range!=[]: + #print("===>",idx_range) self.idx_range=tuple(idx_range) else: + #print("===>",spectral_index_range(name)) self.idx_range=spectral_index_range(name) @@ -173,14 +175,14 @@ def show_pars(self): #---------------------------------------------------- def spectral_index_range(name): spectral_range_dic={} - spectral_range_dic['radio']=[6,10] - spectral_range_dic['radio_mm']=[10,11] - spectral_range_dic['mm_IR']=[11,13] - spectral_range_dic['IR_Opt']=[13,14] - spectral_range_dic['Opt_UV']=[14,16] - spectral_range_dic['UV_X']=[15,17.5] - spectral_range_dic['BBB']=[15,16] - spectral_range_dic['X']=[16,19] + spectral_range_dic['radio']=[6.,10.] + spectral_range_dic['radio_mm']=[10.,11.] + spectral_range_dic['mm_IR']=[11.,13.] + spectral_range_dic['IR_Opt']=[13.,14.] + spectral_range_dic['Opt_UV']=[14.,16.] + spectral_range_dic['UV_X']=[15.,17.5] + spectral_range_dic['BBB']=[15.,16.] + spectral_range_dic['X']=[16.,19.] spectral_range_dic['Fermi']=[22.38,25.38] spectral_range_dic['TeV'] = [25.00, 28.38] return spectral_range_dic[name] @@ -188,10 +190,10 @@ def spectral_index_range(name): def sync_fit_range(name,indices): spectral_range_dic={} - spectral_range_dic['blind']=[9,19] - spectral_range_dic['ISP']=[10,17] - spectral_range_dic['LSP']=[10,17] - spectral_range_dic['HSP']=[11,20] + spectral_range_dic['blind']=[9.,19.] + spectral_range_dic['ISP']=[10.,17.] + spectral_range_dic['LSP']=[10.,17.] + spectral_range_dic['HSP']=[11.,20.] if name=='ISP' and indices.get_by_name('X').val is not None : if indices.get_by_name('X').val.sed<0: @@ -826,10 +828,10 @@ def add_host_template(self,fit_model): fit_model.add_component(host_gal) - #print('nuFnu_p_host',self.S_peak.nuFnu_p_val) + fit_model.set('host_galaxy','nuFnu_p_host',val=(self.S_peak.nuFnu_p_val),fit_range=[self.S_peak.nuFnu_p_val-2,self.S_peak.nuFnu_p_val+2 ]) - fit_model.set('host_galaxy','nu_scale',val=0,fit_range=[-0.5,0.5]) + fit_model.set('host_galaxy','nu_scale',val=0,fit_range=[-np.log10(1+self.sed_data.z*.2),np.log10(1+self.sed_data.z*.2)]) self.host_gal=host_gal diff --git a/jetset/template_2Dmodel.py b/jetset/template_2Dmodel.py index cb296a58..c08e586e 100644 --- a/jetset/template_2Dmodel.py +++ b/jetset/template_2Dmodel.py @@ -116,7 +116,7 @@ def plot_model(self,plot_obj=None,clean=False,label=None,sed_data=None,color=Non if label is None: label=self.name - plot_obj.add_model_plot(self.SED, line_style='-', label=label, flim=self.flux_plot_lim,color=color, density=density,frame=frame) + plot_obj.add_model_plot(self.SED, line_style='-', label=label, flim=self.flux_plot_lim,color=color,frame=frame) return plot_obj diff --git a/jetset/template_model.py b/jetset/template_model.py index 0d7f63c2..7b41b7ec 100644 --- a/jetset/template_model.py +++ b/jetset/template_model.py @@ -121,7 +121,7 @@ def plot_model(self,plot_obj=None,clean=False,label=None,sed_data=None,color=Non if label is None: label=self.name - plot_obj.add_model_plot(self.SED, line_style='-', label=label, flim=self.flux_plot_lim,color=color, density=density, frame=frame) + plot_obj.add_model_plot(self.SED, line_style='-', label=label, flim=self.flux_plot_lim,color=color, frame=frame) return plot_obj diff --git a/requirements.txt b/requirements.txt index bcec690c..11c28dee 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,18 +1,18 @@ setuptools -scipy>=1.5.0,<=1.13.1 -numpy>=1.22,<2.0 -astropy>=5.0.1,<=6 +scipy>=1.5.0,<=1.15.2 +numpy==2 +astropy==7 matplotlib>=3.1.0 future -iminuit>=2.0.0,<=2.17.0 +iminuit>=2.30.0,<=2.31.2 corner six -emcee>=3.0.0,<=3.1.6 +emcee>=3.0.0,<=3.2 pyyaml pytest -numba>0.55,<=0.59.1 +numba>0.55,<=0.61.0 tqdm<5 jupyter ipython dill -swig>=4 \ No newline at end of file +swig>=4