Skip to content

Some updates to read_alf.py #16

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
chloe-mt-cheng wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Some updates to read_alf.py #16

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
309 changes: 260 additions & 49 deletions scripts/read_alf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,15 +12,24 @@
class Alf(object):
def __init__(self, outfiles, read_mcmc=True, info=None, index=False):
self.outfiles = outfiles
## ND: New file reading code, don't need this after rehaul. Delete.
#self.legend = info['label']
#self.imf_type = info['imf_type']
self.nsample = None
#self.nsample = None
self.spectra = None
self.model_info = {}

if read_mcmc:
self.mcmc = np.loadtxt('{0}.mcmc'.format(self.outfiles))
results = ascii.read('{0}.sum'.format(self.outfiles))

# To keep track of labels for backwards compatibility of
# previous versions of alf
self.filetype = len(results.colnames)


## ND: New file reading code, don't need this after rehaul. Delete.
"""
with open('{0}.sum'.format(self.outfiles)) as f:
for line in f:
if line[0] == '#':
Expand All @@ -30,18 +39,57 @@ def __init__(self, outfiles, read_mcmc=True, info=None, index=False):
self.nchain = float(line.split('=')[1].strip())
elif 'Nsample' in line:
self.nsample = float(line.split('=')[1].strip())
"""

self.labels = np.array([
'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a',
'C', 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr',
'Mn', 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff',
'IMF1', 'IMF2', 'logfy', 'sigma2', 'velz2', 'logm7g',
'hotteff', 'loghot', 'fy_logage', 'logemline_h',
'logemline_oii', 'logemline_oiii', 'logemline_sii',
'logemline_ni', 'logemline_nii', 'logtrans', 'jitter',
'logsky', 'IMF3', 'IMF4', 'h3', 'h4',
'ML_v','ML_i','ML_k','MW_v', 'MW_i','MW_k'
])
with open('{0}.sum'.format(self.outfiles)) as f:
for line in f:
if line[0] == '#':
attr = line.replace('#', '').split('=')
if len(attr) == 2:
try:
self.model_info[attr[0].strip()] = int(attr[1].strip())
except:
self.model_info[attr[0].strip()] = attr[1].strip()
else:
break

# Older versions of alf have different labels/different label positions
# Current version has 53 labels
# Older versions have 52 or 50 (i.e. van Dokkum+17 has 50)
if self.filetype == 53:
self.labels = np.array([
'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a',
'C', 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr',
'Mn', 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff',
'IMF1', 'IMF2', 'logfy', 'sigma2', 'velz2', 'logm7g',
'hotteff', 'loghot', 'fy_logage', 'logemline_h',
'logemline_oii', 'logemline_oiii', 'logemline_sii',
'logemline_ni', 'logemline_nii', 'logtrans', 'jitter',
'logsky', 'IMF3', 'IMF4', 'h3', 'h4',
'ML_v','ML_i','ML_k','MW_v', 'MW_i','MW_k'
])
elif self.filetype == 52:
self.labels = np.array([
'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a',
'C', 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr',
'Mn', 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff',
'IMF1', 'IMF2', 'logfy', 'sigma2', 'velz2', 'logm7g',
'hotteff', 'loghot', 'fy_logage', 'logtrans', 'logemline_h',
'logemline_oiii', 'logemline_sii', 'logemline_ni',
'logemline_nii', 'jitter', 'IMF3', 'logsky', 'IMF4',
'h3', 'h4','ML_v','ML_i','ML_k','MW_v', 'MW_i','MW_k'
])
elif self.filetype == 50:
self.labels = np.array([
'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a', 'C',
'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr', 'Mn',
'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff', 'IMF1', 'IMF2',
'logfy', 'sigma2', 'velz2', 'logm7g', 'hotteff', 'loghot',
'fy_logage', 'logtrans', 'logemline_h', 'logemline_oiii',
'logemline_sii', 'logemline_ni', 'logemline_nii', 'jitter',
'IMF3', 'logsky', 'IMF4', 'ML_v', 'ML_i', 'ML_k', 'MW_v',
'MW_i', 'MW_k'
])

results = Table(results, names=self.labels)

Expand Down Expand Up @@ -72,48 +120,67 @@ def __init__(self, outfiles, read_mcmc=True, info=None, index=False):
# is filled in abundance_correct()
self.xFe = {}

