Pulled develop into feature/validation_tables to work with a more upd…

…ated version

Metallicity_Stack_Commons/__init__.py

@@ -5,11 +5,15 @@
 import getpass
 import numpy as np
 
-version = "0.5.0"
+version = "0.6.0"
 
 lambda0   = [3726.18, 4101.73, 4340.46, 4363.21, 4861.32, 4958.91, 5006.84]
 line_type = ['Oxy2', 'Balmer', 'Balmer', 'Single', 'Balmer', 'Single', 'Single']
-line_name = ['OII_3727', 'HDELTA', 'HGAMMA', 'OIII_4363', 'HBETA', 'OIII_4958', 'OIII_5007']
+line_name = ['OII_3727', 'HDELTA', 'HGAMMA', 'OIII_4363', 'HBETA', 'OIII_4958',
+             'OIII_5007']
+
+line_name_short = {"OII": line_name[0], "4363": line_name[3],
+                   "HB": line_name[4], "OIII": line_name[-1]}
 
 fitting_lines_dict = {"lambda0": lambda0, "line_type": line_type, "line_name": line_name}
 

Metallicity_Stack_Commons/attenuation.py → Metallicity_Stack_Commons/analysis/attenuation.py

@@ -2,7 +2,7 @@
 from astropy.table import Table
 import numpy as np
 
-from . import k_dict
+from .. import k_dict
 
 HgHb_CaseB = 0.468 # Hg/Hb ratio for zero reddening
 

Metallicity_Stack_Commons/analysis/composite_indv_detect.py

@@ -5,14 +5,19 @@
 from astropy.io import ascii as asc
 from astropy.table import Table
 
-from ..temp_metallicity_calc import metallicity_calculation
+from .temp_metallicity_calc import metallicity_calculation
 from .. import OIII_r
 from ..column_names import bin_names0, indv_names0, temp_metal_names0
 from ..column_names import filename_dict
 
 ID_name = indv_names0[0]
 bin_ID_name = bin_names0[0]
 
