Skip to content

Commit

Permalink
minor.
Browse files Browse the repository at this point in the history
  • Loading branch information
Takahiro Morishita committed Jul 26, 2020
1 parent 11778f6 commit cf3e93f
Show file tree
Hide file tree
Showing 7 changed files with 50 additions and 40 deletions.
17 changes: 11 additions & 6 deletions gsf/fitting.py
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ def __init__(self, inputs, c=3e18, Mpc_cm=3.08568025e+24, m0set=25.0, pixelscale
def update_input(self, inputs, c=3e18, Mpc_cm=3.08568025e+24, m0set=25.0, pixelscale=0.06, Lsun=3.839*1e33, cosmo=None):
'''
INPUT:
==========
======
parfile: Ascii file that lists parameters for everything.
Mpc_cm : cm/Mpc
Expand Down Expand Up @@ -292,11 +292,13 @@ def get_lines(self, LW0):

def read_data(self, Cz0, Cz1, zgal, add_fir=False):
'''
Input:
======
Cz0, Cz1 : Normalization coeffs for grism spectra.
zgal : Current redshift estimate.
Note:
=======
=====
Can be used for any SFH
'''
Expand Down Expand Up @@ -408,12 +410,12 @@ def read_data(self, Cz0, Cz1, zgal, add_fir=False):
def search_redshift(self, dict, xm_tmp, fm_tmp, zliml=0.01, zlimu=6.0, delzz=0.01, lines=False, prior=None, method='powell'):
'''
Purpose:
=========
========
Search redshift space to find the best redshift and probability distribution.
Input:
=========
======
fm_tmp : a library for various templates. Should be in [ n * len(wavelength)].
xm_tmp : a wavelength array, common for the templates above, at z=0. Should be in [len(wavelength)].
Expand All @@ -425,7 +427,7 @@ def search_redshift(self, dict, xm_tmp, fm_tmp, zliml=0.01, zlimu=6.0, delzz=0.0
method : powell is more accurate. nelder is faster.
Return:
=========
=======
zspace :
chi2s :
Expand Down Expand Up @@ -1041,7 +1043,10 @@ def main(self, cornerplot=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0, f_move

################################
print('\nMinimizer Defined\n')
mini = Minimizer(class_post.lnprob, out.params, fcn_args=[dict['fy'],dict['ey'],dict['wht2'],self.f_dust], f_disp=self.f_disp, moves=emcee.moves.DEMove(sigma=1e-05, gamma0=None)) #, f_move=f_move)
mini = Minimizer(class_post.lnprob, out.params, fcn_args=[dict['fy'],dict['ey'],dict['wht2'],self.f_dust], f_disp=self.f_disp, \
moves=[(emcee.moves.DEMove(), 0.8), (emcee.moves.DESnookerMove(), 0.2),]
)
#moves=emcee.moves.DEMove(sigma=1e-05, gamma0=None))
print('######################')
print('### Starting emcee ###')
print('######################')
Expand Down
2 changes: 1 addition & 1 deletion gsf/function.py
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ def printProgressBar (iteration, total, prefix = '', suffix = '', decimals = 1,
if fill == None:
if float(percent) < 33:
fill = emojis[0]
elif float(percent) < 33:
elif float(percent) < 66:
fill = emojis[1]
elif float(percent) < 99:
fill = emojis[2]
Expand Down
12 changes: 7 additions & 5 deletions gsf/gsf.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
def run_gsf_template(inputs, fplt=0):
'''
Purpose:
==========
========
This is only for 0 and 1, to get templates.
Not for fitting, nor plotting.
Expand Down Expand Up @@ -130,7 +130,7 @@ def run_gsf_all(parfile, fplt, cornerplot=True):
#
MB.zprev = MB.zgal #zrecom # redshift from previous run

flag_suc = MB.main(0, cornerplot=cornerplot)
flag_suc = MB.main(cornerplot=cornerplot)

while (flag_suc and flag_suc!=-1):

Expand All @@ -142,7 +142,7 @@ def run_gsf_all(parfile, fplt, cornerplot=True):
print('Going into another trial with updated templates and redshift.')
print('\n\n')

flag_suc = MB.main(1, cornerplot=cornerplot)
flag_suc = MB.main(cornerplot=cornerplot)

# Total calculation time
stop = timeit.default_timer()
Expand All @@ -152,10 +152,12 @@ def run_gsf_all(parfile, fplt, cornerplot=True):
if fplt <= 3 and flag_suc != -1:
from .plot_sfh import plot_sfh
from .plot_sed import plot_sed
plot_sfh(MB, f_comp=MB.ftaucomp, fil_path=MB.DIR_FILT,

plot_sfh(MB, f_comp=MB.ftaucomp, fil_path=MB.DIR_FILT, mmax=100,
inputs=MB.inputs, dust_model=MB.dust_model, DIR_TMP=MB.DIR_TMP, f_SFMS=True)

plot_sed(MB, fil_path=MB.DIR_FILT,
figpdf=False, save_sed=True, inputs=MB.inputs, nmc_rand=1000,
figpdf=False, save_sed=True, inputs=MB.inputs, mmax=30,
dust_model=MB.dust_model, DIR_TMP=MB.DIR_TMP, f_label=True)


Expand Down
1 change: 1 addition & 0 deletions gsf/maketmp_filt.py
Original file line number Diff line number Diff line change
Expand Up @@ -585,6 +585,7 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=20000., ncolbb=10000):
for ii in range(len(ltmpbb[0,:])):
if SFILT[ii] in SKIPFILT:# data point to be skiped;
fw.write('%d %.5f %.5e %.5e\n'%(ii+ncolbb, ltmpbb[0,ii], 0.0, fbb[ii]))
#fw.write('%d %.5f %.5e %.5e\n'%(ii+ncolbb, ltmpbb[0,ii], 0.0, 1000))
elif ebb[ii]>ebblim:
fw.write('%d %.5f 0 1000\n'%(ii+ncolbb, ltmpbb[0,ii]))
else:
Expand Down
23 changes: 11 additions & 12 deletions gsf/plot_sed.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@
col = ['violet', 'indigo', 'b', 'lightblue', 'lightgreen', 'g', 'orange', 'coral', 'r', 'darkred']#, 'k']
#col = ['darkred', 'r', 'coral','orange','g','lightgreen', 'lightblue', 'b','indigo','violet','k']

def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=False, save_sed=True, inputs=False, nmc_rand=300, dust_model=0, DIR_TMP='./templates/', f_label=False, f_bbbox=False, verbose=False, f_silence=True, f_fill=False, f_fancyplot=False):
def plot_sed(MB, flim=0.01, fil_path='./', scale=1e-19, f_chind=True, figpdf=False, save_sed=True, inputs=False, mmax=300, dust_model=0, DIR_TMP='./templates/', f_label=False, f_bbbox=False, verbose=False, f_silence=True, f_fill=False, f_fancyplot=False):
'''
Input:
======
Expand Down Expand Up @@ -550,27 +550,26 @@ def gaus(x,a,x0,sigma):
samples = res

# Saved template;
ytmp = np.zeros((nmc_rand,len(ysum)), dtype='float64')
ytmp_each = np.zeros((nmc_rand,len(ysum),len(age)), dtype='float64')
ytmp = np.zeros((mmax,len(ysum)), dtype='float64')
ytmp_each = np.zeros((mmax,len(ysum),len(age)), dtype='float64')

ytmpmax = np.zeros(len(ysum), dtype='float64')
ytmpmin = np.zeros(len(ysum), dtype='float64')

# MUV;
DL = MB.cosmo.luminosity_distance(zbes).value * Mpc_cm # Luminositydistance in cm
DL10 = Mpc_cm/1e6 * 10 # 10pc in cm
Fuv = np.zeros(nmc_rand, dtype='float64') # For Muv
Fuv28 = np.zeros(nmc_rand, dtype='float64') # For Fuv(1500-2800)
Lir = np.zeros(nmc_rand, dtype='float64') # For L(8-1000um)
UVJ = np.zeros((nmc_rand,4), dtype='float64') # For UVJ color;
Fuv = np.zeros(mmax, dtype='float64') # For Muv
Fuv28 = np.zeros(mmax, dtype='float64') # For Fuv(1500-2800)
Lir = np.zeros(mmax, dtype='float64') # For L(8-1000um)
UVJ = np.zeros((mmax,4), dtype='float64') # For UVJ color;

Cmznu = 10**((48.6+m0set)/(-2.5)) # Conversion from m0_25 to fnu

# From random chain;
alp=0.02
for kk in range(0,nmc_rand,1):
for kk in range(0,mmax,1):
nr = np.random.randint(len(samples['A0']))
#nr = np.random.randint(int(len(samples['A0'])*0.99), len(samples['A0']))
try:
Av_tmp = samples['Av'][nr]
except:
Expand Down Expand Up @@ -655,7 +654,7 @@ def gaus(x,a,x0,sigma):
# Do stuff...
time.sleep(0.01)
# Update Progress Bar
printProgressBar(kk, nmc_rand, prefix = 'Progress:', suffix = 'Complete', length = 40)
printProgressBar(kk, mmax, prefix = 'Progress:', suffix = 'Complete', length = 40)


#
Expand Down Expand Up @@ -726,13 +725,13 @@ def func_tmp(xint,eobs,fmodel):
ndim_eff -= 1
'''

con_up = (ey>0) & (fy/ey<=SNlim)
con_up = (fy==0) & (ey>0) & (fy/ey<=SNlim)
chi_nd = 0.0
if f_chind:
# Chi2 for non detection;
for nn in range(len(ey[con_up])):
#result = integrate.quad(lambda xint: func_tmp(xint, ey[con_up][nn]/SNlim, ysump[con_up][nn]), -ey[con_up][nn]/SNlim, ey[con_up][nn]/SNlim, limit=100)
result = integrate.quad(lambda xint: func_tmp(xint, ey[con_up][nn], ysump[con_up][nn]), -ey[con_up][nn], ey[con_up][nn], limit=100)
result = integrate.quad(lambda xint: func_tmp(xint, ey[con_up][nn], ysump[con_up][nn]), -ey[con_up][nn]*100, ey[con_up][nn], limit=100)
chi_nd += np.log(result[0])

# Number of degree;
Expand Down
22 changes: 12 additions & 10 deletions gsf/plot_sfh.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@

lcb = '#4682b4' # line color, blue

def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, fil_path = './FILT/', inputs=None, dust_model=0, DIR_TMP='./templates/',f_SFMS=False, verbose=False):
def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmin=0.08, Txmax=4, lmmin=9.5, fil_path = './FILT/', inputs=None, dust_model=0, DIR_TMP='./templates/',f_SFMS=False, verbose=False):
'''
Purpose:
========
Expand All @@ -47,9 +47,11 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
age = MB.age
nage = MB.nage #np.arange(0,len(age),1)
tau0 = MB.tau0 #[0.1,0.2,0.3]

age = np.asarray(age)

if Txmin > np.min(age):
Txmin = np.min(age) * 0.8

################
# RF colors.
home = os.path.expanduser('~')
Expand Down Expand Up @@ -482,13 +484,13 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
if f_comp == 1:
lsfru = np.max([lsfru, np.log10(np.max(sfr_exp*C_exp))])

ax1.set_xlim(0.008, Txmax)
ax1.set_xlim(Txmin, Txmax)
ax1.set_ylim(lsfrl, lsfru)
ax1.set_xscale('log')

ax2.set_ylabel('$\log M_*/M_\odot$', fontsize=12)

ax2.set_xlim(0.008, Txmax)
ax2.set_xlim(Txmin, Txmax)
ax2.set_ylim(y2min, y2max)
ax2.set_xscale('log')

Expand Down Expand Up @@ -591,7 +593,7 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
#ax3.set_xlabel('$t$ (Gyr)', fontsize=12)
#ax3.set_ylabel('$\log Z_*/Z_\odot$', fontsize=12)
y3min, y3max = np.min(Z), np.max(Z)
#ax3.set_xlim(0.008, Txmax)
#ax3.set_xlim(Txmin, Txmax)
#ax3.set_ylim(y3min, y3max)
#ax3.set_xscale('log')
#ax3.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))
Expand Down Expand Up @@ -628,21 +630,21 @@ def plot_sfh(MB, f_comp=0, flim=0.01, lsfrl=-1, mmax=1000, Txmax=4, lmmin=9.5, f
Tzz[zz] = 0.01

#ax3t.set_xscale('log')
#ax3t.set_xlim(0.008, Txmax)
#ax3t.set_xlim(Txmin, Txmax)

ax1.set_xlabel('$t_\mathrm{lookback}$/Gyr', fontsize=12)
ax2.set_xlabel('$t_\mathrm{lookback}$/Gyr', fontsize=12)
ax4.set_xlabel('$t_\mathrm{lookback}$/Gyr', fontsize=12)
ax4.set_ylabel('$\log Z_*/Z_\odot$', fontsize=12)

ax1t.set_xscale('log')
ax1t.set_xlim(0.008, Txmax)
ax1t.set_xlim(Txmin, Txmax)
ax2t.set_xscale('log')
ax2t.set_xlim(0.008, Txmax)
ax2t.set_xlim(Txmin, Txmax)
ax4t.set_xscale('log')
ax4t.set_xlim(0.008, Txmax)
ax4t.set_xlim(Txmin, Txmax)

ax4.set_xlim(0.008, Txmax)
ax4.set_xlim(Txmin, Txmax)
ax4.set_ylim(y3min-0.05, y3max)
ax4.set_xscale('log')

Expand Down
13 changes: 7 additions & 6 deletions gsf/posterior_flexible.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,12 +18,12 @@ def __init__(self, mainbody):
def residual(self, pars, fy, ey, wht, f_fir=False, out=False):
'''
Input:
==========
======
out : model as second output. For lnprob func.
f_fir : Bool. If dust component is on or off.
Returns:
==========
========
residual of model and data.
'''
Expand Down Expand Up @@ -76,7 +76,7 @@ def lnprob(self, pars, fy, ey, wht, f_fir, f_chind=True, SNlim=1.0):
'''
Returns:
=========
========
log posterior
'''
Expand All @@ -89,16 +89,17 @@ def lnprob(self, pars, fy, ey, wht, f_fir, f_chind=True, SNlim=1.0):

resid, model = self.residual(pars, fy, ey, wht, f_fir, out=True)
sig = np.sqrt(1./wht+f**2*model**2)
chi_nd = 0.0

chi_nd = 0
con_up = (fy==0) & (ey>0)
if f_chind:
# This does not improve, but costs time;
con_up = (fy==0) & (fy/ey<=SNlim) & (ey>0)
# This may be a bit cost of time;
for nn in range(len(ey[con_up])):
result = integrate.quad(lambda xint: self.func_tmp(xint, ey[con_up][nn]/SNlim, model[con_up][nn]), -ey[con_up][nn], ey[con_up][nn], limit=100)
chi_nd += np.log(result[0])

con_res = (model>=0) & (wht>0) & (fy>0) # Instead of model>0, model>=0 is for Lyman limit where flux=0.
#print(len(fy[con_up]),len(fy[con_res]))
lnlike = -0.5 * (np.sum(resid[con_res]**2 + np.log(2 * 3.14 * sig[con_res]**2)) - 2 * chi_nd)

else:
Expand Down

0 comments on commit cf3e93f

Please sign in to comment.