self.results = results['Type',
'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a',
'C', 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr',
'Mn', 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff',
'IMF1', 'IMF2', 'logfy', 'sigma2', 'velz2', 'logm7g',
'hotteff', 'loghot', 'fy_logage', 'logemline_h',
'logemline_oii', 'logemline_oiii', 'logemline_sii',
'logemline_ni', 'logemline_nii', 'logtrans', 'jitter',
'logsky', 'IMF3', 'IMF4', 'h3', 'h4',
'ML_v','ML_i','ML_k','MW_v', 'MW_i','MW_k' ]
# Same filetype differentiations as above
if self.filetype == 53:
self.results = results['Type',
'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a',
'C', 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr',
'Mn', 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff',
'IMF1', 'IMF2', 'logfy', 'sigma2', 'velz2', 'logm7g',
'hotteff', 'loghot', 'fy_logage', 'logemline_h',
'logemline_oii', 'logemline_oiii', 'logemline_sii',
'logemline_ni', 'logemline_nii', 'logtrans', 'jitter',
'logsky', 'IMF3', 'IMF4', 'h3', 'h4',
'ML_v','ML_i','ML_k','MW_v', 'MW_i','MW_k' ]
elif self.filetype == 52:
self.results = results['Type',
'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a',
'C', 'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr',
'Mn', 'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff',
'IMF1', 'IMF2', 'logfy', 'sigma2', 'velz2', 'logm7g',
'hotteff', 'loghot', 'fy_logage', 'logtrans', 'logemline_h',
'logemline_oiii', 'logemline_sii', 'logemline_ni',
'logemline_nii', 'jitter', 'IMF3', 'logsky', 'IMF4',
'h3', 'h4','ML_v','ML_i','ML_k','MW_v', 'MW_i','MW_k']
elif self.filetype == 50:
self.results = results['Type',
'chi2', 'velz', 'sigma', 'logage', 'zH', 'FeH', 'a', 'C',
'N', 'Na', 'Mg', 'Si', 'K', 'Ca', 'Ti', 'V', 'Cr', 'Mn',
'Co', 'Ni', 'Cu', 'Sr', 'Ba', 'Eu', 'Teff', 'IMF1', 'IMF2',
'logfy', 'sigma2', 'velz2', 'logm7g', 'hotteff', 'loghot',
'fy_logage', 'logtrans', 'logemline_h', 'logemline_oiii',
'logemline_sii', 'logemline_ni', 'logemline_nii', 'jitter',
'IMF3', 'logsky', 'IMF4', 'ML_v', 'ML_i', 'ML_k', 'MW_v',
'MW_i', 'MW_k']

"""
Read in input data and best fit model

This isn't going to work correctly if the file
doesn't exist
"""
#try:
m = np.loadtxt('{0}.bestspec'.format(self.outfiles))
#except:
# warning = ('Do not have the *.bestspec file')
# warnings.warn(warning)
data = {}
data['wave'] = m[:,0]/(1.+self.results['velz'][5]*1e3/constants.c)
data['m_flux'] = m[:,1] # Model spectrum, normalization applied
data['d_flux'] = m[:,2] # Data spectrum
data['snr'] = m[:,3] # Including jitter and inflated errors
data['unc'] = 1/m[:,3]
if not index:
data['poly'] = m[:,4] # Polynomial used to create m_flux
data['residual'] = (m[:,1] - m[:,2])/m[:,1] * 1e2
self.spectra = data

try:
m = np.loadtxt('{0}.bestspec2'.format(self.outfiles))
model = {}
model['wave'] = m[:,0]
#model['wave'] = m[:,0]/(1.+self.results['velz'][5]*1e3/constants.c)
model['flux'] = m[:,1]
self.ext_model = model
except:
self.ext_model = None
m = np.loadtxt('{0}.bestspec'.format(self.outfiles))
data = {}
data['wave'] = m[:,0]/(1.+self.results['velz'][5]*1e3/constants.c)
data['m_flux'] = m[:,1] # Model spectrum, normalization applied
data['d_flux'] = m[:,2] # Data spectrum
data['snr'] = m[:,3] # Including jitter and inflated errors
data['unc'] = 1/m[:,3]
if not index:
data['poly'] = m[:,4] # Polynomial used to create m_flux
data['residual'] = (m[:,1] - m[:,2])/m[:,1] * 1e2
self.spectra = data

try:
m = np.loadtxt('{0}.bestspec2'.format(self.outfiles))
model = {}
model['wave'] = m[:,0]
#model['wave'] = m[:,0]/(1.+self.results['velz'][5]*1e3/constants.c)
model['flux'] = m[:,1]
self.ext_model = model
except:
self.ext_model = None