+logR23_name = indv_names0[1]
+logO32_name = indv_names0[2]
+two_beta_name = indv_names0[5]
+three_beta_name = indv_names0[6]
+
 
 def main(fitspath, dataset, revised=False, det3=True):
     """
@@ -57,6 +62,13 @@ def main(fitspath, dataset, revised=False, det3=True):
     bin_id = composite_table['bin_ID'].data
     bin_temp = composite_table['T_e'].data
 
+    # Read in validation table
+    valid_file = join(fitspath, dataset, filename_dict['bin_valid'])
+    if not exists(valid_file):
+        print("ERROR: File not found! "+valid_file)
+        return
+    valid_table = asc.read(valid_file)
+
     # Define [indv_em_line_file]
     indv_em_line_file = join(fitspath, dataset, filename_dict['indv_prop'])
     if not exists(indv_em_line_file):
@@ -78,34 +90,36 @@ def main(fitspath, dataset, revised=False, det3=True):
     # Populate composite temperature for individual galaxies
     adopted_temp = np.zeros(len(indv_em_line_table))
     bin_id_indv = np.zeros(len(indv_em_line_table))
-    for comp_bin, comp_temp in zip(bin_id, bin_temp):
+    detect_indv = np.zeros(len(indv_em_line_table))
+    for comp_bin, comp_temp, detect in zip(bin_id, bin_temp, valid_table['Detection']):
         bin_idx = np.where(indv_bin_info_table['bin_ID'].data == comp_bin)[0]
         adopted_temp[bin_idx] = comp_temp
         bin_id_indv[bin_idx]  = comp_bin
+        detect_indv[bin_idx]  = detect
+
+    O2 = indv_em_line_table['OII_3727_Flux_Gaussian'].data   # [OII]3726,3728 fluxes
+    O3 = indv_em_line_table['OIII_5007_Flux_Gaussian'].data  # [OIII]5007 fluxes
+    O3 = O3 * (1+1/OIII_r)  # Scale to include OIII4959; Assume 3.1:1 ratio
+    Hb = indv_em_line_table['HBETA_Flux_Gaussian'].data      # H-beta fluxes
 
-    O2 = indv_em_line_table['OII_3727_Flux_Gaussian'].data            # [OII]3726,3728 fluxes
-    O3 = indv_em_line_table['OIII_5007_Flux_Gaussian'].data * OIII_r  # [OIII]4959,5007 fluxes (Assume 3.1:1 ratio)
-    Hb = indv_em_line_table['HBETA_Flux_Gaussian'].data               # H-beta fluxes
+    two_beta = O2/Hb
+    three_beta = O3/Hb
+    logR23 = np.log10(two_beta + three_beta)
+    logO32 = np.log10(O3/O2)
 
-    if det3:
-        com_O_log, metal_dict = metallicity_calculation(adopted_temp, O2/Hb, O3/Hb)
+    if not det3:
+        metal_dict = metallicity_calculation(adopted_temp, O2/Hb, O3/Hb)
     else:
-        det3 = np.where((indv_bin_info_table['Detection'] == 1.0) | (indv_bin_info_table['Detection'] == 0.5))[0]
-        temp_com_O_log, temp_metal_dict = \
-            metallicity_calculation(adopted_temp[det3], O2[det3]/Hb[det3],
-                                    O3[det3]/Hb[det3])
-        com_O_log = np.zeros(len(indv_em_line_table))
-        com_O_log[det3] = temp_com_O_log
-
-        metal_dict = dict()
-        for key0 in temp_metal_dict.keys():
-            metal_dict[key0] = np.zeros(len(indv_em_line_table))
-            metal_dict[key0][det3] = temp_metal_dict[key0]
+        det3_idx = np.where((detect_indv == 1.0) | (detect_indv == 0.5))[0]
+        metal_dict = \
+            metallicity_calculation(adopted_temp, O2/Hb, O3/Hb, det3=det3_idx)
 
     # Define [indv_derived_prop_table] to include ID, bin_ID, composite T_e,
     # and 12+log(O/H)
-    arr0 = [indv_em_line_table[ID_name], bin_id_indv, adopted_temp, com_O_log]
-    names0 = [ID_name, bin_ID_name] + temp_metal_names0[:2]
+    arr0 = [indv_em_line_table[ID_name], bin_id_indv, logR23, logO32,
+            two_beta, three_beta, adopted_temp]
+    names0 = [ID_name, bin_ID_name, logR23_name, logO32_name, two_beta_name,
+              three_beta_name, temp_metal_names0[0]]
 
     # Include other metallicities
     arr0 += list(metal_dict.values())

Metallicity_Stack_Commons/analysis/error_prop.py

@@ -0,0 +1,161 @@
+from os.path import join, exists
+
+from chun_codes import random_pdf, compute_onesig_pdf
+from . import line_name
+from astropy.io import ascii as asc
+import numpy as np
+
+from ..column_names import filename_dict, temp_metal_names0, npz_filename_dict
+from .ratios import error_prop_flux_ratios
+from .temp_metallicity_calc import temp_calculation, metallicity_calculation
+
+
+def write_npz(path, npz_files, dict_list):
+    """
+    Purpose:
+      Write numpy files with provided dictionaries
+
+    :param path: str - prefix for filename output
+    :param npz_files: list - contains npz file names
+    :param dict_list: list - contains dictionaries for each corresponding npz file
+
+    :return: Write npz files
+    """
+    for file, dict_input in zip(npz_files, dict_list):
+        npz_outfile = join(path, file)
+        if exists(npz_outfile):
+            print("Overwriting : "+npz_outfile)
+        else:
+            print("Writing : "+npz_outfile)
+        np.savez(npz_outfile, **dict_input)
+
+
+def fluxes_derived_prop(path, binned_data=True):
+    """
+    Purpose:
+      Use measurements and their uncertainties to perform a randomization
+      approach. The randomization is performed on individual emission lines.
+      It carries that information to derived flux ratios and then
+      determines electron temperature and metallicity
+
+    :param path: str of full path
+    :param binned_data: bool for whether to analysis binned data. Default: True
+
+    :return:
+    """
+
+    # Define files to read in for binned data
+    if binned_data:
+        flux_file = join(path, filename_dict['bin_fit'])
+        prop_file = join(path, filename_dict['bin_derived_prop'])
+        verify_file = join(path, filename_dict['bin_valid'])
+
+    flux_tab0 = asc.read(flux_file)
+    prop_tab0 = asc.read(prop_file)
+
+    if binned_data:
+        verify_tab = asc.read(verify_file)
+        detection = verify_tab['Detection'].data
+
+        # For now we are only considering those with reliable detection and
+        # excluding those with reliable non-detections (detect = 0.5)
+        detect_idx = np.where((detection == 1))[0]
+
+        ID = verify_tab['bin_ID'].data
+        ID_detect = ID[detect_idx]
+        print(ID_detect)
+
+        flux_tab = flux_tab0[detect_idx]
+        prop_tab = prop_tab0[detect_idx]
+
+    flux_cols     = [str0+'_Flux_Gaussian' for str0 in line_name]
+    flux_rms_cols = [str0+'_RMS' for str0 in line_name]
+
+    #
+    # EMISSION-LINE SECTION
+    #
+
+    # Initialize dictionaries
+    flux_pdf_dict = dict()
+    flux_peak = dict()
+    flux_lowhigh = dict()
+
+    # Randomization for emission-line fluxes
+    for aa, flux, rms in zip(range(len(flux_cols)), flux_cols, flux_rms_cols):
+        flux_pdf = random_pdf(flux_tab[flux], flux_tab[rms], seed_i=aa,
+                              n_iter=1000)
+        err, peak = compute_onesig_pdf(flux_pdf, flux_tab[flux], usepeak=True)
+
+        # Fill in dictionary
+        flux_pdf_dict[line_name[aa]] = flux_pdf
+        flux_peak[line_name[aa] + '_peak'] = peak
+        flux_lowhigh[line_name[aa] + '_lowhigh_error'] = err
+
+        # Update values
+        flux_tab0[line_name[aa] + '_Flux_Gaussian'][detect_idx] = peak
+
+    # Write revised emission-line fit ASCII table
+    new_flux_file = join(path, filename_dict['bin_fit_rev'])
+    if exists(new_flux_file):
+        print("Overwriting: "+new_flux_file)
+    else:
+        print("Writing: "+new_flux_file)
+    asc.write(flux_tab0, new_flux_file, overwrite=True, format='fixed_width_two_line')
+
+    # Save flux npz files
+    npz_files = [npz_filename_dict['flux_pdf'],
+                 npz_filename_dict['flux_errors'],
+                 npz_filename_dict['flux_peak']]
+    dict_list = [flux_pdf_dict, flux_lowhigh, flux_peak]
+    write_npz(path, npz_files, dict_list)
+
+    # Obtain distributions of line ratios: logR23, logO32, two_beta, three_beta, R
+    flux_ratios_dict = error_prop_flux_ratios(flux_pdf_dict)
+
+    #
+    # DERIVED PROPERTIES SECTION
+    #
+
+    # Initialize dictionaries
+    derived_prop_pdf_dict = dict()
+    derived_prop_error_dict = dict()
+    derived_prop_peak_dict = dict()
+
+    # Calculate temperature distribution
+    Te_pdf = temp_calculation(flux_ratios_dict['R'])
+    derived_prop_pdf_dict[temp_metal_names0[0]] = Te_pdf
+
+    # Calculate metallicity distribution
+    metal_dict = metallicity_calculation(Te_pdf, flux_ratios_dict['two_beta'],
+                                         flux_ratios_dict['three_beta'])
+    derived_prop_pdf_dict.update(metal_dict)
+
+    # Loop for each derived properties (T_e, metallicity, etc.)
+    for names0 in temp_metal_names0:
+        arr0 = prop_tab[names0].data
+
+        err_prop, peak_prop = \
+            compute_onesig_pdf(derived_prop_pdf_dict[names0], arr0, usepeak=True)
+
+        # Fill in dictionary
+        derived_prop_error_dict[names0+'_lowhigh_error'] = err_prop
+        derived_prop_peak_dict[names0+'_peak'] = peak_prop
+
+        # Update values
+        prop_tab0[names0][detect_idx] = peak_prop
+
+    # Write revised properties ASCII table
+    new_prop_file = join(path, filename_dict['bin_derived_prop_rev'])
+    if exists(new_prop_file):
+        print("Overwriting: "+new_prop_file)
+    else:
+        print("Writing: "+new_prop_file)
+    asc.write(prop_tab0, new_prop_file, overwrite=True, format='fixed_width_two_line')
+
+    # Save derived properties npz files
+    npz_files = [npz_filename_dict['der_prop_pdf'],
+                 npz_filename_dict['der_prop_errors'],
+                 npz_filename_dict['der_prop_peak']]
+    dict_list = [derived_prop_pdf_dict, derived_prop_error_dict,
+                 derived_prop_peak_dict]
+    write_npz(path, npz_files, dict_list)

Metallicity_Stack_Commons/fitting.py → Metallicity_Stack_Commons/analysis/fitting.py

@@ -2,7 +2,7 @@
 from astropy.convolution import Box1DKernel, convolve
 from astropy.io import ascii as asc
 
-from . import scalefact, wavelength_dict
+from .. import scalefact, wavelength_dict
 
 con1 = wavelength_dict['OII_3729'] / wavelength_dict['OII_3726'] # Ratio of OII doublet line
 

Metallicity_Stack_Commons/analysis/ratios.py

@@ -0,0 +1,44 @@
+import numpy as np
+
+from . import line_name_short, OIII_r
+from .temp_metallicity_calc import R_calculation
+
+
+def error_prop_flux_ratios(flux_dict, EBV=None):
+    """
+    Purpose:
+      Primary code to determine a variety of line ratios based on a dictionary
+      containing emission-line fluxes
+
+    :param flux_dict: dictionary containing line ratios
+    :param EBV: array of E(B-V). Same dimensions as contents of flux_dict
+
+    :return flux_ratios_dict: dictionary containing flux ratios
+    """
+
+    # Retrieve emission line fluxes
+    OII  = flux_dict[line_name_short['OII']]
+    OIII = flux_dict[line_name_short['OIII']]
+    Hb   = flux_dict[line_name_short['HB']]
+    OIII4363 = flux_dict[line_name_short['4363']]
+
+    # Define flux ratios
+    two_beta = OII/Hb
+    three_beta = (1+1/OIII_r) * OIII/Hb
+    logR23 = np.log10(two_beta + three_beta)
+    logO32 = np.log10((1+1/OIII_r) * OIII/OII)
+
+    # Define dictionary of flux ratios
+    flux_ratios_dict = dict()
+    flux_ratios_dict['two_beta'] = two_beta
+    flux_ratios_dict['three_beta'] = three_beta
+    flux_ratios_dict['logR23'] = logR23
+    flux_ratios_dict['logO32'] = logO32
+
+    if EBV is None:
+        print("Not applying dust attenuation correction")
+        EBV = np.zeros(OIII4363.shape)
+
+    flux_ratios_dict['R'] = R_calculation(OIII4363, OIII, EBV)
+
+    return flux_ratios_dict