Skip to content

Commit

Permalink
Merge pull request #45 from mtakahiro/nirspec
Browse files Browse the repository at this point in the history 
Nirspec
  • Loading branch information
@mtakahiro
mtakahiro authored Feb 1, 2023
2 parents 0b5fbaf + 27d5d4d commit e87b7cd
Show file tree
Hide file tree
Showing 9 changed files with 701 additions and 921 deletions.
451 changes: 261 additions & 190 deletions gsf/fitting.py
View file

Large diffs are not rendered by default.

29 changes: 28 additions & 1 deletion gsf/function.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@
import scipy.interpolate as interpolate
from scipy.interpolate import interp1d
import logging
from colorama import Fore, Back, Style
from datetime import datetime

################
# Line library
Expand All @@ -17,6 +19,32 @@
c = 3.e18 # A/s


def print_err(msg, exit=False, details=None):
'''
'''
now = datetime.now()
print(Fore.RED)
print('$$$ =================== $$$')
print('$$$ gsf error message $$$')
print('%s'%now)
if not details == None:
# @@@ This does not work??
exc_type, exc_obj, exc_tb = details
fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
print(exc_type, fname, exc_tb.tb_lineno)

print('$$$ =================== $$$')
print(Fore.RED)
print(msg)
print(Fore.RED)
print(Style.RESET_ALL)

if exit:
print(Fore.CYAN)
print('Exiting.')
print(Style.RESET_ALL)
sys.exit()