"""
Check the values of the nuisance parameters
Expand All @@ -125,7 +192,7 @@ def __init__(self, outfiles, read_mcmc=True, info=None, index=False):
# warnings.warn(warning.format(self.path, 'loghot',
# self.results['loghot'][0]))

def get_total_met(self):
def get_total_met(self, return_posterior=False):

zh = np.where(self.labels == 'zH')
feh = np.where(self.labels == 'FeH')
Expand All @@ -134,6 +201,9 @@ def get_total_met(self):
#Computing errors directly from the chains.
self.tmet = self.get_cls(total_met)

if return_posterior is True:
return total_met

def normalize_spectra(self):
"""
Normalize the data and model spectra
Expand All @@ -146,6 +216,13 @@ def normalize_spectra(self):
min_ = min(self.spectra['wave'])
max_ = max(self.spectra['wave'])
num = int((max_ - min_)/chunks) + 1

# If you have a model in .bestspec2, need to normalize this as well
if self.ext_model != None:
self.ext_model['flux'] = deepcopy(self.ext_model['flux'])
min_ext = min(self.ext_model['wave'])
max_ext = max(self.ext_model['wave'])
num_ext = int((max_ext - min_ext)/chunks) + 1

for i in range(num):
k = ((self.spectra['wave'] >= min_ + chunks*i) &
Expand All @@ -164,6 +241,18 @@ def normalize_spectra(self):
self.spectra['m_flux_norm'][k], 2)
poly = chebval(self.spectra['wave'][k], coeffs)
self.spectra['m_flux_norm'][k] = self.spectra['m_flux_norm'][k]/poly

# .bestspec2 model normalization
if self.ext_model != None:
for i in range(num_ext):
k_ext = ((self.ext_model['wave'] >= min_ext + chunks*i) & (self.ext_model['wave'] <= min_ext + chunks*(i+1)))

if len(self.ext_model['flux'][k_ext]) < 10:
continue

coeffs_ext = chebfit(self.ext_model['wave'][k_ext], self.ext_model['flux'][k_ext], 2)
poly_ext = chebval(self.ext_model['wave'][k_ext], coeffs_ext)
self.ext_model['flux'][k_ext] = self.ext_model['flux'][k_ext]/poly_ext

def abundance_correct(self, s07=False, b14=False, m11=True):
"""
Expand Down Expand Up @@ -444,16 +533,138 @@ def plot_posterior(self, fname):
def get_cls(self, distribution):
distribution = np.sort(np.squeeze(distribution))

num = self.nwalkers*self.nchain/self.nsample
if self.filetype == 50:
num = self.model_info['Nwalkers']*self.model_info['Nchain']
else:
num = self.model_info['Nwalkers']*self.model_info['Nchain']/self.model_info['Nsample']
## ND: New file reading code, don't need this after rehaul. Delete.
#num = self.nwalkers*self.nchain/self.nsample
lower = distribution[int(0.160*num)]
median = distribution[int(0.500*num)]
upper = distribution[int(0.840*num)]
cl95 = distribution[int(0.950*num)]
std = np.std(distribution)

return {'cl50': median, 'cl84': upper, 'cl16': lower, 'std': std}
return {'cl50': median, 'cl84': upper, 'cl16': lower,
'cl95': cl95, 'std': std}

def write_params(self):
pass

def get_mass(self):
"""
Based on the alf function.

Compute amss in stars and remnants (normalized to 1 Msun at t=0)
for different IMF parameterizations (including a non-parameterization)
"""

#try: # initial sanity check
# mlo < m2
#except:
# print("Something is wrong: ml>m2")

"""
mass = 0.0

if self.results['IMF4'] is None:
# normalize the weights so that 1 Msun formed at t=0
imf_norm = (m2**(-self.results['IMF1']+2)-mlo**(-self.results['IMF1']+2))/(-self.results['IMF1']+2) +
m2**(-self.results['IMF1']+self.results['IMF2'])*(m3**(-self.results['IMF2']+2)-
m2**(-self.results['IMF2']+2))/(-self.['IMF2']+2) +
m2**(-self.['IMF1']+self.results['IMF2'])*(imf_hi**(-imf_up+2)-m3**(-imf_up+2))/(-imf_up+2)

