diff --git a/gsf/fitting.py b/gsf/fitting.py index d9a368e..8089683 100644 --- a/gsf/fitting.py +++ b/gsf/fitting.py @@ -137,9 +137,9 @@ def __init__(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set:floa 'CAT_BB', 'CAT_BB_DUST', 'SNLIM', 'MORP', 'MORP_FILE', 'SPEC_FILE', 'DIR_EXTR', 'MAGZP_SPEC', 'UNIT_SPEC', 'DIFF_CONV', - 'CZ0', 'CZ1', 'CZ2', 'LINE', ], + 'CZ0', 'CZ1', 'CZ2', 'LINE', 'PA', ], - 'Misc' : ['CONFIG', 'DIR_OUT', 'FILTER', 'SKIPFILT', 'FIR_FILTER'] + 'Misc' : ['CONFIG', 'DIR_OUT', 'FILTER', 'SKIPFILT', 'FIR_FILTER', 'DIR_FILT'] } @@ -600,7 +600,7 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: self.fzmc = int(inputs['F_ZMC']) except: self.fzmc = 0 - self.logger.warning('Cannot find ZMC. Set to %d.'%(self.fzmc)) + self.logger.warning('Cannot find F_ZMC. Set to %d.'%(self.fzmc)) # Metallicity self.has_ZFIX = False @@ -667,6 +667,12 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: self.Zmax, self.Zmin = float(inputs['ZMAX']), float(inputs['ZMIN']) con_z = np.where((Zbpass >= self.Zmin) & (Zbpass <= self.Zmax)) self.Zall = Zbpass[con_z] + + if len(self.Zall) == 0: + self.logger.warning('No metallicity can be found. Available Zs are:') + self.logger.info(Zbpass) + sys.exit() + self.delZ = 0.0001 self.Zmax,self.Zmin = np.max(self.Zall), np.min(self.Zall) self.logger.info('Final list for log(Z_BPASS/Zsun) is : ') @@ -785,15 +791,16 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set: elif self.dust_model == 4: self.dust_model_name = 'KriekConroy' else: - self.logger.warning('Unknown index number for dust attenuation. Setting to Calzetti.') + self.logger.warning('Unknown index number for dust model. Setting to Calzetti.') self.dust_model = 0 self.dust_model_name = 'Calz' except: + self.logger.warning('Index number for dust model (DUST_MODEL) cannot be found. Setting to Calzetti.') self.dust_model = 0 self.dust_model_name = 'Calz' if self.verbose: - self.logger.info('Dust attenuation is set to %s\n'%self.dust_model_name) + self.logger.info('Dust model is set to %s\n'%self.dust_model_name) # If FIR data; try: @@ -1145,7 +1152,7 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, #zliml=0.01, zlim x_cz = self.dict['x'][con_cz] # Observed range NR_cz = self.dict['NR'][con_cz] if len(NR_cz) == 0: - self.logger.error('No data point exists at SN>%.1f'%(snlim)) + self.logger.error('No data point exists at SNR > `snlim` (%.1f)'%(snlim)) self.logger.error('Decrease `snlim` or turn `f_exclude_negative` to False') return False @@ -2408,6 +2415,7 @@ def search_redshift(self, dict, xm_tmp, fm_tmp, zliml=0.01, zlimu=6.0, delzz=0.0 fit_par_cz.add('C%d'%nn, value=1., min=0., max=1e5) def residual_z(pars,z): + '''''' vals = pars.valuesdict() xm_s = xm_tmp * (1+z) fm_s = np.zeros(len(xm_tmp),float) @@ -2416,7 +2424,6 @@ def residual_z(pars,z): fm_s += fm_tmp[nn,:] * pars['C%d'%nn] fint = interpolate.interp1d(xm_s, fm_s, kind='nearest', fill_value="extrapolate") - #fm_int = np.interp(xobs, xm_s, fm_s) fm_int = fint(xobs) if fcon is None: @@ -2427,7 +2434,6 @@ def residual_z(pars,z): # Start redshift search; for zz in range(len(zspace)): - # Best fit out_cz = minimize(residual_z, fit_par_cz, args=([zspace[zz]]), method=method) keys = fit_report(out_cz).split('\n') diff --git a/gsf/function.py b/gsf/function.py index c9eed43..67ec8fc 100644 --- a/gsf/function.py +++ b/gsf/function.py @@ -739,7 +739,7 @@ def apply_dust(yy, xx, nr, Av, dust_model=0): elif dust_model == 2: # LMC yyd, xxd, nrd = dust_gen(xx, yy, Av, nr, Rv=4.05, gamma=-0.06, Eb=2.8) elif dust_model == 3: # SMC - yyd, xxd, nrd = dust_smc(xx, yy, Av, nr, Rv=2.74, x0=4.703, gamma=1.212, lmlimu=3.115, f_Alam=False) + yyd, xxd, nrd = dust_smc(xx, yy, Av, nr, Rv=2.74, x0=4.703, gamma=1.212, f_Alam=False) elif dust_model == 4: # Kriek&Conroy with gamma=-0.2 yyd, xxd, nrd = dust_kc(xx, yy, Av, nr, Rv=4.05, gamma=-0.2) else: @@ -747,16 +747,14 @@ def apply_dust(yy, xx, nr, Av, dust_model=0): return yyd, xxd, nrd -def dust_smc(lm, fl, Av, nr, Rv=2.74, x0=4.703, gamma=1.212, lmlimu=3.115, f_Alam=False): +def dust_smc(lm, fl, Av, nr, Rv=2.74, x0=4.6, gamma=1.0, f_Alam=False): ''' - For general purpose (Noll+09). - This function is much better than previous, but is hard to impliment for the current version. - A difference from dust_gen is Eb is defined as a function of gamma. + Dust law for SMC (Gordon+03). Parameters ---------- lm : float array - wavelength, at RF. + wavelength in AA, at RF. fl : float array fnu Av : float @@ -772,42 +770,42 @@ def dust_smc(lm, fl, Av, nr, Rv=2.74, x0=4.703, gamma=1.212, lmlimu=3.115, f_Ala # if any(np.diff(lm)<0): # print('Something is wrong in lm: dust_smc of function.py') - lmm = lm/10000. # in micron - con1 = (lmm<=0.63) - con2 = (lmm>0.63) & (lmm<=lmlimu) - con3 = (lmm>lmlimu) - - #nr0 = nr[con0] - nr1 = nr[con1] - nr2 = nr[con2] - nr3 = nr[con3] - - #lmm0 = lmm[con0] - lmm1 = lmm[con1] - lmm2 = lmm[con2] - lmm3 = lmm[con3] - - #fl0 = fl[con0] - fl1 = fl[con1] - fl2 = fl[con2] - fl3 = fl[con3] - - nrd = np.concatenate([nr1,nr2,nr3]) - lmmc = np.concatenate([lmm1,lmm2,lmm3]) - flc = np.concatenate([fl1,fl2,fl3]) + lmm = lm/10000. # into micron + nrd = nr #np.concatenate([nr1,nr2,nr3]) + lmmc = lmm #np.concatenate([lmm1,lmm2,lmm3]) + flc = fl #np.concatenate([fl1,fl2,fl3]) # Using average for SMC; - c1,c2,c3,c4 = -0.856, 1.038, 3.215, 0.107 + c1,c2,c3,c4 = -4.959, 2.264, 0.389, 0.461 + # SMC Wing Sample; + # c1,c2,c3,c4 = -0.856, 1.038, 3.215, 0.107 + x = 1./lmmc - Dx = x**2 / ((x**2-x0**2)**2 +x**2*gamma**2) - Fx = 0.5329 * (x - 5.9)**2 + 0.05644 * (x-5.9)**3 + Dx = x**2 / ((x**2-x0**2)**2 + x**2*gamma**2) + Fx = 0.5392 * (x - 5.9)**2 + 0.05644 * (x-5.9)**3 con_fx = (x<5.9) Fx[con_fx] = 0 - EBlam_to_EB = c1+c2*x+c3*Dx+c4*Fx + EBlam_to_EB = c1 + c2*x + c3*Dx + c4*Fx Alam = Av / Rv * EBlam_to_EB fl_cor = flc[:] * 10**(-0.4*Alam[:]) + if False: + import matplotlib.pyplot as plt + plt.close() + xs = np.arange(0.5, 10, 0.01) + x = xs + Dx = x**2 / ((x**2-x0**2)**2 + x**2*gamma**2) + Fx = 0.5392 * (x - 5.9)**2 + 0.05644 * (x-5.9)**3 + con_fx = (x<5.9) + Fx[con_fx] = 0 + EBlam_to_EB = c1 + c2*x + c3*Dx + c4*Fx + Av = 2.0 + Alam = Av / Rv * EBlam_to_EB + plt.scatter(x, EBlam_to_EB, ) + plt.show() + hoge + if f_Alam: return fl_cor, lmmc*10000., nrd, Alam else: diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py index 04cabe7..a58ead7 100644 --- a/gsf/maketmp_filt.py +++ b/gsf/maketmp_filt.py @@ -1762,7 +1762,10 @@ def show_template_parameters(af): print('nimf :',af['nimf']) print('age/Gyr :',af['age'][:]) print('logZ/Zsun:',af['Z'][:]) - print('logU :',np.arange(af['logUMIN'], af['logUMAX']+0.01, af['DELlogU'])) + try: + print('logU :',np.arange(af['logUMIN'], af['logUMAX']+0.01, af['DELlogU'])) + except: + pass print('\n') return diff --git a/gsf/maketmp_z0.py b/gsf/maketmp_z0.py index 777daa5..905e616 100644 --- a/gsf/maketmp_z0.py +++ b/gsf/maketmp_z0.py @@ -553,6 +553,17 @@ def make_tmp_z0_bpass(MB, lammin=100, lammax=160000, Zforce=None, Zsun=0.02): 'version_gsf': gsf.__version__ } + tree.update({'age': MB.age}) + tree.update({'Z': MB.Zall}) + if fneb: + tree.update({'logUMIN': MB.logUMIN}) + tree.update({'logUMAX': MB.logUMAX}) + tree.update({'DELlogU': MB.DELlogU}) + if MB.fagn: + tree.update({'AGNTAUMIN': MB.AGNTAUMIN}) + tree.update({'AGNTAUMAX': MB.AGNTAUMAX}) + tree.update({'DELAGNTAU': MB.DELAGNTAU}) + # ASDF tree_spec.update({'wavelength': wave}) flagz = True @@ -579,6 +590,7 @@ def make_tmp_z0_bpass(MB, lammin=100, lammax=160000, Zforce=None, Zsun=0.02): col4 = fits.Column(name='tau_'+str(zz), format='E', unit='Gyr', array=tau_age) tree_ML.update({'realtau_'+str(zz): ms}) + # Write; for aa in range(len(age)): tree.update({'realtau%d(Gyr)'%(aa): tau_age[aa]}) diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py index 0598122..9b561fd 100644 --- a/gsf/plot_sed.py +++ b/gsf/plot_sed.py @@ -1077,6 +1077,8 @@ def plot_sed(MB, flim=0.01, fil_path='./', scale=None, f_chind=True, figpdf=Fals hdr['library'] = LIBRARY hdr['nimf'] = nimf hdr['scale'] = scale + hdr['dust model'] = MB.dust_model_name + hdr['ndust model'] = MB.dust_model try: # Chi square: @@ -1184,7 +1186,6 @@ def plot_sed(MB, flim=0.01, fil_path='./', scale=None, f_chind=True, figpdf=Fals tree_spec['header'].update({'%s'%key: hdr[key] * u.solMass / u.yr}) else: tree_spec['header'].update({'%s'%key: hdr[key]}) - # BB; Cnu_to_Jy = 10**((23.9-m0set)) # fnu_mzpset to microJy. So the final output SED library has uJy. # tree_spec['model'].update({'wave_bb': lbb * u.AA}) @@ -2761,12 +2762,17 @@ def plot_filter(MB, ax, ymax, scl=0.3, cmap='gist_rainbow', alp=0.4, def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax:int=1000, TMIN=0.0001, tau_lim=0.01, f_plot_filter=True, - scale=1e-19, NRbb_lim=10000, save_pcl=True, return_figure=False, SNlim=1, tset_SFR_SED=0.1): + scale=1e-19, NRbb_lim=10000, save_pcl=True, return_figure=False, SNlim=1, tset_SFR_SED=0.1, use_SFR_UV=True): ''' Purpose ------- For summary. In the same format as plot_corner_physparam_frame. + Parameters + ---------- + use_SFR_UV : bool + if True, SFR_UV will be used, instead of SFR_SFH. + Notes ----- Tau model not supported. @@ -2784,10 +2790,6 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax:in age = MB.age c = MB.c - tau0 = MB.tau0 - dust_model = MB.dust_model - DIR_TMP = MB.DIR_TMP - Txmax = np.max(age) + 1.0 ########################### @@ -2828,21 +2830,12 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax:in Z16 = np.zeros(len(age), dtype='float') Z84 = np.zeros(len(age), dtype='float') - NZbest = np.zeros(len(age), dtype='int') for aa in range(len(age)): Z50[aa] = hdul[1].data['Z'+str(aa)][1] Z16[aa] = hdul[1].data['Z'+str(aa)][0] Z84[aa] = hdul[1].data['Z'+str(aa)][2] - ZZ50 = np.sum(Z50*A50)/np.sum(A50) # Light weighted Z. - chi = hdul[1].data['chi'][0] - chin = hdul[1].data['chi'][1] - fitc = chin - Cz0 = hdul[0].header['Cz0'] - Cz1 = hdul[0].header['Cz1'] - Cz2 = hdul[0].header['Cz2'] zbes = hdul[0].header['z'] - zscl = (1.+zbes) # plot Configuration if MB.fzmc == 1: @@ -3105,17 +3098,6 @@ def plot_corner_physparam_summary(MB, fig=None, out_ind=0, DIR_OUT='./', mmax:in if SFmax > 0.5e4: SFmax = 0.5e4 - if MB.fzmc == 1: - if use_pickl: - NPARmin = [np.log10(M16)-.1, np.log10(SFminmc)-.1, np.log10(Tsmin/AMtmp16)-0.1, Av16-0.1, np.log10(Zsmin/AMtmp16)-0.2, np.nanpercentile(samples['zmc'],1)-0.1] - NPARmax = [np.log10(M84)+.1, np.log10(SFmaxmc)+.1, np.log10(Tsmax/AMtmp84)+0.2, Av84+0.1, np.log10(Zsmax/AMtmp84)+0.2, np.nanpercentile(samples['zmc'],99)+0.1] - else: - NPARmin = [np.log10(M16)-.1, np.log10(SFminmc)-.1, np.log10(Tsmin/AMtmp16)-0.1, Av16-0.1, np.log10(Zsmin/AMtmp16)-0.2, np.nanpercentile(list(samples['zmc'].values()),1)-0.1] - NPARmax = [np.log10(M84)+.1, np.log10(SFmaxmc)+.1, np.log10(Tsmax/AMtmp84)+0.2, Av84+0.1, np.log10(Zsmax/AMtmp84)+0.2, np.nanpercentile(list(samples['zmc'].values()),99)+0.1] - else: - NPARmin = [np.log10(M16)-.1, np.log10(SFminmc)-.1, np.log10(Tsmin/AMtmp16)-0.1, Av16-0.1, np.log10(Zsmin/AMtmp16)-0.2] - NPARmax = [np.log10(M84)+.1, np.log10(SFmaxmc)+.1, np.log10(Tsmax/AMtmp84)+0.2, Av84+0.1, np.log10(Zsmax/AMtmp84)+0.2] - # For redshift if zbes<2: zred = [zbes, 2, 3, 6] @@ -3167,6 +3149,7 @@ def density_estimation(m1, m2): MUV = np.zeros(mmax, dtype=float) for kk in range(0,mmax,1): + delt_tot = 0 nr = np.random.randint(nshape_sample) try: @@ -3283,11 +3266,6 @@ def density_estimation(m1, m2): Ztmp[kk] = np.log10(Ztmp[kk]) Ttmp[kk] = np.log10(Ttmp[kk]) - if MB.fzmc == 1: - NPAR = [lmtmp[:kk+1], SFR_SED[:kk+1], Ttmp[:kk+1], Avtmp[:kk+1], Ztmp[:kk+1], redshifttmp[:kk+1]] - else: - NPAR = [lmtmp[:kk+1], SFR_SED[:kk+1], Ttmp[:kk+1], Avtmp[:kk+1], Ztmp[:kk+1]] - # store other params; # Get FUV flux density at 10pc; if MB.fzmc == 1: @@ -3318,112 +3296,139 @@ def density_estimation(m1, m2): UVJ[kk,2] = -2.5*np.log10(fconv[2]/fconv[3]) UVJ[kk,3] = -2.5*np.log10(fconv[4]/fconv[3]) - # This should happen at the last kk; - if kk == mmax-1: - # Histogram - for i, x in enumerate(Par): - ax = axes[i, i] - x1min, x1max = NPARmin[i], NPARmax[i] - nbin = 50 - binwidth1 = (x1max-x1min)/nbin - bins1 = np.arange(x1min, x1max + binwidth1, binwidth1) - n, bins, patches = ax.hist(NPAR[i], bins=bins1, orientation='vertical', color='b', histtype='stepfilled', alpha=0.6) - yy = np.arange(0,np.max(n)*1.3,1) + # Summary; + kk = mmax - 1 + if MB.fzmc == 1: + NPAR = [lmtmp[:kk+1], SFR_SED[:kk+1], Ttmp[:kk+1], Avtmp[:kk+1], Ztmp[:kk+1], redshifttmp[:kk+1]] + else: + NPAR = [lmtmp[:kk+1], SFR_SED[:kk+1], Ttmp[:kk+1], Avtmp[:kk+1], Ztmp[:kk+1]] - try: - ax.plot(yy*0+np.percentile(NPAR[i],16), yy, linestyle='--', color='gray', lw=1) - ax.plot(yy*0+np.percentile(NPAR[i],84), yy, linestyle='--', color='gray', lw=1) - ax.plot(yy*0+np.percentile(NPAR[i],50), yy, linestyle='-', color='gray', lw=1) - except: - MB.logger.warning('Failed at i,x=%d,%d'%(i,x)) + # Define range; + if MB.fzmc == 1: + if use_pickl: + NPARmin = [np.log10(M16)-.1, np.log10(SFminmc)-.1, np.log10(Tsmin/AMtmp16)-0.1, Av16-0.1, np.log10(Zsmin/AMtmp16)-0.2, np.nanpercentile(samples['zmc'],1)-0.1] + NPARmax = [np.log10(M84)+.1, np.log10(SFmaxmc)+.1, np.log10(Tsmax/AMtmp84)+0.2, Av84+0.1, np.log10(Zsmax/AMtmp84)+0.2, np.nanpercentile(samples['zmc'],99)+0.1] + else: + NPARmin = [np.log10(M16)-.1, np.log10(SFminmc)-.1, np.log10(Tsmin/AMtmp16)-0.1, Av16-0.1, np.log10(Zsmin/AMtmp16)-0.2, np.nanpercentile(list(samples['zmc'].values()),1)-0.1] + NPARmax = [np.log10(M84)+.1, np.log10(SFmaxmc)+.1, np.log10(Tsmax/AMtmp84)+0.2, Av84+0.1, np.log10(Zsmax/AMtmp84)+0.2, np.nanpercentile(list(samples['zmc'].values()),99)+0.1] + else: + NPARmin = [np.log10(M16)-.1, np.log10(SFminmc)-.1, np.log10(Tsmin/AMtmp16)-0.1, Av16-0.1, np.log10(Zsmin/AMtmp16)-0.2] + NPARmax = [np.log10(M84)+.1, np.log10(SFmaxmc)+.1, np.log10(Tsmax/AMtmp84)+0.2, Av84+0.1, np.log10(Zsmax/AMtmp84)+0.2] + + if use_SFR_UV: + NPAR[1] = np.log10(SFRUV) + NPARmin[1] = np.nanmin(NPAR[1]) - 0.1 + NPARmax[1] = np.nanmax(NPAR[1]) + 0.1 + Par[1] = '$\log \mathrm{SFR_{UV}}/M_\odot \mathrm{yr}^{-1}$' + + # This should happen at the last kk; + if kk == mmax-1: + # Histogram + for i, x in enumerate(Par): + ax = axes[i, i] + x1min, x1max = NPARmin[i], NPARmax[i] + nbin = 50 + binwidth1 = (x1max-x1min)/nbin + bins1 = np.arange(x1min, x1max + binwidth1, binwidth1) + n, bins, patches = ax.hist(NPAR[i], bins=bins1, orientation='vertical', color='b', histtype='stepfilled', alpha=0.6) + yy = np.arange(0,np.max(n)*1.3,1) + + try: + ax.plot(yy*0+np.percentile(NPAR[i],16), yy, linestyle='--', color='gray', lw=1) + ax.plot(yy*0+np.percentile(NPAR[i],50), yy, linestyle='-', color='gray', lw=1) + ax.plot(yy*0+np.percentile(NPAR[i],84), yy, linestyle='--', color='gray', lw=1) + except: + MB.logger.warning('Failed at i,x=%d,%d'%(i,x)) + + ax.set_xlim(x1min, x1max) + ax.set_yticklabels([]) + if i == K-1: + ax.set_xlabel('%s'%(Par[i]), fontsize=12) + if i < K-1: + ax.set_xticklabels([]) + + # save pck; + if save_pcl: + if MB.fzmc == 1: + NPAR_LIB = {'logM_stel':lmtmp[:kk+1], 'logSFR':SFR_SED[:kk+1], 'logT_MW':Ttmp[:kk+1], 'AV':Avtmp[:kk+1], 'logZ_MW':Ztmp[:kk+1], 'z':redshifttmp[:kk+1], + 'MUV':MUV[:kk+1], 'Luv1600':Luv16[:kk+1], 'beta_UV':betas[:kk+1], 'SFRUV':SFRUV[:kk+1], 'SFRUV_UNCOR':SFRUV_UNCOR[:kk+1] + } + else: + NPAR_LIB = {'logM_stel':lmtmp[:kk+1], 'logSFR':SFR_SED[:kk+1], 'logT_MW':Ttmp[:kk+1], 'AV':Avtmp[:kk+1], 'logZ_MW':Ztmp[:kk+1], + 'MUV':MUV[:kk+1], 'Luv1600':Luv16[:kk+1], 'beta_UV':betas[:kk+1], 'SFRUV':SFRUV[:kk+1], 'SFRUV_UNCOR':SFRUV_UNCOR[:kk+1] + } + + # UVJ; + for cc in range(len(UVJ[0,:])): + NPAR_LIB['COR_RF_%d'%cc] = UVJ[:kk+1,cc] + + use_pickl = False + if use_pickl: + cpklname = os.path.join(MB.DIR_OUT, 'chain_' + MB.ID + '_phys.cpkl') + savecpkl({'chain':NPAR_LIB, + 'burnin':burnin, 'nwalkers':nwalk,'niter':nmc,'ndim':ndim}, + cpklname) # Already burn in + else: + cpklname = os.path.join(MB.DIR_OUT, 'chain_' + MB.ID + '_phys.asdf') + tree = {'chain':NPAR_LIB, 'burnin':burnin, 'nwalkers':nwalk,'niter':nmc,'ndim':ndim} + af = asdf.AsdfFile(tree) + af.write_to(cpklname, all_array_compression='zlib') + + + # Scatter and contour plot; + alp_sct = 0.1 / (mmax/1000) + for i, x in enumerate(Par): + for j, _ in enumerate(Par): + if i > j: + ax = axes[i, j] + ax.scatter(NPAR[j], NPAR[i], c='b', s=1, marker='.', alpha=alp_sct) + ax.set_xlabel('%s'%(Par[j]), fontsize=12) + + if kk == mmax-1: + try: + Xcont, Ycont, Zcont = density_estimation(NPAR[j], NPAR[i]) + mZ = np.max(Zcont) + ax.contour(Xcont, Ycont, Zcont, levels=[0.68*mZ, 0.95*mZ, 0.99*mZ], linewidths=[0.8,0.5,0.3], colors='orange') + except: + print('Error occurs when density estimation. Maybe because some params are fixed.') + pass + + x1min, x1max = NPARmin[j], NPARmax[j] + y1min, y1max = NPARmin[i], NPARmax[i] ax.set_xlim(x1min, x1max) + ax.set_ylim(y1min, y1max) + + if j==0: + ax.set_ylabel('%s'%(Par[i]), fontsize=12) + if j>0: + ax.set_yticklabels([]) + if i j: - ax = axes[i, j] - ax.scatter(NPAR[j], NPAR[i], c='b', s=1, marker='.', alpha=0.01) - ax.set_xlabel('%s'%(Par[j]), fontsize=12) - - if kk == mmax-1: - try: - Xcont, Ycont, Zcont = density_estimation(NPAR[j], NPAR[i]) - mZ = np.max(Zcont) - ax.contour(Xcont, Ycont, Zcont, levels=[0.68*mZ, 0.95*mZ, 0.99*mZ], linewidths=[0.8,0.5,0.3], colors='orange') - except: - print('Error occurs when density estimation. Maybe because some params are fixed.') - pass - - x1min, x1max = NPARmin[j], NPARmax[j] - y1min, y1max = NPARmin[i], NPARmax[i] - ax.set_xlim(x1min, x1max) - ax.set_ylim(y1min, y1max) - - if j==0: - ax.set_ylabel('%s'%(Par[i]), fontsize=12) - if j>0: - ax.set_yticklabels([]) - if i