def str2bool(v):
'''
'''
Expand Down Expand Up @@ -285,7 +313,6 @@ def get_leastsq(MB, ZZtmp, fneld, age, fit_params, residual, fy, ey, wht, ID0,
fit_params['Z'+str(aa)].value = ZZ

f_fir = False

out_tmp = minimize(residual, fit_params, args=(fy, ey, wht, f_fir),
method=fit_name, kws={'f_only_spec':f_only_spec})

Expand Down
111 changes: 43 additions & 68 deletions gsf/gsf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,19 +4,22 @@
import os.path
import string
from astropy.cosmology import WMAP9 as cosmo
from astropy.io import fits
import asdf

# From gsf
from .fitting import Mainbody
from .maketmp_filt import maketemp,maketemp_tau
from .function_class import Func,Func_tau
from .basic_func import Basic,Basic_tau
from .function import read_input

import timeit
start = timeit.default_timer()


def run_gsf_template(inputs, fplt=0, tau_lim=0.001, idman=None, nthin=1, delwave=10):
def run_gsf_template(inputs, fplt=0, tau_lim=0.001, idman=None, nthin=1, delwave=10,
f_IGM=True):
'''
Purpose
-------
Expand Down Expand Up @@ -44,7 +47,7 @@ def run_gsf_template(inputs, fplt=0, tau_lim=0.001, idman=None, nthin=1, delwave
#
# 0. Make basic templates
#
if fplt == 0 or fplt == 1:
if fplt == 0 or fplt == 1 or fplt == 2:
#
# 0. Make basic templates
#
Expand Down Expand Up @@ -72,44 +75,21 @@ def run_gsf_template(inputs, fplt=0, tau_lim=0.001, idman=None, nthin=1, delwave
# 1. Start making redshifted templates.
#
if MB.SFH_FORM == -99:
maketemp(MB, tau_lim=tau_lim, nthin=nthin, delwave=delwave)
maketemp(MB, tau_lim=tau_lim, nthin=nthin, delwave=delwave, f_IGM=f_IGM)
else:
maketemp_tau(MB, tau_lim=tau_lim, nthin=nthin, delwave=delwave)

# Read temp from asdf;
# This has to happend after fplt==1 and before fplt>=2.
# MB.af = asdf.open(MB.DIR_TMP + 'spec_all_' + MB.ID + '.asdf')

'''
#
# 2. Load templates
#
MB.lib = MB.fnc.open_spec_fits(MB, fall=0)
MB.lib_all = MB.fnc.open_spec_fits(MB, fall=1)
if MB.f_dust:
MB.lib_dust = MB.fnc.open_spec_dust_fits(MB, fall=0)
MB.lib_dust_all = MB.fnc.open_spec_dust_fits(MB, fall=1)
'''

# How to get SED?
if False:
import matplotlib.pyplot as plt
for T in MB.age[:1]:
y0, x0 = MB.fnc.get_template(MB.lib_all, Amp=1.0, T=T, Av=0.0, Z=-1.0, zgal=MB.zgal)
plt.plot(x0/(1.+MB.zgal),y0,linestyle='-',lw=1.0, label='$T=%.2f$\n$z=%.2f$'%(T,MB.zgal))
plt.xlim(600,20000)
plt.xlabel('Rest frame wavelength')
plt.xscale('log')
plt.legend(loc=0)
plt.show()
maketemp_tau(MB, tau_lim=tau_lim, nthin=nthin, delwave=delwave, f_IGM=f_IGM)

return MB


def run_gsf_all(parfile, fplt, cornerplot=True, f_Alog=True, idman=None, zman=None, f_label=True, f_symbol=True,
f_SFMS=False, f_fill=True, save_sed=True, figpdf=False, mmax=300, skip_sfh=False, f_fancyplot=False,
skip_zhist=False, tau_lim=0.001, tset_SFR_SED=0.1, f_shuffle=False, amp_shuffle=1e-2, Zini=None,
nthin=1, delwave=1, f_plot_resid=False, scale=1e-19, f_plot_filter=True, f_prior_sfh=False, norder_sfh_prior=3):
def run_gsf_all(parfile, fplt, cornerplot=True, f_plot_chain=True, f_Alog=True, idman:str=None,
zman=None, zman_min=None, zman_max=None, f_label=True, f_symbol=True,
f_SFMS=False, f_fill=True, save_sed=True, figpdf=False, mmax=300,
f_prior_sfh=False, norder_sfh_prior=3,
f_shuffle=False, amp_shuffle=1e-2, Zini=None, tau_lim=0.001,
skip_sfh=False, f_fancyplot=False, skip_zhist=False, f_sfh_yaxis_force=True, tset_SFR_SED=0.1,
nthin=1, delwave=1, f_plot_resid=False, scale=1e-19, f_plot_filter=True
):
'''
Purpose
-------
Expand All @@ -125,10 +105,9 @@ def run_gsf_all(parfile, fplt, cornerplot=True, f_Alog=True, idman=None, zman=No
######################
# Read from Input file
######################
from .function import read_input
inputs = read_input(parfile)

MB = Mainbody(inputs, c=3e18, Mpc_cm=3.08568025e+24, m0set=25.0, pixelscale=0.06, cosmo=cosmo, idman=idman, zman=zman)
MB = Mainbody(inputs, c=3e18, Mpc_cm=3.08568025e+24, m0set=25.0, pixelscale=0.06,
cosmo=cosmo, idman=idman, zman=zman, zman_min=zman_min, zman_max=zman_max)

# Register some params;
MB.tau_lim = tau_lim
Expand Down Expand Up @@ -175,17 +154,9 @@ def run_gsf_all(parfile, fplt, cornerplot=True, f_Alog=True, idman=None, zman=No
if not flag_suc:
sys.exit()

# Read temp from asdf;
# Template must be registered before fplt>=2.
# try:
# aftmp = MB.af
# except:
# MB.af = asdf.open(os.path.join(MB.DIR_TMP, 'spec_all_' + MB.ID + '.asdf'))

if fplt <= 2:
MB.zprev = MB.zgal
MB.ndim_keep = MB.ndim

#
# 1. Start making redshifted templates, at z=MB.zgal.
#
Expand All @@ -194,10 +165,13 @@ def run_gsf_all(parfile, fplt, cornerplot=True, f_Alog=True, idman=None, zman=No
else:
flag_suc = maketemp_tau(MB, tau_lim=MB.tau_lim, nthin=MB.nthin, delwave=MB.delwave)

if not flag_suc:
return False

#
# 2. Main fitting part.
#
flag_suc = MB.main(cornerplot=cornerplot, f_shuffle=f_shuffle, amp_shuffle=amp_shuffle, Zini=Zini,
flag_suc = MB.main(cornerplot=cornerplot, f_plot_chain=f_plot_chain, f_shuffle=f_shuffle, amp_shuffle=amp_shuffle, Zini=Zini,
f_prior_sfh=f_prior_sfh, norder_sfh_prior=norder_sfh_prior)

while (flag_suc and flag_suc!=2):
Expand All @@ -217,7 +191,7 @@ def run_gsf_all(parfile, fplt, cornerplot=True, f_Alog=True, idman=None, zman=No
print('Going into another round with updated templates and redshift.')
print('\n\n')

flag_suc = MB.main(cornerplot=cornerplot, f_shuffle=f_shuffle, amp_shuffle=amp_shuffle, Zini=Zini,
flag_suc = MB.main(cornerplot=cornerplot, f_plot_chain=f_plot_chain, f_shuffle=f_shuffle, amp_shuffle=amp_shuffle, Zini=Zini,
f_prior_sfh=f_prior_sfh, norder_sfh_prior=norder_sfh_prior)

# Total calculation time
Expand All @@ -229,18 +203,17 @@ def run_gsf_all(parfile, fplt, cornerplot=True, f_Alog=True, idman=None, zman=No

if fplt <= 3 and flag_suc:

if True: #fplt >= 3: # If main fitting has not been done;
# make redshifted templates and register data;

# Use the final redshift;
from astropy.io import fits
hd_sum = fits.open(os.path.join(MB.DIR_OUT, 'summary_%s.fits'%MB.ID))[0].header
MB.zgal = hd_sum['ZMC']

# Use the final redshift;
hd_sum = fits.open(os.path.join(MB.DIR_OUT, 'summary_%s.fits'%MB.ID))[0].header
MB.zgal = hd_sum['ZMC']

if not MB.ztemplate:
if MB.SFH_FORM == -99:
flag_suc = maketemp(MB, tau_lim=tau_lim, nthin=nthin, delwave=delwave)
else:
flag_suc = maketemp_tau(MB, tau_lim=tau_lim, nthin=nthin, delwave=delwave)
if not flag_suc:
return False

if MB.SFH_FORM == -99:
from .plot_sfh import plot_sfh
Expand All @@ -252,32 +225,34 @@ def run_gsf_all(parfile, fplt, cornerplot=True, f_Alog=True, idman=None, zman=No
if not skip_sfh:
plot_sfh(MB, fil_path=MB.DIR_FILT, mmax=mmax,
dust_model=MB.dust_model, DIR_TMP=MB.DIR_TMP, f_silence=True,
f_SFMS=f_SFMS, f_symbol=f_symbol, skip_zhist=skip_zhist, tau_lim=tau_lim, tset_SFR_SED=tset_SFR_SED)
f_SFMS=f_SFMS, f_symbol=f_symbol, skip_zhist=skip_zhist,
tau_lim=tau_lim, tset_SFR_SED=tset_SFR_SED, f_sfh_yaxis_force=f_sfh_yaxis_force)

plot_sed(MB, fil_path=MB.DIR_FILT,
figpdf=figpdf, save_sed=save_sed, mmax=mmax,
dust_model=MB.dust_model, DIR_TMP=MB.DIR_TMP, f_label=f_label, f_fill=f_fill,
f_fancyplot=f_fancyplot, f_plot_resid=f_plot_resid, scale=scale, f_plot_filter=f_plot_filter)

'''
if fplt == 4:
from .plot_sfh import get_evolv
get_evolv(MB, MB.ID, MB.PA, MB.Zall, MB.age, f_comp=MB.ftaucomp, fil_path=MB.DIR_FILT, inputs=MB.inputs, dust_model=MB.dust_model, DIR_TMP=MB.DIR_TMP)
if fplt == 5:
from .plot_sfh import plot_evolv
plot_evolv(MB, MB.ID, MB.PA, MB.Zall, MB.age, f_comp=MB.ftaucomp, fil_path=MB.DIR_FILT, inputs=MB.inputs, dust_model=MB.dust_model, DIR_TMP=MB.DIR_TMP, nmc=10)
'''

if fplt == 6:
# Use the final redshift;
hd_sum = fits.open(os.path.join(MB.DIR_OUT, 'summary_%s.fits'%MB.ID))[0].header
MB.zgal = hd_sum['ZMC']

if not MB.ztemplate:
if MB.SFH_FORM == -99:
flag_suc = maketemp(MB, tau_lim=tau_lim, nthin=nthin, delwave=delwave)
else:
flag_suc = maketemp_tau(MB, tau_lim=tau_lim, nthin=nthin, delwave=delwave)

if MB.SFH_FORM == -99:
from .plot_sed import plot_corner_physparam_frame,plot_corner_physparam_summary
plot_corner_physparam_summary(MB)
else:
#from .plot_sed_logA import plot_corner_physparam_summary_tau as plot_corner_physparam_summary
print('One for Tau model is TBD...')

return MB


if __name__ == "__main__":
'''
Expand Down
Loading

0 comments on commit e87b7cd

Please sign in to comment.