# mass from stars that are still alive
mass = (m2**(-self.results['IMF1']+2)-m_lo**(-self.results['IMF1']+2))/(-self.results['IMF1']+2)
if msto < m3:
mass = mass + m2**(-self.results['IMF1'+self.results['IMF2'])*
(msto**(-self.results['IMF2']+2)-m2**(-self.results['IMF2']+2))/(-self.['IMF2']+2)
else:
mass = mass +
m2**(-self.reslts['IMF1']+self.results['IMF2'])*(m3**(-self.results['IMF2']+2)-m2**(-self.results['IMF2']+2))/(-self.results['IMF2']+2)
+ m2**(-self.results['IMF1']+self.results['IMF2'])*(msto**(-imf_up+2)-m3**(-imf_up+2))/(-imf_up+2)

mass = mass/imf_norm


# Mass from blackhole remnants
# 40 < M < imf_up, leave behind a 0.5*M BH
mass = mass + 0.5*m2**(-self.results['IMF1']+self.results['IMF2'])*(imf_hi**(-imf_up+2)-bh_lim**(-imf_up+2))/(-imf_up+2)/imf_norm

# Neutron star remnants
# 8.5 < M < 40, leave behind 1.4 MSun NS
mass = mass + 1.4*m2**(-self.results['IMF1']+self.results['IMF2'])*(bh_lim**(-imf_up+1)-ns_lim**(-imf_up+1))/(-imf_up+1)/imf_norm

# White dwarf remnants
# M < 8.5, leave behind 0.077*M+0.58 WD
if msto < m3:
mass = mass + 0.48*m2**(-self.results['IMF1']+self.results['IMF2'])*
(ns_lim**(-imf_up+1)-m3**(-imf_up+1))/(-imf_up+1)/imf+norm
mass = mass + 0.48*m2**(-self.results['IMF1']+self.results['IMF2'])*
(m3**(-self.results['IMF2']+1)-msto**(-self.results['IMF2']+1))/(-self.results['IMF2']+1)/imf_norm
mass = mass + 0.077*m2**(-self.results['IMF1']+self.results['IMF2'])*(ns_lim**(-imf_up+2)-m3**(-imf_up+2))/(-imf_up+2)/imfnorm
mass = mass +
0.077*m2**(-self.results['IMF1']+self.results['IMF2'])*(m3**(-self.results['IMF2']+2)-mtso**(-self.results['IMF2']+2))/(-self.results['IMF2']+2)/imf_norm
else:
mass = mass + 0.48*m2**(-self.results['IMF1']+self.results['IMF2'])*(ns_lim**(-imf_up+1)-mtso**(-imf_up+1))/(-imf_up+1)/imfnorm
mass = mass + 0.077*m2**(-self.results['IMF1']+self.results['IMF2'])*(ns_lim**(-imf_up+2)-msto**(-imf_up+2))/(-imf_up+2)/imfnorm


#else: Non-parametric IMF fit


#IF (PRESENT(timfnorm)) timfnorm = imfnorm

return norm
"""
return None

def get_m2l(self, filters):
"""
Based on the alf code.

Computes mass-to-light ratios (AB mags) over user given filters
"""

# convert to the "proper" units
aspec = self.data['d_flux']*lsun/1E6*self.data['wave']**2/clight/1E8/4/mypi/pc2cm**2

if self.model_info['mwimf'] == 1: # model was forced to fit for a MW (Kroupa) IMF
mass = get_mass(imf_lo, msto, kroupa_imf1, kroupa_imf2, kroupa_imf3)
"""
else:
if imf_type == 0: # single powerlaw IMF with a fixed lower-mass cutoff
mass = get_mass(imf_lo, msto, self.results['IMF1'], self.results['IMF1'], kroupa_imf3)
elif imf_type == 1: # double powerlaw with a fixed lower-mass cutoff
mass = get_mass(imf_lo, msto, self.results['IMF1'], self.results['IMF2'], kroupa_imf3)
elif imf_type == 2: # single powerlaw index with variable lower-mass cutoff
mass = get_mass(self.results['IMF3'], msto, self.results['IMF1'], self.results['IMF1'], kroupa_imf3)
elif imf_type == 3: # double powerlaw index with variable lower-mass cutoff
mass = get_mass(self.results['IMF3'], msto, self.results['IMF1'], self.results['IMF2'], kroupa_imf3)
elif imf_type == 4: # non-parametric IMF for 0.08-1.0 Msun; Salpeter slope at >1 Msun
mass = get_mass(imf_lo, msto, self.results['IMF1'], self.results['IMF2'], kroupa_imf3,
self.results['IMF3'], self.results['IMF4'])
"""


m2l = {}
for name, bandpass in filters:
mag = tsum(wave, aspec*f/lam)
if mag < 0.0:
m2l[name] == 0.0
else:
mag = -2.5*np.log10(mag) - 48.60
m2l[name] = mass//10**(2./5*(magsun[name]-mag))
# if m2l > 100 then m2l = 0.0

return m2l

if __name__=='__main__':
pass