diff --git a/CHANGES.rst b/CHANGES.rst index 7d75f0a8..c94f28eb 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -12,4 +12,4 @@ ================== - IDL versions -- http://tir.astro.utoledo.edu/jdsmith/research/pahfit.php +- https://github.com/PAHFIT/pahfit_classic diff --git a/README.rst b/README.rst index dc92abed..1a02367c 100644 --- a/README.rst +++ b/README.rst @@ -2,7 +2,7 @@ PAHFIT ====== -``PAHFIT`` is a decomposition model and tool for astronomical infrared spectra, focusing on dust and gas emission features from the interstellar medium (see `Smith, J.D.T., Draine B.T., et al., 2007, ApJ, 656, 770 `_). +``PAHFIT`` is a decomposition model and tool for astronomical infrared spectra, focusing on dust and gas emission features from the interstellar medium (see `Smith, J.D.T., Draine B.T., et al., 2007, ApJ, 656, 770 `_). This package provides an updated python implementation of ``PAHFIT``. While the original versions of ``PAHFIT`` (``v1.x``) were written in IDL and focused mainly on Spitzer/IRS spectroscopic observations, the newer python-based versions (``>=v2.0``) will expand instrument coverage to other existing (e.g., AKARI) and planned (e.g., JWST) facilities, and will offer a more flexible modeling framework suitable for modeling a wider range of astrophysical sources. diff --git a/docs/index.rst b/docs/index.rst index 21352c0d..e6072367 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -14,7 +14,7 @@ and a more flexible modeling framework suitable for modeling a wider range of astrophysical sources. For details for the IDL version of PAHFIT see -`Smith, J.D.T., Draine B.T., et al., 2007, ApJ, 656, 770 `_. +`Smith, J.D.T., Draine B.T., et al., 2007, ApJ, 656, 770 `_. This package is potentially an `astropy affiliated package `_ @@ -77,7 +77,7 @@ contributors who will abide by the `Python Software Foundation Code of Conduct `Astropy`_. The following pages will help you get started with contributing fixes, code, or documentation (no git or GitHub experience necessary): -* `How to make a code contribution `_ +* `How to make a code contribution `_ * `Coding Guidelines `_ @@ -92,10 +92,6 @@ contributors page on Github Reference API ============= -Base class for PAHFIT - -.. automodapi:: pahfit.base - Component models not provided by astropy.models. -.. automodapi:: pahfit.component_models +.. automodapi:: pahfit.fitters.ap_components diff --git a/docs/version_description/version_description.rst b/docs/version_description/version_description.rst index 609592da..882f12a5 100644 --- a/docs/version_description/version_description.rst +++ b/docs/version_description/version_description.rst @@ -5,7 +5,7 @@ Version Description v1.0 - v1.4 ------------ -The original IDL version of PAHFIT can be found in: `http://tir.astro.utoledo.edu/jdsmith/research/pahfit.php `_. For more details and background of ``PAHFIT``, see the `background `_ page. +The original IDL version of PAHFIT can be found in: `https://github.com/PAHFIT/pahfit_classic `_. For more details and background of ``PAHFIT``, see the `background `_ page. v2.0 ------------ diff --git a/pahfit/base.py b/pahfit/base.py deleted file mode 100644 index a0ff47e9..00000000 --- a/pahfit/base.py +++ /dev/null @@ -1,758 +0,0 @@ -import numpy as np -import matplotlib as mpl - -from pahfit.instrument import within_segment, fwhm -from pahfit.errors import PAHFITModelError -from pahfit.component_models import ( - BlackBody1D, - ModifiedBlackBody1D, - S07_attenuation, - att_Drude1D, -) -from astropy.modeling.physical_models import Drude1D -from astropy.modeling.functional_models import Gaussian1D - -__all__ = ["PAHFITBase"] - - -def _ingest_limits(min_vals, max_vals): - """ - Ingest the limits read from yaml file and generate the appropriate - internal format (list of tuples). Needed to handle the case - where a limit is not desired as numpy arrays cannot have elements - of None, instead a value of nan is used. - - Limits that are not set are designated as 'nan' in files and - these are changed to the python None to be compatible with - the astropy.modeling convention. - - Parameters - ---------- - min_vals, - max_vals : numpy.array (masked arrays) - min/max values of the limits for a parameter - nan designates no limit - - Returns - ------- - plimits : list of tuples - tuples give the min/max limits for a parameter - """ - plimits = [] - mask_min = min_vals.mask - data_min = min_vals.data - mask_max = max_vals.mask - data_max = max_vals.data - - mask_min_ind = np.where(np.logical_not(mask_min))[0] - mask_max_ind = np.where(np.logical_not(mask_max))[0] - - min_vals = np.zeros(len(mask_min)) - min_vals[mask_min_ind] = data_min[mask_min_ind] - - max_vals = np.zeros(len(mask_max)) - max_vals[mask_max_ind] = data_max[mask_max_ind] - - plimits = [] - for cmin, cmax in zip(min_vals, max_vals): - if np.isinf(cmin): - cmin = None - if np.isinf(cmax): - cmax = None - plimits.append((cmin, cmax)) - - return plimits - - -def _ingest_fixed(fixed_vals): - """ - Ingest the fixed value read from a file and generate the appropriate - internal format (list of booleans). Since this information is indirectly - hidden in the parameter of a feature, this function is needed to - extract that. - - Parameters - ---------- - min_vals : numpy.array (masked array) - fixed designations - - Returns - ------- - pfixed : list (boolean) - True/False designation for parameters - """ - mask_false_ind = np.where(np.logical_not(fixed_vals.mask))[0] - fixed_vals = ["True"] * len(fixed_vals.mask) - for i in range(0, len(mask_false_ind)): - ind = mask_false_ind[i] - fixed_vals[ind] = "False" - - pfixed = [] - for cfixed in fixed_vals: - if cfixed == "True": - cfixed = True - if cfixed == "False": - cfixed = False - pfixed.append(cfixed) - - return pfixed - - -class PAHFITBase: - """ - Old implementation. Some functions are still used by the new Model - class. The unused functionality has been removed. - - Construct that is still used for now - - param_info: tuple of dicts called (bb_info, df_info, h2_info, ion_info, abs_info, att_info) - The dictionaries contain info for each type of component. Each - component of the dictonaries is a vector. - bb_info - - dict with {name, temps, temps_limits, temps_fixed, - amps, amps_limits, amps_fixed}, each a vector - dust_features, h2_features, ion_features - - dict with {name amps, amps_limits, amps_fixed, - x_0, x_0_limits, x_0_fixed, fwhms, fwhms_limits, fwhm_fixed}. - """ - @staticmethod - def model_from_param_info(param_info): - # setup the model - bb_info = param_info[0] - dust_features = param_info[1] - h2_features = param_info[2] - ion_features = param_info[3] - abs_info = param_info[4] - att_info = param_info[5] - - model = None - if bb_info is not None: - bbs = [] - for k in range(len(bb_info["names"])): - BBClass = ModifiedBlackBody1D if bb_info["modified"][ - k] else BlackBody1D - bbs.append( - BBClass( - name=bb_info["names"][k], - temperature=bb_info["temps"][k], - amplitude=bb_info["amps"][k], - bounds={ - "temperature": bb_info["temps_limits"][k], - "amplitude": bb_info["amps_limits"][k], - }, - fixed={ - "temperature": bb_info["temps_fixed"][k], - "amplitude": bb_info["amps_fixed"][k], - }, - ) - ) - model = sum(bbs[1:], bbs[0]) - - if dust_features is not None: - df = [] - for k in range(len(dust_features["names"])): - df.append( - Drude1D( - name=dust_features["names"][k], - amplitude=dust_features["amps"][k], - x_0=dust_features["x_0"][k], - fwhm=dust_features["fwhms"][k], - bounds={ - "amplitude": dust_features["amps_limits"][k], - "x_0": dust_features["x_0_limits"][k], - "fwhm": dust_features["fwhms_limits"][k], - }, - fixed={ - "amplitude": dust_features["amps_fixed"][k], - "x_0": dust_features["x_0_fixed"][k], - "fwhm": dust_features["fwhms_fixed"][k], - }, - )) - - df = sum(df[1:], df[0]) - if model: - model += df - else: - model = df - - if h2_features is not None: - h2 = [] - for k in range(len(h2_features["names"])): - h2.append( - Gaussian1D( - name=h2_features["names"][k], - amplitude=h2_features["amps"][k], - mean=h2_features["x_0"][k], - stddev=h2_features["fwhms"][k] / 2.355, - bounds={ - "amplitude": - h2_features["amps_limits"][k], - "mean": - h2_features["x_0_limits"][k], - "stddev": ( - h2_features["fwhms"][k] * 0.9 / 2.355, - h2_features["fwhms"][k] * 1.1 / 2.355, - ), - }, - fixed={ - "amplitude": h2_features["amps_fixed"][k], - "mean": h2_features["x_0_fixed"][k], - "stddev": h2_features["fwhms_fixed"][k], - }, - )) - h2 = sum(h2[1:], h2[0]) - if model: - model += h2 - else: - model = h2 - - if ion_features is not None: - ions = [] - for k in range(len(ion_features["names"])): - ions.append( - Gaussian1D( - name=ion_features["names"][k], - amplitude=ion_features["amps"][k], - mean=ion_features["x_0"][k], - stddev=ion_features["fwhms"][k] / 2.355, - bounds={ - "amplitude": - ion_features["amps_limits"][k], - "mean": - ion_features["x_0_limits"][k], - "stddev": ( - ion_features["fwhms"][k] * 0.9 / 2.355, - ion_features["fwhms"][k] * 1.1 / 2.355, - ), - }, - fixed={ - "amplitude": ion_features["amps_fixed"][k], - "mean": ion_features["x_0_fixed"][k], - "stddev": ion_features["fwhms_fixed"][k], - }, - )) - ions = sum(ions[1:], ions[0]) - if model: - model += ions - else: - model = ions - - # add additional att components to the model if necessary - if not model: - raise PAHFITModelError("No model components found") - - if abs_info is not None: - for k in range(len(abs_info["names"])): - model *= att_Drude1D( - name=abs_info["names"][k], - tau=abs_info["amps"][k], - x_0=abs_info["x_0"][k], - fwhm=abs_info["fwhms"][k], - bounds={ - "tau": abs_info["amps_limits"][k], - "fwhm": abs_info["fwhms_limits"][k], - }, - fixed={"x_0": abs_info["x_0_fixed"][k]}, - ) - - if att_info is not None: - model *= S07_attenuation( - name=att_info["name"], - tau_sil=att_info["tau_sil"], - bounds={"tau_sil": att_info["tau_sil_limits"]}, - fixed={"tau_sil": att_info["tau_sil_fixed"]}, - ) - - return model - - @staticmethod - def plot(axs, x, y, yerr, model, model_samples=1000, scalefac_resid=2): - """ - Plot model using axis object. - - Parameters - ---------- - axs : matplotlib.axis objects - where to put the plot - x : floats - wavelength points - y : floats - observed spectrum - yerr: floats - observed spectrum uncertainties - model : PAHFITBase model (astropy modeling CompoundModel) - model giving all the components and parameters - model_samples : int - Total number of wavelength points to allocate to the model display - scalefac_resid : float - Factor multiplying the standard deviation of the residuals to adjust plot limits - """ - # remove units if they are present - if hasattr(x, "value"): - x = x.value - if hasattr(y, "value"): - y = y.value - if hasattr(yerr, "value"): - yerr = yerr.value - - # Fine x samples for model fit - x_mod = np.logspace(np.log10(min(x)), np.log10(max(x)), model_samples) - - # spectrum and best fit model - ax = axs[0] - ax.set_yscale("linear") - ax.set_xscale("log") - ax.minorticks_on() - ax.tick_params(axis="both", - which="major", - top="on", - right="on", - direction="in", - length=10) - ax.tick_params(axis="both", - which="minor", - top="on", - right="on", - direction="in", - length=5) - - ax_att = ax.twinx() # axis for plotting the extinction curve - ax_att.tick_params(which="minor", direction="in", length=5) - ax_att.tick_params(which="major", direction="in", length=10) - ax_att.minorticks_on() - - # get the extinction model (probably a better way to do this) - ext_model = None - for cmodel in model: - if isinstance(cmodel, S07_attenuation): - ext_model = cmodel(x_mod) - - # get additional extinction components that can be - # characterized by functional forms (Drude profile in this case) - for cmodel in model: - if isinstance(cmodel, att_Drude1D): - if ext_model is not None: - ext_model *= cmodel(x_mod) - else: - ext_model = cmodel(x_mod) - ax_att.plot(x_mod, ext_model, "k--", alpha=0.5) - ax_att.set_ylabel("Attenuation") - ax_att.set_ylim(0, 1.1) - - # Define legend lines - Leg_lines = [ - mpl.lines.Line2D([0], [0], color="k", linestyle="--", lw=2), - mpl.lines.Line2D([0], [0], color="#FE6100", lw=2), - mpl.lines.Line2D([0], [0], color="#648FFF", lw=2, alpha=0.5), - mpl.lines.Line2D([0], [0], color="#DC267F", lw=2, alpha=0.5), - mpl.lines.Line2D([0], [0], color="#785EF0", lw=2, alpha=1), - mpl.lines.Line2D([0], [0], color="#FFB000", lw=2, alpha=0.5), - ] - - # create the continum compound model (base for plotting lines) - cont_components = [] - - for cmodel in model: - if isinstance(cmodel, BlackBody1D): - cont_components.append(cmodel) - # plot as we go - ax.plot(x_mod, - cmodel(x_mod) * ext_model / x_mod, - "#FFB000", - alpha=0.5) - cont_model = cont_components[0] - for cmodel in cont_components[1:]: - cont_model += cmodel - cont_y = cont_model(x_mod) - - # now plot the dust bands and lines - for cmodel in model: - if isinstance(cmodel, Gaussian1D): - ax.plot( - x_mod, - (cont_y + cmodel(x_mod)) * ext_model / x_mod, - "#DC267F", - alpha=0.5, - ) - if isinstance(cmodel, Drude1D): - ax.plot( - x_mod, - (cont_y + cmodel(x_mod)) * ext_model / x_mod, - "#648FFF", - alpha=0.5, - ) - - ax.plot(x_mod, cont_y * ext_model / x_mod, "#785EF0", alpha=1) - - ax.plot(x_mod, model(x_mod) / x_mod, "#FE6100", alpha=1) - ax.errorbar( - x, - y / x, - yerr=yerr / x, - fmt="o", - markeredgecolor="k", - markerfacecolor="none", - ecolor="k", - elinewidth=0.2, - capsize=0.5, - markersize=6, - ) - - ax.set_ylim(0) - ax.set_ylabel(r"$\nu F_{\nu}$") - - ax.legend( - Leg_lines, - [ - "S07_attenuation", - "Spectrum Fit", - "Dust Features", - r"Atomic and $H_2$ Lines", - "Total Continuum Emissions", - "Continuum Components", - ], - prop={"size": 10}, - loc="best", - facecolor="white", - framealpha=1, - ncol=3, - ) - - # residuals, lower sub-figure - res = (y - model(x)) / x - std = np.std(res) - ax = axs[1] - - ax.set_yscale("linear") - ax.set_xscale("log") - ax.tick_params(axis="both", - which="major", - top="on", - right="on", - direction="in", - length=10) - ax.tick_params(axis="both", - which="minor", - top="on", - right="on", - direction="in", - length=5) - ax.minorticks_on() - - # Custom X axis ticks - ax.xaxis.set_ticks( - [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 20, 25, 30, 40]) - - ax.axhline(0, linestyle="--", color="gray", zorder=0) - ax.plot(x, res, "ko-", fillstyle="none", zorder=1) - ax.set_ylim(-scalefac_resid * std, scalefac_resid * std) - ax.set_xlim(np.floor(np.amin(x)), np.ceil(np.amax(x))) - ax.set_xlabel(r"$\lambda$ [$\mu m$]") - ax.set_ylabel("Residuals [%]") - - # scalar x-axis marks - ax.xaxis.set_minor_formatter(mpl.ticker.ScalarFormatter()) - ax.xaxis.set_major_formatter(mpl.ticker.ScalarFormatter()) - - @staticmethod - def update_dictionary(feature_dict, instrumentname, update_fwhms=False, redshift=0): - """ - Update parameter dictionary based on the instrument name. - Based on the instrument name, this function removes the - features outside of the wavelength range and - updates the FWHMs of the lines. - - - Parameters - ---------- - feature_dict : dictionary - Dictionary created by reading in a science pack. - - instrumentname : string - Name of the instrument with which the input spectrum - is observed. - - update_fwhms = Boolean - True for h2_info and ion_info - False for df_info - - Returns - ------- - updated feature_dict - """ - if feature_dict is None: - return None - - # convert from physical feature, to observed wavelength - def redshifted_waves(): - return feature_dict["x_0"] * (1 + redshift) - - ind = np.nonzero(within_segment(redshifted_waves(), instrumentname))[0] - - # select the valid entries in these arrays - array_keys = ("x_0", "amps", "fwhms", "names") - new_values_1 = {key: feature_dict[key][ind] for key in array_keys} - - # these are lists instead - list_keys = ( - "amps_fixed", - "fwhms_fixed", - "x_0_fixed", - "x_0_limits", - "amps_limits", - "fwhms_limits", - ) - new_values_2 = { - key: [feature_dict[key][i] for i in ind] for key in list_keys - } - - feature_dict.update(new_values_1) - feature_dict.update(new_values_2) - - if len(feature_dict['names']) == 0: - # if we removed all the things, be careful here. Setting to - # None should make the model construction function behave - # normally. - feature_dict = None - return feature_dict - - if update_fwhms: - # observe the lines at the redshifted wavelength - fwhm_min_max = fwhm(instrumentname, redshifted_waves(), as_bounded=True) - # shift the observed fwhm back to the rest frame (where the - # observed data will be moved, and its features will become - # narrower) - fwhm_min_max /= (1 + redshift) - # For astropy a numpy.bool does not work for the 'fixed' - # parameter. It needs to be a regular bool. Doing tolist() - # instead of using the array mask directly solves this. - feature_dict.update( - { - "fwhms": fwhm_min_max[:, 0], - # masked means there is no min/max, i.e. they need to be fixed - "fwhms_fixed": fwhm_min_max[:, 1].mask.tolist(), - "fwhms_limits": fwhm_min_max[:, 1:].tolist(), - } - ) - - return feature_dict - - @staticmethod - def parse_table(pack_table): - """ - Load the model parameters from a Table - - Parameters - ---------- - pack_table : Table - Table created by reading in a science pack. - - Returns - ------- - readout : tuple - Tuple containing dictionaries of all components from the - input Table. Can be used to create PAHFITBase instance using - param_info argument. Dictionary in tuple is None if no - components of that type were specified. - """ - # Getting indices for the different components - bb_ind = np.where((pack_table["kind"] == "starlight") - | (pack_table["kind"] == "dust_continuum"))[0] - df_ind = np.where(pack_table["kind"] == "dust_feature")[0] - ga_ind = np.where(pack_table["kind"] == "line")[0] - at_ind = np.where(pack_table["kind"] == "attenuation")[0] - ab_ind = np.where(pack_table["kind"] == "absorption")[0] - - # now split the gas emission lines between H2 and ions - names = [str(i) for i in pack_table["name"][ga_ind]] - if len(names) > 0: - # this has trouble with empty list - h2_temp = np.char.find(names, "H2") >= 0 - ion_temp = np.char.find(names, "H2") == -1 - h2_ind = ga_ind[h2_temp] - ion_ind = ga_ind[ion_temp] - else: - h2_ind = [] - ion_ind = [] - # the rest works fine with empty list - - # Creating the blackbody dict - bb_info = None - if len(bb_ind) > 0: - bb_info = { - "names": - np.array(pack_table["name"][bb_ind].data), - "temps": - np.array(pack_table["temperature"][:, 0][bb_ind].data), - "temps_limits": - _ingest_limits( - pack_table["temperature"][:, 1][bb_ind], - pack_table["temperature"][:, 2][bb_ind], - ), - "temps_fixed": - _ingest_fixed(pack_table["temperature"][:, 1][bb_ind]), - "amps": - np.array(pack_table["tau"][:, 0][bb_ind].data), - "amps_limits": - _ingest_limits( - pack_table["tau"][:, 1][bb_ind], - pack_table["tau"][:, 2][bb_ind], - ), - "amps_fixed": - _ingest_fixed(pack_table["tau"][:, 1][bb_ind]), - "modified": - np.array(pack_table["kind"][bb_ind] == "dust_continuum"), - } - - # Creating the dust_features dict - df_info = None - if len(df_ind) > 0: - df_info = { - "names": - np.array(pack_table["name"][df_ind].data), - "x_0": - np.array(pack_table["wavelength"][:, 0][df_ind].data), - "x_0_limits": - _ingest_limits( - pack_table["wavelength"][:, 1][df_ind], - pack_table["wavelength"][:, 2][df_ind], - ), - "x_0_fixed": - _ingest_fixed(pack_table["wavelength"][:, 1][df_ind]), - "amps": - np.array(pack_table["power"][:, 0][df_ind].data), - "amps_limits": - _ingest_limits( - pack_table["power"][:, 1][df_ind], - pack_table["power"][:, 2][df_ind], - ), - "amps_fixed": - _ingest_fixed(pack_table["power"][:, 1][df_ind]), - "fwhms": - np.array(pack_table["fwhm"][:, 0][df_ind].data), - "fwhms_limits": - _ingest_limits( - pack_table["fwhm"][:, 1][df_ind], - pack_table["fwhm"][:, 2][df_ind], - ), - "fwhms_fixed": - _ingest_fixed(pack_table["fwhm"][:, 1][df_ind]), - } - - # Creating the H2 dict - h2_info = None - if len(h2_ind) > 0: - h2_info = { - "names": - np.array(pack_table["name"][h2_ind].data), - "x_0": - np.array(pack_table["wavelength"][:, 0][h2_ind].data), - "x_0_limits": - _ingest_limits( - pack_table["wavelength"][:, 1][h2_ind], - pack_table["wavelength"][:, 2][h2_ind], - ), - "x_0_fixed": - _ingest_fixed(pack_table["wavelength"][:, 1][h2_ind]), - "amps": - np.array(pack_table["power"][:, 0][h2_ind].data), - "amps_limits": - _ingest_limits( - pack_table["power"][:, 1][h2_ind], - pack_table["power"][:, 2][h2_ind], - ), - "amps_fixed": - _ingest_fixed(pack_table["power"][:, 1][h2_ind]), - "fwhms": - np.array(pack_table["fwhm"][:, 0][h2_ind].data), - "fwhms_limits": - _ingest_limits( - pack_table["fwhm"][:, 1][h2_ind], - pack_table["fwhm"][:, 2][h2_ind], - ), - "fwhms_fixed": - _ingest_fixed(pack_table["fwhm"][:, 1][h2_ind]), - } - - # Creating the ion dict - ion_info = None - if len(ion_ind) > 0: - ion_info = { - "names": - np.array(pack_table["name"][ion_ind].data), - "x_0": - np.array(pack_table["wavelength"][:, 0][ion_ind].data), - "x_0_limits": - _ingest_limits( - pack_table["wavelength"][:, 1][ion_ind], - pack_table["wavelength"][:, 2][ion_ind], - ), - "x_0_fixed": - _ingest_fixed(pack_table["wavelength"][:, 1][ion_ind]), - "amps": - np.array(pack_table["power"][:, 0][ion_ind].data), - "amps_limits": - _ingest_limits( - pack_table["power"][:, 1][ion_ind], - pack_table["power"][:, 2][ion_ind], - ), - "amps_fixed": - _ingest_fixed(pack_table["power"][:, 1][ion_ind]), - "fwhms": - np.array(pack_table["fwhm"][:, 0][ion_ind].data), - "fwhms_limits": - _ingest_limits( - pack_table["fwhm"][:, 1][ion_ind].data, - pack_table["fwhm"][:, 2][ion_ind].data, - ), - "fwhms_fixed": - _ingest_fixed(pack_table["fwhm"][:, 1][ion_ind].data), - } - - # Create the attenuation dict (could be absorption drudes - # and S07 model) - abs_info = None - if len(ab_ind) > 0: - abs_info = { - "names": - np.array(pack_table["name"][at_ind].data), - "x_0": - np.array(pack_table["wavelength"][:, 0][at_ind].data), - "x_0_limits": - _ingest_limits( - pack_table["wavelength"][:, 1][at_ind], - pack_table["wavelength"][:, 2][at_ind], - ), - "x_0_fixed": - _ingest_fixed(pack_table["wavelength"][:, 1][at_ind]), - "amps": - np.array(pack_table["tau"][:, 0][at_ind].data), - "amps_limits": - _ingest_limits( - pack_table["tau"][:, 0][at_ind], - pack_table["tau"][:, 1][at_ind], - ), - "amps_fixed": - _ingest_fixed(pack_table["tau"][:, 1][at_ind]), - "fwhms": - np.array(pack_table["fwhm"][:, 0][at_ind].data), - "fwhms_limits": - _ingest_limits( - pack_table["fwhm"][:, 1][at_ind], - pack_table["fwhm"][:, 2][at_ind], - ), - "fwhms_fixed": - _ingest_fixed(pack_table["fwhm"][:, 1][at_ind]), - } - - att_info = None - if len(at_ind) > 1: - raise NotImplementedError("More than one attenuation component not supported") - elif len(at_ind) == 1: - i = at_ind[0] - att_info = {"name": pack_table["name"][i], - "tau_sil": pack_table["tau"][i][0], - "tau_sil_limits": pack_table["tau"][i][1:], - "tau_sil_fixed": True if pack_table["tau"][i].mask[1] else False} - - return [bb_info, df_info, h2_info, ion_info, abs_info, att_info] diff --git a/pahfit/component_models.py b/pahfit/component_models.py deleted file mode 100644 index db719492..00000000 --- a/pahfit/component_models.py +++ /dev/null @@ -1,136 +0,0 @@ -import numpy as np -from scipy import interpolate -from astropy.modeling.physical_models import Drude1D -from astropy.modeling import Fittable1DModel -from astropy.modeling import Parameter - -__all__ = ["BlackBody1D", "ModifiedBlackBody1D", "S07_attenuation", "att_Drude1D"] - - -class BlackBody1D(Fittable1DModel): - """ - A blackbody component. - - Current astropy BlackBody1D does not play well with Lorentz1D and Gauss1D - maybe, need to check again, possibly a units issue - """ - - amplitude = Parameter() - temperature = Parameter() - - @staticmethod - def evaluate(x, amplitude, temperature): - """ - """ - return ( - amplitude - * 3.97289e13 - / x ** 3 - / (np.exp(1.4387752e4 / x / temperature) - 1.0) - ) - - -class ModifiedBlackBody1D(BlackBody1D): - """ - Modified blackbody with an emissivity propoportional to nu^2 - """ - - @staticmethod - def evaluate(x, amplitude, temperature): - return BlackBody1D.evaluate(x, amplitude, temperature) * ((9.7 / x) ** 2) - - -class S07_attenuation(Fittable1DModel): - """ - Smith, Draine, et al. (2007) kvt attenuation model calculation. - Calculation is for a fully mixed geometrically model. - Uses an extinction curve based on the silicate profiles from - Kemper, Vriend, & Tielens (2004, apJ, 609, 826). - Constructed as a weighted sum of two components: silicate profiles, - and an exponent 1.7 power-law. - - Attenuation curve for a mixed case calculated from - .. math:: - - Att(x) = \\frac{1 - e^{-\\tau_{x}}}{\\tau_{x}} - - Parameters - ---------- - kvt_amp : float - amplitude - """ - - # Attenuation tau - tau_sil = Parameter(description="kvt term: amplitude", default=1.0, min=0.0, max=10) - - @staticmethod - def kvt(in_x): - """ - Output the kvt extinction curve - """ - # fmt: off - kvt_wav = np.array([8.0, 8.2, 8.4, 8.6, 8.8, 9.0, 9.2, 9.4, 9.6, 9.7, - 9.75, 9.8, 10.0, 10.2, 10.4, 10.6, 10.8, 11.0, - 11.2, 11.4, 11.6, 11.8, 12.0, 12.2, 12.4, 12.6, - 12.7]) - - kvt_int = np.array([.06, .09, .16, .275, .415, .575, .755, .895, .98, - .99, 1.0, .99, .94, .83, .745, .655, .58, .525, - .43, .35, .27, .20, .13, .09, .06, .045, .04314]) - # fmt: on - - # Extend kvt profile to shorter wavelengths - if min(in_x) < min(kvt_wav): - kvt_wav_short = in_x[in_x < min(kvt_wav)] - kvt_int_short_tmp = min(kvt_int) * np.exp(2.03 * (kvt_wav_short - min(kvt_wav))) - # Since kvt_int_shoft_tmp does not reach min(kvt_int), - # we scale it to stitch it. - kvt_int_short = kvt_int_short_tmp * (kvt_int[0] / max(kvt_int_short_tmp)) - - spline_x = np.concatenate([kvt_wav_short, kvt_wav]) - spline_y = np.concatenate([kvt_int_short, kvt_int]) - else: - spline_x = kvt_wav - spline_y = kvt_int - - intfunc = interpolate.interp1d(spline_x, spline_y) - in_x_spline = in_x[in_x < max(kvt_wav)] - new_spline_y = intfunc(in_x_spline) - - nf = Drude1D(amplitude=0.4, x_0=18.0, fwhm=0.247 * 18.0) - in_x_drude = in_x[in_x >= max(kvt_wav)] - - ext = np.concatenate([new_spline_y, nf(in_x_drude)]) - - # Extend to ~2 um - # assuming beta is 0.1 - beta = 0.1 - y = (1.0 - beta) * ext + beta * (9.7 / in_x) ** 1.7 - - return y - - def evaluate(self, in_x, tau_si): - if tau_si == 0.0: - return np.full((len(in_x)), 1.0) - else: - tau_x = tau_si * self.kvt(in_x) - return (1.0 - np.exp(-1.0 * tau_x)) / tau_x - - -class att_Drude1D(Fittable1DModel): - """ - Attenuation components that can be parameterized by Drude profiles. - """ - - tau = Parameter() - x_0 = Parameter() - fwhm = Parameter() - - @staticmethod - def evaluate(x, tau, x_0, fwhm): - if tau == 0.0: - return np.full((len(x)), 0.0) - else: - profile = Drude1D(amplitude=1.0, fwhm=fwhm, x_0=x_0) - tau_x = tau * profile(x) - return (1.0 - np.exp(-1.0 * tau_x)) / tau_x diff --git a/pahfit/feature_strengths.py b/pahfit/feature_strengths.py index f1c67228..40eba460 100644 --- a/pahfit/feature_strengths.py +++ b/pahfit/feature_strengths.py @@ -1,5 +1,5 @@ from astropy.modeling.functional_models import Gaussian1D -from pahfit.component_models import BlackBody1D +from pahfit.fitters.ap_components import BlackBody1D import numpy as np diff --git a/pahfit/features/features.py b/pahfit/features/features.py index da5550e5..b563bbef 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -21,16 +21,33 @@ import numpy as np from astropy.table import vstack, Table, TableAttribute from astropy.io.misc.yaml import yaml -import astropy.units as u -from pkg_resources import resource_filename +from importlib import resources from pahfit.errors import PAHFITFeatureError -from pahfit.features.features_format import BoundedMaskedColumn, BoundedParTableFormatter +from pahfit.features.features_format import BoundedParTableFormatter +import pahfit.units + +# Feature kinds and associated parameters +KIND_PARAMS = {'starlight': {'temperature', 'tau'}, + 'dust_continuum': {'temperature', 'tau'}, + 'line': {'wavelength', 'power'}, # 'fwhm', Instrument Pack detail! + 'dust_feature': {'wavelength', 'fwhm', 'power'}, + 'attenuation': {'model', 'tau', 'geometry'}, + 'absorption': {'wavelength', 'fwhm', 'tau', 'geometry'}} + +# Parameter default units: flux density/intensity/power (other units determined on fit) +# Note: power is actually amplitude for now. Change this to +# intensity_power when the power features are implemented. + +PARAM_UNITS = {'temperature': pahfit.units.temperature, + 'wavelength': pahfit.units.wavelength, + 'fwhm': pahfit.units.wavelength, + 'power': pahfit.units.intensity_power} class UniqueKeyLoader(yaml.SafeLoader): def construct_mapping(self, node, deep=False): mapping = set() - for key_node, value_node in node.value: + for key_node, _ in node.value: key = self.construct_object(key_node, deep=deep) if key in mapping: raise PAHFITFeatureError(f"Duplicate {key!r} key found in YAML.") @@ -69,8 +86,8 @@ def value_bounds(val, bounds): Returns: ------- - The value, if unbounded, or a 3 element tuple (value, min, max). - Any missing bound is replaced with the numpy `masked' value. + A 3 element tuple (value, min, max). + Any missing bound is replaced with the numpy.nan value. Raises: ------- @@ -82,7 +99,7 @@ def value_bounds(val, bounds): if val is None: val = np.ma.masked if not bounds: - return (val,) + 2 * (np.ma.masked,) # Fixed + return (val,) + 2 * (np.nan,) # (val,nan,nan) indicates fixed ret = [val] for i, b in enumerate(bounds): if isinstance(b, str): @@ -102,27 +119,25 @@ def value_bounds(val, bounds): class Features(Table): - """A class for holding PAHFIT features and their associated - parameter information. Note that each parameter has an associated - `kind', and that each kind has an associated set of allowable - parameters (see _kind_params, below). + """A class for holding a table of PAHFIT features and associated + parameter information. + + Note that each parameter has an associated `kind', and that each + kind has an associated set of allowable parameters (see + `KIND_PARAMS`). + + See Also + -------- + `~astropy.table.Table`: The parent table class. """ TableFormatter = BoundedParTableFormatter - MaskedColumn = BoundedMaskedColumn param_covar = TableAttribute(default=[]) - _kind_params = {'starlight': {'temperature', 'tau'}, - 'dust_continuum': {'temperature', 'tau'}, - 'line': {'wavelength', 'power'}, # 'fwhm', Instrument Pack detail! - 'dust_feature': {'wavelength', 'fwhm', 'power'}, - 'attenuation': {'model', 'tau', 'geometry'}, - 'absorption': {'wavelength', 'fwhm', 'tau', 'geometry'}} - - _units = {'temperature': u.K, 'wavelength': u.um, 'fwhm': u.um} _group_attrs = set(('bounds', 'features', 'kind')) # group-level attributes _param_attrs = set(('value', 'bounds')) # Each parameter can have these attributes - _no_bounds = set(('name', 'group', 'geometry', 'model')) # String attributes (no bounds) + _no_bounds = set(('name', 'group', 'kind', 'geometry', 'model')) # str attributes (no bounds) + _bounds_dtype = np.dtype([("val", float), ("min", float), ("max", float)]) # bounded param type @classmethod def read(cls, file, *args, **kwargs): @@ -158,8 +173,7 @@ def _read_scipack(cls, file): feat_tables = dict() if not os.path.isfile(file): - pack_path = resource_filename("pahfit", "packs/science") - file = os.path.join(pack_path, file) + file = resources.files("pahfit") / "packs/science" / file try: with open(file) as fd: scipack = yaml.load(fd, Loader=UniqueKeyLoader) @@ -179,7 +193,7 @@ def _read_scipack(cls, file): raise PAHFITFeatureError(f"No kind found for {name}\n\t{file}") try: - valid_params = cls._kind_params[kind] + valid_params = KIND_PARAMS[kind] except KeyError: raise PAHFITFeatureError(f"Unknown kind {kind} for {name}\n\t{file}") unknown_params = [x for x in keys @@ -264,7 +278,7 @@ def _add_feature(cls, kind: str, t: dict, name: str, *, t[kind][name]['group'] = group t[kind][name]['kind'] = kind for (param, val) in pars.items(): - if param not in cls._kind_params[kind]: + if param not in KIND_PARAMS[kind]: continue if isinstance(val, dict): # A param attribute dictionary unknown_attrs = [x for x in val.keys() if x not in cls._param_attrs] @@ -307,28 +321,36 @@ def _construct_table(cls, inp: dict): """ tables = [] for (kind, features) in inp.items(): - kind_params = cls._kind_params[kind] # All params for this kind + if kind == "_ratios": + continue + kp = KIND_PARAMS[kind] # All params for this kind rows = [] for (name, params) in features.items(): - for missing in kind_params - params.keys(): + for missing in kp - params.keys(): if missing in cls._no_bounds: params[missing] = 0.0 else: params[missing] = value_bounds(0.0, bounds=(0.0, None)) rows.append(dict(name=name, **params)) - table_columns = rows[0].keys() - t = cls(rows, names=table_columns) - for p in cls._kind_params[kind]: - if p not in cls._no_bounds: - t[p].info.format = "0.4g" # Nice format (customized by Formatter) + param_names = rows[0].keys() + dtypes = [str if x in cls._no_bounds else cls._bounds_dtype for x in param_names] + t = cls(rows, names=param_names, dtype=dtypes) tables.append(t) tables = vstack(tables) for cn, col in tables.columns.items(): - if cn in cls._units: - col.unit = cls._units[cn] - tables.add_index('name') + if cn in PARAM_UNITS: + col.unit = PARAM_UNITS[cn] + cls._index_table(tables) + + if '_ratios' in inp: + tables.meta['_ratios'] = inp['_ratios'] return tables + @staticmethod + def _index_table(tbl): + for indx in ('name', 'group'): + tbl.add_index(indx) + def mask_feature(self, name, mask_value=True): """Mask all the parameters of a feature. @@ -345,15 +367,19 @@ def mask_feature(self, name, mask_value=True): """ row = self.loc[name] - relevant_params = self._kind_params[row['kind']] + relevant_params = KIND_PARAMS[row['kind']] for col_name in relevant_params: if col_name in self._no_bounds: # these are all strings, so can't mask pass else: # mask only the value, not the bounds - row[col_name].mask[0] = mask_value + row[col_name].mask['val'] = mask_value def unmask_feature(self, name): """Remove the mask for all parameters of a feature.""" self.mask_feature(name, mask_value=False) + + def _base_repr_(self, *args, **kwargs): + """Omit dtype on self-print.""" + return super()._base_repr_(*args, ** kwargs | dict(show_dtype=False)) diff --git a/pahfit/features/features_format.py b/pahfit/features/features_format.py index 78369727..2871e936 100644 --- a/pahfit/features/features_format.py +++ b/pahfit/features/features_format.py @@ -1,53 +1,45 @@ -import numpy.ma as ma -from astropy.table import MaskedColumn +import numpy as np from astropy.table.pprint import TableFormatter # * Special table formatting for bounded (val, min, max) values -def fmt_func(fmt): - def _fmt(v): - if ma.is_masked(v[0]): - return " " - if ma.is_masked(v[1]): - return f"{v[0]:{fmt}} (Fixed)" - return f"{v[0]:{fmt}} ({v[1]:{fmt}}, {v[2]:{fmt}})" - +def fmt_func(fmt: str): + """Format bounded variables specially.""" + if fmt.startswith('%'): + fmt = fmt[1:] + + def _fmt(x): + ret = f"{x['val']:{fmt}}" + if np.isnan(x['min']) and np.isnan(x['max']): + return ret + " (fixed)" + else: + mn = ("-∞" if np.isnan(x['min']) or x['min'] == -np.inf + else f"{x['min']:{fmt}}") + mx = ("∞" if np.isnan(x['max']) or x['max'] == np.inf + else f"{x['max']:{fmt}}") + return f"{ret} ({mn}, {mx})" return _fmt -class BoundedMaskedColumn(MaskedColumn): - """Masked column which can be toggled to group rows into one item - for formatting. To be set as Table's `MaskedColumn'. - """ - - _omit_shape = False - - @property - def shape(self): - sh = super().shape - return sh[0:-1] if self._omit_shape and len(sh) > 1 else sh - - def is_fixed(self): - return ma.getmask(self)[:, 1:].all(1) - - class BoundedParTableFormatter(TableFormatter): """Format bounded parameters. Bounded parameters are 3-field structured arrays, with fields - 'var', 'min', and 'max'. To be set as Table's `TableFormatter'. + 'val', 'min', and 'max'. To be set as Table's `TableFormatter'. """ - def _pformat_table(self, table, *args, **kwargs): bpcols = [] + tlfmt = table.meta.get('pahfit_format') try: - colsh = [(col, col.shape) for col in table.columns.values()] - BoundedMaskedColumn._omit_shape = True - for col, sh in colsh: - if len(sh) == 2 and sh[1] == 3: + for col in table.columns.values(): + if len(col.dtype) == 3: # bounded! bpcols.append((col, col.info.format)) - col.info.format = fmt_func(col.info.format or "g") + fmt = col.meta.get('pahfit_format') or tlfmt or "g" + col.info.format = fmt_func(fmt) return super()._pformat_table(table, *args, **kwargs) finally: - BoundedMaskedColumn._omit_shape = False for col, fmt in bpcols: col.info.format = fmt + + def _name_and_structure(self, name, *args): + "Simplified column name: no val, min, max needed." + return name diff --git a/pahfit/features/util.py b/pahfit/features/util.py index 68513aa0..d9493f94 100644 --- a/pahfit/features/util.py +++ b/pahfit/features/util.py @@ -1,29 +1,28 @@ """pahfit.util General pahfit.features utility functions.""" import numpy as np -import numpy.ma as ma def bounded_is_missing(val): """Return a mask array indicating which of the bounded values are missing. A missing bounded value has a masked value.""" - return ma.getmask(val)[..., 0] + return getattr(val['val'], 'mask', None) or np.zeros_like(val['val'], dtype=bool) def bounded_is_fixed(val): """Return a mask array indicating which of the bounded values are fixed. A fixed bounded value has masked bounds.""" - return ma.getmask(val)[..., -2:].all(-1) + return np.isnan(val['min']) & np.isnan(val['max']) def bounded_min(val): """Return the minimum of each bounded value passed. Either the lower bound, or, if no such bound is set, the value itself.""" - lower = val[..., 1] - return np.where(lower, lower, val[..., 0]) + lower = val['min'] + return np.where(lower, lower, val['val']) def bounded_max(val): """Return the maximum of each bounded value passed. Either the upper bound, or, if no such bound is set, the value itself.""" - upper = val[..., 2] - return np.where(upper, upper, val[..., 0]) + upper = val['max'] + return np.where(upper, upper, val['val']) diff --git a/pahfit/fitters/__init__.py b/pahfit/fitters/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/pahfit/fitters/ap_components.py b/pahfit/fitters/ap_components.py new file mode 100644 index 00000000..921e23b3 --- /dev/null +++ b/pahfit/fitters/ap_components.py @@ -0,0 +1,283 @@ +import numpy as np +from scipy import interpolate +from astropy.modeling.physical_models import Drude1D +from astropy.modeling import Fittable1DModel +from astropy.modeling import Parameter +from astropy import constants +from pahfit import units + +__all__ = ["BlackBody1D", "ModifiedBlackBody1D", "S07_attenuation", "att_Drude1D"] + + +class BlackBody1D(Fittable1DModel): + """ + A blackbody component. + + Current astropy BlackBody1D does not play well with Lorentz1D and Gauss1D + maybe, need to check again, possibly a units issue + """ + + amplitude = Parameter() + temperature = Parameter() + + @staticmethod + def evaluate(x, amplitude, temperature): + """ """ + return ( + amplitude + * 3.97289e13 + / x**3 + / (np.exp(1.4387752e4 / x / temperature) - 1.0) + ) + + +class ModifiedBlackBody1D(BlackBody1D): + """ + Modified blackbody with an emissivity propoportional to nu^2 + """ + + @staticmethod + def evaluate(x, amplitude, temperature): + return BlackBody1D.evaluate(x, amplitude, temperature) * ((9.7 / x) ** 2) + + +class S07_attenuation(Fittable1DModel): + """ + Smith, Draine, et al. (2007) kvt attenuation model calculation. + Calculation is for a fully mixed geometrically model. + Uses an extinction curve based on the silicate profiles from + Kemper, Vriend, & Tielens (2004, apJ, 609, 826). + Constructed as a weighted sum of two components: silicate profiles, + and an exponent 1.7 power-law. + + Attenuation curve for a mixed case calculated from + .. math:: + + Att(x) = \\frac{1 - e^{-\\tau_{x}}}{\\tau_{x}} + + Parameters + ---------- + kvt_amp : float + amplitude + """ + + # Attenuation tau + tau_sil = Parameter(description="kvt term: amplitude", default=1.0, min=0.0, max=10) + + @staticmethod + def kvt(in_x): + """ + Output the kvt extinction curve + """ + # fmt: off + kvt_wav = np.array([8.0, 8.2, 8.4, 8.6, 8.8, 9.0, 9.2, 9.4, 9.6, 9.7, + 9.75, 9.8, 10.0, 10.2, 10.4, 10.6, 10.8, 11.0, + 11.2, 11.4, 11.6, 11.8, 12.0, 12.2, 12.4, 12.6, + 12.7]) + + kvt_int = np.array([.06, .09, .16, .275, .415, .575, .755, .895, .98, + .99, 1.0, .99, .94, .83, .745, .655, .58, .525, + .43, .35, .27, .20, .13, .09, .06, .045, .04314]) + # fmt: on + + # Extend kvt profile to shorter wavelengths + if min(in_x) < min(kvt_wav): + kvt_wav_short = in_x[in_x < min(kvt_wav)] + kvt_int_short_tmp = min(kvt_int) * np.exp( + 2.03 * (kvt_wav_short - min(kvt_wav)) + ) + # Since kvt_int_shoft_tmp does not reach min(kvt_int), + # we scale it to stitch it. + kvt_int_short = kvt_int_short_tmp * (kvt_int[0] / max(kvt_int_short_tmp)) + + spline_x = np.concatenate([kvt_wav_short, kvt_wav]) + spline_y = np.concatenate([kvt_int_short, kvt_int]) + else: + spline_x = kvt_wav + spline_y = kvt_int + + intfunc = interpolate.interp1d(spline_x, spline_y) + in_x_spline = in_x[in_x < max(kvt_wav)] + new_spline_y = intfunc(in_x_spline) + + nf = Drude1D(amplitude=0.4, x_0=18.0, fwhm=0.247 * 18.0) + in_x_drude = in_x[in_x >= max(kvt_wav)] + + ext = np.concatenate([new_spline_y, nf(in_x_drude)]) + + # Extend to ~2 um + # assuming beta is 0.1 + beta = 0.1 + y = (1.0 - beta) * ext + beta * (9.7 / in_x) ** 1.7 + + return y + + def evaluate(self, in_x, tau_si): + if tau_si == 0.0: + return np.full((len(in_x)), 1.0) + else: + tau_x = tau_si * self.kvt(in_x) + return (1.0 - np.exp(-1.0 * tau_x)) / tau_x + + +class att_Drude1D(Fittable1DModel): + """ + Attenuation components that can be parameterized by Drude profiles. + """ + + tau = Parameter() + x_0 = Parameter() + fwhm = Parameter() + + @staticmethod + def evaluate(x, tau, x_0, fwhm): + if tau == 0.0: + return np.full((len(x)), 0.0) + else: + profile = Drude1D(amplitude=1.0, fwhm=fwhm, x_0=x_0) + tau_x = tau * profile(x) + return (1.0 - np.exp(-1.0 * tau_x)) / tau_x + + +class PowerDrude1D(Fittable1DModel): + """ + Drude profile with amplitude determined by power. + + Special attention is needed for the units. An implementation for + this is 'unitful' because 'power' is defined as integral over + frequency, while the profile is formulated as a function of + wavelength. If we assume the output has unit(flux), then without + additional conversions the power unit will be unit(flux) * frequency + = unit(flux) * unit(c) / unit(lam). + + Example: if x_0 and fwhm are in micron, and the flux is in MJy / sr, + then the unit of the fitted power will be MJy sr-1 * unit(c) / + micron, which differs by a constant factor from MJy sr-1 Hz, + depending on the chosen unit of c. + + For efficiency and to prevent ambiguity, we assume that all units + are the internal pahfit units defined in pahfit.units, and + precalculate a conversion factor. + + TODO: We need to check if the flux is 'intensity' or 'flux_density', + and assume the power parameter has 'intensity_power' or + 'flux_density_power' units respectively. For now, only intensity is + supported. + + """ + + power = Parameter(min=0.0) + x_0 = Parameter(min=0.0) + fwhm = Parameter(default=1, min=0.0) + + # constant factors in the equation to convert power to amplitude of + # the profile. + intensity_amplitude_factor = ( + (2 * units.intensity_power * units.wavelength / (constants.c * np.pi)) + .to(units.intensity) + .value + ) + + def evaluate(self, x, power, x_0, fwhm): + """ + Smith, et al. (2007) dust features model. Calculation is for a + Drude profile (equation in section 4.1.4). + + The intensity profile as a function of wavelength is + + Inu(lambda) = (b * g**2) / ((lambda / x0 - x0 / lambda)**2 + g**2) + + With + b = amplitude (has same unit as flux) + g = fwhm / x0 + x0 = central wavelength + + The integrated power (Fnu dnu) of the profile is + + P = (pi * c / 2) * (b * g / x0) + + Which can be solved for the amplitude b. + + b = (P * 2 * x0) / (pi * c * g) + + According to the above equations, without additional + conversions, the resulting amplitude unit will be unit(P) * + Hz-1. The factor x0 / c = 1 / nu0 (Hz-1) will result in very + small values for for Inu(lambda), or very large values for P. To + avoid numerical problems, we apply a conversion that ensures + Inu(lambda) and P are in internal units. + + Parameters + ---------- + power : float + fwhm : float + central intensity (x_0) : float + + """ + # The equation and unit conversion for the amplitude: + # b = (2 * P * x_0 / (pi * c * g)).to(output_unit).value + + # Use predetermined factor that deals with units and constants. + # factor = (2 * unit(power) * unit(wavelength) / (pi * c)).to(unit(intensity)) + + g = fwhm / x_0 + b = power * x_0 / g * self.intensity_amplitude_factor + return b * g**2 / ((x / x_0 - x_0 / x) ** 2 + g**2) + + +class PowerGaussian1D(Fittable1DModel): + """ + Gaussian profile with amplitude derived from power. + + Implementation analogous to PowerDrude1D. + + The amplitude of a gaussian profile given its power P, is P / + (stddev sqrt(2 pi)). Since stddev is given in wavelength units, this + equation gives the peak density per wavelength interval, Alambda. + The profile Flambda is then + + Flambda(lambda) = Alambda * G(lambda; mean, stddev) + + where G is a gaussian profile with amplitude 1. Converting this to + Fnu units yields + + Fnu(lambda) = lambda**2 / c * Flambda(lambda) + + Approximating the lambda**2 factor as a constant (the central + wavelength = 'mean'), then yields + + Fnu(lambda) = mean**2 / c * Alambda * G(lambda; mean, stddev) + + In other words, for narrow lines, the per-frequency profile is + approximately a gaussian with amplitude + + Anu = P * mean**2 / (c * stddev sqrt(2 pi)). + + So the constant factor we can set is + (unit(power) * unit(wavelength)**2 / (c * unit(wavelength) * sqrt(2 pi))).to(intensity) + + """ + + power = Parameter(min=0.0) + mean = Parameter() + stddev = Parameter(default=1, min=0.0) + + intensity_amplitude_factor = ( + ( + units.intensity_power + * (units.wavelength) ** 2 + / (constants.c * units.wavelength * np.sqrt(2 * np.pi)) + ) + .to(units.intensity) + .value + ) + + def evaluate(self, x, power, mean, stddev): + """ + Evaluate F_nu(lambda) given the power. + + See class description for equations and unit notes.""" + + # amplitude in intensity units + Anu = power * mean**2 / stddev * self.intensity_amplitude_factor + return Anu * np.exp(-0.5 * np.square((x - mean) / stddev)) diff --git a/pahfit/fitters/ap_fitter.py b/pahfit/fitters/ap_fitter.py new file mode 100644 index 00000000..af8be31a --- /dev/null +++ b/pahfit/fitters/ap_fitter.py @@ -0,0 +1,406 @@ +from pahfit.fitters.fitter import Fitter +from pahfit.errors import PAHFITModelError +from pahfit.fitters.ap_components import ( + BlackBody1D, + ModifiedBlackBody1D, + S07_attenuation, + att_Drude1D, + PowerDrude1D, + PowerGaussian1D, +) +from astropy.modeling.fitting import LevMarLSQFitter +import numpy as np + + +class APFitter(Fitter): + """Astropy fitting implementation using Fitter API. + + Fitter API is still subject to change. This draft was written with + what is needed, and the Fitter class was set up based on this draft. + + APFitter implements the Fitter methods using actions that involve an + astropy-based CompoundModel. Some more implementation-specific + details for these tasks are given in the documentation of the + functions. + + Multi-segment fitting was not implemented for this subclass, but + here is a sketch for how it could work: + 1. During set up, an extra argument is passed, indicating the + segment number. Each segment gets a different list of components + to use. + 2. During finalize the joint model / something else for joint + fitting is set up. Alternatively, all the separate models are + constructed, and then a set of "unique" parameters is gathered. + An appropriate name would be the "parameter union" (as it is like + taking the union of the parameters sets for the individual + segments). + 3. During the fitting, the objective function (chi2) is the sum of + the objective functions for the individual models. The arguments + of this objective function, are the unique parameters that were + derived above. At every fitting step, these parameters will + change. Inside the objective function, these changes should be + propagated to the individual models, after which they can be + evaluated again. All the fitting algorithm needs, is the + parameter space, and the objective function value. + + Alternative idea, based on concatenated segments (not spectrally + merged! Just a raw data concatenation). + - concatenate x + - concatenate y + - bookkeep the boundaries between the concatenations + - when evaluating, fill the concatenated model spectrum using a + model that depends on the index (according to the boundaries that + were set up) + - Use the 'tied' option for parameters that should be equivalent + (e.g. starlight in segment 1 should have same temperature as + starlight in segment 2). + + """ + + def __init__(self): + """Construct a new fitter. + + After construction, use the add_feature_() functions to start + setting up a model, then call finalize(). + + """ + self.additive_components = [] + self.multiplicative_components = [] + self.feature_types = {} + self.model = None + self.message = None + + def finalize(self): + """Sum the registered components into one CompoundModel. + + To be called after a series of "add_feature_()" calls, and before + using the other functionality (fitting and evaluating). + + """ + if len(self.additive_components) > 1: + self.model = sum(self.additive_components[1:], self.additive_components[0]) + elif len(self.additive_components) == 1: + self.model = self.additive_components[0] + else: + raise PAHFITModelError("No components were set up for APFitter!") + + for c in self.multiplicative_components: + self.model *= c + + def _add_component(self, astropy_model_class, multiplicative=False, **kwargs): + """Generically add any feature as an astropy model component. + + To be finalized with finalize() + + Parameters + ---------- + astropy_model_class : class + Class symbol, e.g. Blackbody1D. + + multiplicative : bool + During finalize_model(), all components with multiplicative + == False will be added together, and then the model will be + multiplied with the components with multiplicative == True. + + kwargs : dict + Arguments for the astropy model constructor, including a + unique value for "name". Should be generated with + self._astropy_model_kwargs; the add_feature_() functions show how + to do this for each type of feature. + + """ + if multiplicative: + self.multiplicative_components.append(astropy_model_class(**kwargs)) + else: + self.additive_components.append(astropy_model_class(**kwargs)) + + def add_feature_starlight(self, name, temperature, tau): + """Register a BlackBody1D. + + Parameters + ---------- + name : str + Unique name. + + temperature : array of size 3 + Temperature for the blackbody, given as [value, lower_bound, + upper_bound]. The bounds assume the same system as the + features table, where masked == fixed, and +inf or -inf mean + unbounded. + + tau : analogous. Used as amplitude. + + """ + self.feature_types[name] = "starlight" + kwargs = self._astropy_model_kwargs( + name, ["temperature", "amplitude"], [temperature, tau] + ) + self._add_component(BlackBody1D, **kwargs) + + def add_feature_dust_continuum(self, name, temperature, tau): + """Register a ModifiedBlackBody1D. + + Analogous. Temperature and tau are used as temperature and + amplitude + + """ + self.feature_types[name] = "dust_continuum" + kwargs = self._astropy_model_kwargs( + name, ["temperature", "amplitude"], [temperature, tau] + ) + self._add_component(ModifiedBlackBody1D, **kwargs) + + def add_feature_line(self, name, power, wavelength, fwhm): + """Register a PowerGaussian1D + + Analogous. Uses an implementation of the Gaussian profile, that + directly fits the power based on the internal PAHFIT units. + + """ + self.feature_types[name] = "line" + + kwargs = self._astropy_model_kwargs( + name, + ["power", "mean", "stddev"], + [power, wavelength, fwhm / 2.355], + ) + self._add_component(PowerGaussian1D, **kwargs) + + def add_feature_dust_feature(self, name, power, wavelength, fwhm): + """Register a PowerDrude1D. + + Analogous. Uses an implementation of the Drude profile that + directly fits the power based on the internal PAHFIT units. + + """ + self.feature_types[name] = "dust_feature" + kwargs = self._astropy_model_kwargs( + name, ["power", "x_0", "fwhm"], [power, wavelength, fwhm] + ) + self._add_component(PowerDrude1D, **kwargs) + + def add_feature_attenuation(self, name, tau, model="S07", geometry="screen"): + """Register the S07 attenuation component. + + Analogous. Uses tau as tau_sil for S07_attenuation. Is + multiplicative. + + model : string to select attenuation shape. Only 'S07' is supported for now. + + geometry : string to select different geometries. Only 'screen' + is available for now. + + """ + self.feature_types[name] = "attenuation" + kwargs = self._astropy_model_kwargs(name, ["tau_sil"], [tau]) + self._add_component(S07_attenuation, multiplicative=True, **kwargs) + + def add_feature_absorption(self, name, tau, wavelength, fwhm, geometry="screen"): + """Register an absorbing Drude1D component. + + Analogous. Is multiplicative. + + """ + self.feature_types[name] = "absorption" + kwargs = self._astropy_model_kwargs( + name, ["tau", "x_0", "fwhm"], [tau, wavelength, fwhm] + ) + self._add_component(att_Drude1D, multiplicative=True, **kwargs) + + def evaluate(self, lam): + """Evaluate internal astropy model with its current parameters. + + Parameters + ---------- + lam : array + Rest frame wavelengths in micron + + Returns + ------- + flux : array + Rest frame flux in internal units + """ + return self.model(lam) + + def fit(self, lam, flux, unc, maxiter=10000): + """Fit the internal model using the astropy fitter. + + The fitter class is unit agnostic, and deal with the numbers the + Model tells it to deal with. Internal renormalizations could be + good to consider, as long as any values are converted back to + the original system before returning them. In practice, the + input spectrum is expected to be in internal units, and orrected + for redshift (models operate in the rest frame). + + After the fit, the results can be retrieved via get_result(). + + Retrieval of uncertainties and fit details is yet to be + implemented. + + CAVEAT: flux unit (flux) is still ambiguous, since it can be + flux density or intensity, according to the options defined in + pahfit.units. After the fit, the return units of "power" in + get_results depend on the given spectrum (they will be flux unit + times wavelength unit). + + Parameters + ---------- + lam : array + Rest frame wavelengths in micron + + flux : array + Rest frame flux in internal units. + + unc : array + Uncertainty on rest frame flux. Same units as flux. + + """ + # clean, because astropy does not like nan + w = 1 / unc + + # make sure there are no zero uncertainties either + mask = np.isfinite(lam) & np.isfinite(flux) & np.isfinite(w) + + self.fit_info = [] + + fit = LevMarLSQFitter(calc_uncertainties=True) + astropy_result = fit( + self.model, + lam[mask], + flux[mask], + weights=w[mask], + maxiter=maxiter, + epsilon=1e-10, + acc=1e-10, + ) + self.fit_info = fit.fit_info + self.model = astropy_result + self.message = fit.fit_info["message"] + + def get_result(self, component_name): + """Retrieve results from astropy model component. + + Every relevant value inside the astropy model, need to be + written to the right position in the features table. For some + cases (amplitude/power, fwhm/stddev), conversions are necessary + (generally the inverse conversions of what was done in the + register function). + + Parameters + ---------- + component_name : str + One of the names provided to any of the add_feature_*() calls + made during setup. + + Returns + ------- + dict with Parameters according to the PAHFIT definitions. + + e.g., for a feature with amplitude, stddev, and mean parameters: + {'power': converted from amplitude, 'fwhm': converted from + stddev, 'mean': wavelength} + + """ + if self.model is None: + raise PAHFITModelError("Model not finalized yet.") + + if hasattr(self.model, "submodel_names"): + component = self.model[component_name] + else: + # deals with edge case with single component, so is not + # CompoundModel but normal single-component model. + component = self.model + + c_type = self.feature_types[component_name] + if c_type == "starlight" or c_type == "dust_continuum": + return { + "temperature": component.temperature.value, + "tau": component.amplitude.value, + } + elif c_type == "line": + return { + "power": component.power.value, + "wavelength": component.mean.value, + "fwhm": component.stddev.value * 2.355, + } + elif c_type == "dust_feature": + return { + "power": component.power.value, + "wavelength": component.x_0.value, + "fwhm": component.fwhm.value, + } + elif c_type == "attenuation": + return {"tau": component.tau_sil.value} + elif c_type == "absorption": + return { + "tau": component.tau.value, + "wavelength": component.x_0.value, + "fwhm": component.fwhm.value, + } + else: + raise PAHFITModelError(f"Unsupported component type: {c_type}") + + @staticmethod + def _astropy_model_kwargs(component_name, param_names, param_values): + """Create kwargs for an astropy model constructor. + + This is a utility that deduplicates the logic for going from + (value, min, max) tuples, to astropy model constructor keyword + arguments as in the following example: + + AstropyModelClass(name="feature name", + param1=value1, + param2=value2, + bounds={param1: (min,max), param2:(min,max)}, + fixed={param1: True, param2: False}) + + The returned arguments are in a dict that looks as follows, and + can be passed to the appropriate astropy model constructor using + **kwargs. + + {"name": "feature name" + param_name: double, ..., + "bounds": {param_name: array of size 2, ...}, + "fixed": {param_name: True or False, ...}} + + Parameters: + ----------- + component_name : str + Unique name for the component. Will later be used for + indexing the components in the astropy model. + + param_names : list of str + Names of the parameters for the astropy model, e.g. + ["dust_feature1", "dust_feature2"] + + param_values : list of (array of size 3 OR scalar) + One for each param name, each in the format of [value, + min_bound, max_bound] for variable parameters, or a scalar + (single float) for fixed parameters. + + Returns + ------- + dict : kwargs to be used in an astropy model constructor + + """ + # basic format of constructor parameters of astropy model + kwargs = {"name": component_name, "bounds": {}, "fixed": {}} + + for param_name, tuple_or_scalar in zip(param_names, param_values): + if np.isscalar(tuple_or_scalar): + is_fixed = True + value, lo_bound, up_bound = tuple_or_scalar, 0, 0 + else: + is_fixed = False + value, lo_bound, up_bound = tuple_or_scalar + + # For the limits, use 0 if fixed, the raw values if + # variable, but None if infinite (this is the convention for + # astropy modeling) + kwargs[param_name] = value + kwargs["fixed"][param_name] = is_fixed + kwargs["bounds"][param_name] = [ + None if np.isinf(x) else x for x in (lo_bound, up_bound) + ] + + return kwargs diff --git a/pahfit/fitters/fitter.py b/pahfit/fitters/fitter.py new file mode 100644 index 00000000..a75f2bde --- /dev/null +++ b/pahfit/fitters/fitter.py @@ -0,0 +1,189 @@ +from abc import ABC, abstractmethod + + +class Fitter(ABC): + """Abstract base class for interal Fitter API. + + All shared methods should have the same arguments, enforced by this + abstract class. Any API-specific options preferably go into the + constructor of the subclass, although some general-purpose + dictionaries could also be used if absolutely necessary. + + The main functionalities of a Fitter subclass: + 1. Convert the numbers that are in the Features table to a fittable + model configuration for a certain framework. The details of the + fitting framework are hidden behind the respective subclass. + 2. Fit the model to the spectrum without any additional assumptions. + The Fitter will fit the given data using the given model without + thinking about redshift, units, instrumental effects. + 3. Retrieve the fitted quantities, which are the values that were + passed during step 1. When fit result uncertainties are + implemented, they will also need to be retrieved through this + API. + 4. Access to the evaluation of the underlying model (again with no + assumptions like in step 2.). + + A few notes on how the above is achieved: + + For the model setup, there is one function per type of component + supported by PAHFIT, and the arguments of these functions will ask + for certain standard PAHFIT quantities (in practice, these are the + values stored in the Features table). The abstract Fitter class + ensure that the function signatures are the same between different + Fitter implementations, so that only a single logic has to be + implemented to up the Fitter (in practice, this is a loop over the + Features table implemented in Model). + + During the Fitter setup, the initial values, bounds, and "fixed" + flags are passed using one function call for each component, e.g. + add_feature_line()). Once all components have been added, the + finalize() function should be called; some subclasses (e.g. + APFitter) need to consolidate the registered components to prepare + the model that they manage for fitting. After this, fit() can be + called to apply the model and the astropy fitter to the data. The + results will then be retrievable for one component at a time, by + passing the component name to get_result(). + + """ + + @abstractmethod + def finalize(self): + """Process the registered features and prepare for fitting. + + The register functions below allow adding individual features. + The exact implementation of how features are added, and + finalized in to a single fittable model, depend on the + underlying implementation. + + """ + pass + + @abstractmethod + def add_feature_starlight(self, name, temperature, tau): + """Register a starlight feature. + + The exact representation depends on the implementation, but the + meaning of the parameters should be equivalent. + + Parameters + ---------- + name : str + Unique name. Will be used to allow retrieval of the results + after the fitting. + + temperature : array of size 3 or scalar + Blackbody temperature. Given as [value, lower_bound, + upper_bound] if the parameter should be variable (and + fitted). Given as scalar if parameter should be fixed. + + tau : array of size 3 + Analogously, used as power. + + """ + pass + + @abstractmethod + def add_feature_dust_continuum(self, name, temperature, tau): + """Register a dust continuum feature.""" + pass + + @abstractmethod + def add_feature_line(self, name, power, wavelength, fwhm): + """Register an emission line feature. + + Typically a Gaussian profile. + + """ + pass + + @abstractmethod + def add_feature_dust_feature(self, name, power, wavelength, fwhm): + """Register a dust feature. + + Typically a Drude profile. + + """ + pass + + @abstractmethod + def add_feature_attenuation(self, name, tau, model="S07", geometry="screen"): + """Register the S07 attenuation component. + + Other types of attenuation might be possible in the future. Is + multiplicative. + + """ + pass + + @abstractmethod + def add_feature_absorption(self, name, tau, wavelength, fwhm, geometry="screen"): + """Register an absorption feature. + + Modeled by a Drude profile. Is multiplicative. + + """ + pass + + @abstractmethod + def evaluate(self, lam): + """Evaluate the fitting function at the given wavelengths. + + Parameters + ---------- + lam : array + Rest frame wavelengths in micron + + Returns + ------- + flux : array + Rest frame flux in internal units + + """ + pass + + @abstractmethod + def fit(self, lam, flux, unc, maxiter=1000): + """Perform the fit using the framework of the subclass. + + Fitter is unit agnostic, and deals with the numbers the Model + tells it to deal with. In practice, the input spectrum is + expected to be in internal units, and corrected for redshift + (models operate in the rest frame). + + After the fit, the results can be retrieved via get_result(). + + Parameters + ---------- + lam : array + Rest frame wavelengths in micron + + flux : array + Rest frame flux in internal units. + + unc : array + Uncertainty on rest frame flux. Same units as flux. + + """ + pass + + @abstractmethod + def get_result(self, feature_name): + """Retrieve results from underlying model after fit. + + Parameters + ---------- + component_name : str + One of the names provided to any of the add_feature_() calls + made during setup. + + Returns + ------- + dict : parameters according to the PAHFIT definitions. Keys are + the same as the function signature of the relevant register + function. Values are in the same format as Features, and can + therefore be directly filled in. + + e.g. {'name': 'line0', 'power': value, 'fwhm': value, 'wavelength': value} + + """ + pass diff --git a/pahfit/helpers.py b/pahfit/helpers.py index cc42836e..f6bcabc2 100644 --- a/pahfit/helpers.py +++ b/pahfit/helpers.py @@ -1,13 +1,13 @@ import os -import pkg_resources +from importlib import resources -from specutils import Spectrum1D +from pahfit import units -from pahfit.component_models import BlackBody1D, S07_attenuation -from astropy.modeling.physical_models import Drude1D -from astropy.modeling.functional_models import Gaussian1D +from specutils import Spectrum1D +from astropy import units as u +from astropy.nddata import StdDevUncertainty -__all__ = ["read_spectrum", "calculate_compounds"] +__all__ = ["find_packfile", "read_spectrum"] def find_packfile(packfile): @@ -32,14 +32,13 @@ def find_packfile(packfile): if os.path.isfile(packfile): packfile_found = packfile else: - pack_path = pkg_resources.resource_filename("pahfit", "packs/science/") - test_packfile = "{}/{}".format(pack_path, packfile) + test_packfile = resources.files("pahfit") / "packs/science" / packfile if os.path.isfile(test_packfile): packfile_found = test_packfile else: raise ValueError("Input packfile {} not found".format(packfile)) - return packfile_found + return str(packfile_found) def read_spectrum(specfile, format=None): @@ -63,8 +62,7 @@ def read_spectrum(specfile, format=None): """ # resolve filename if not os.path.isfile(specfile): - pack_path = pkg_resources.resource_filename("pahfit", "data/") - test_specfile = "{}/{}".format(pack_path, specfile) + test_specfile = resources.files("pahfit") / "data" / specfile if os.path.isfile(test_specfile): specfile = test_specfile else: @@ -75,7 +73,7 @@ def read_spectrum(specfile, format=None): tformat = None # process user-specified or filename extension based format if format is None: - suffix = specfile.split(".")[-1].lower() + suffix = specfile.name.split(".")[-1].lower() if suffix == "ecsv": tformat = "ECSV" elif suffix == "ipac": @@ -83,104 +81,19 @@ def read_spectrum(specfile, format=None): else: tformat = format - return Spectrum1D.read(specfile, format=tformat) - - -def calculate_compounds(obsdata, pmodel): - """ - Determine model compounds for total continuum, stellar continuum, - total dust continuum, combined dust features, - combined atomic and H2 lines, combined H2 lines, - combined atomic lines, and extinction model - - Parameters - ---------- - obsdata : dict - observed data where x = wavelength, y = SED, and unc = uncertainties - - pmodel : PAHFITBase model - model giving all the components and parameters - - Returns - ------- - compounds : dict - x = wavelength in microns; - tot_cont = total continuum; - stellar_cont = stellar continuum; - dust_cont = total dust continuum; - dust_features = combined dust features; - tot_lines = combined atomic and H2 emission lines; - h2_lines = combined H2 lines; - atomic_lines = combined atomic lines; - extinction_model = extinction model - """ - - # get wavelength array - x = obsdata["x"].value - - # calculate total dust continuum and total continuum (including stellar continuum) - # v2.0: first BlackBody1D is stellar continuum - cont_components = [] - - for cmodel in pmodel.model: - if isinstance(cmodel, BlackBody1D): - cont_components.append(cmodel) - stellar_cont_model = cont_components[0] - dust_cont_model = cont_components[1] - for cmodel in cont_components[2:]: - dust_cont_model += cmodel - totcont = dust_cont_model(x) + stellar_cont_model(x) - - # calculate total dust features - dust_features = [] - - for cmodel in pmodel.model: - if isinstance(cmodel, Drude1D): - dust_features.append(cmodel) - dust_features_model = dust_features[0] - for cmodel in dust_features[1:]: - dust_features_model += cmodel - - # calculate H2 spectrum - h2_features = [] - - for cmodel in pmodel.model: - if isinstance(cmodel, Gaussian1D): - if cmodel.name[0:2] == "H2": - h2_features.append(cmodel) - h2_features_model = h2_features[0] - for cmodel in h2_features[1:]: - h2_features_model += cmodel - - # calculate atomic line spectrum - atomic_features = [] - - for cmodel in pmodel.model: - if isinstance(cmodel, Gaussian1D): - if cmodel.name[0:2] != "H2": - atomic_features.append(cmodel) - atomic_features_model = atomic_features[0] - for cmodel in atomic_features[1:]: - atomic_features_model += cmodel - - # all atomic and H2 lines - totlines = h2_features_model(x) + atomic_features_model(x) - - # get extinction model - for cmodel in pmodel.model: - if isinstance(cmodel, S07_attenuation): - ext_model = cmodel(x) - - # save compounds in dictionary - compounds = {} - compounds["x"] = x - compounds["tot_cont"] = totcont - compounds["stellar_cont"] = stellar_cont_model(x) - compounds["dust_cont"] = dust_cont_model(x) - compounds["dust_features"] = dust_features_model(x) - compounds["tot_lines"] = totlines - compounds["h2_lines"] = h2_features_model(x) - compounds["atomic_lines"] = atomic_features_model(x) - compounds["extinction_model"] = ext_model - - return compounds + s = Spectrum1D.read(specfile, format=tformat) + + # Convert to intensity units by assuming an arbitrary solid angle + # for now. To be removed when dual unit support (intensity and flux + # density) is supported. + if s.flux.unit.is_equivalent(units.flux_density): + solid_angle = (3 * u.arcsec) ** 2 + alt_flux = (s.flux / solid_angle).to(units.intensity) + alt_unc_array = (s.uncertainty.array * s.flux.unit / solid_angle).to( + units.intensity + ) + s = Spectrum1D( + alt_flux, s.spectral_axis, uncertainty=StdDevUncertainty(alt_unc_array) + ) + + return s diff --git a/pahfit/instrument.py b/pahfit/instrument.py index 3828d9c7..c55bdb8a 100644 --- a/pahfit/instrument.py +++ b/pahfit/instrument.py @@ -11,11 +11,10 @@ import os from pathlib import Path -import glob import numpy as np from numpy.polynomial import Polynomial from astropy.io.misc import yaml -from pkg_resources import resource_filename +from importlib import resources from pahfit.errors import PAHFITPackError, PAHFITWarning from warnings import warn @@ -25,7 +24,7 @@ def read_instrument_packs(): """Read all instrument packs into the 'packs' variable.""" global packs - for pack in glob.glob(resource_filename("pahfit", "packs/instrument/*.yaml")): + for pack in (resources.files("pahfit") / "packs/instrument").glob("*.yaml"): try: with open(pack) as fd: p = yaml.load(fd) diff --git a/pahfit/model.py b/pahfit/model.py index af4382fb..70b55ae1 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -1,17 +1,19 @@ from specutils import Spectrum1D from astropy import units as u +from astropy import constants import copy -from astropy.modeling.fitting import LevMarLSQFitter +import matplotlib as mpl from matplotlib import pyplot as plt import numpy as np from scipy import interpolate, integrate -from pahfit.features.util import bounded_is_fixed +from pahfit import units +from pahfit.features.util import bounded_is_fixed, bounded_is_missing from pahfit.features import Features -from pahfit.base import PAHFITBase from pahfit import instrument from pahfit.errors import PAHFITModelError -from pahfit.component_models import BlackBody1D +from pahfit.fitters.ap_components import BlackBody1D, S07_attenuation +from pahfit.fitters.ap_fitter import APFitter class Model: @@ -180,17 +182,22 @@ def guess( """ inst, z = self._parse_instrument_and_redshift(spec, redshift) # save these as part of the model (will be written to disk too) - self.features.meta["redshift"] = inst - self.features.meta["instrument"] = z + self.features.meta["redshift"] = z + self.features.meta["instrument"] = inst # parse spectral data self.features.meta["user_unit"]["flux"] = spec.flux.unit - _, _, _, xz, yz, _ = self._convert_spec_data(spec, z) - wmin = min(xz) - wmax = max(xz) + _, _, _, lam, flux, _ = self._convert_spec_data(spec, z) + lam_min = min(lam) + lam_max = max(lam) + + # Some useful quantities for guessing + median_flux = np.median(flux) + Flambda = flux * units.intensity * (lam * units.wavelength) ** -2 * constants.c + total_power = integrate.trapezoid(Flambda, lam * units.wavelength) # simple linear interpolation function for spectrum - sp = interpolate.interp1d(xz, yz) + sp = interpolate.interp1d(lam, flux) # we will repeat this loop logic several times def loop_over_non_fixed(kind, parameter, estimate_function, force=False): @@ -204,7 +211,7 @@ def loop_over_non_fixed(kind, parameter, estimate_function, force=False): # guess starting point of bb def starlight_guess(row): bb = BlackBody1D(1, row["temperature"][0]) - w = wmin + 0.1 # the wavelength used to compare + w = lam_min + 0.1 # the wavelength used to compare if w < 5: # wavelength is short enough to not have numerical # issues. Evaluate both at w. @@ -226,71 +233,100 @@ def dust_continuum_guess(row): temp = row["temperature"][0] fmax_lam = 2898.0 / temp bb = BlackBody1D(1, temp) - if fmax_lam >= wmin and fmax_lam <= wmax: - w_ref = fmax_lam - elif fmax_lam > wmax: - w_ref = wmax + if fmax_lam >= lam_min and fmax_lam <= lam_max: + lam_ref = fmax_lam + elif fmax_lam > lam_max: + lam_ref = lam_max else: - w_ref = wmin + lam_ref = lam_min - flux_ref = np.median(yz[(xz > w_ref - 0.2) & (xz < w_ref + 0.2)]) - amp_guess = flux_ref / bb(w_ref) + flux_ref = np.median(flux[(lam > lam_ref - 0.2) & (lam < lam_ref + 0.2)]) + amp_guess = flux_ref / bb(lam_ref) return amp_guess / nbb loop_over_non_fixed("dust_continuum", "tau", dust_continuum_guess) def line_fwhm_guess(row): - w = row["wavelength"][0] - if not instrument.within_segment(w, inst): + lam_line = row["wavelength"][0] + if not instrument.within_segment(lam_line, inst): return 0 - fwhm = instrument.fwhm(inst, w, as_bounded=True)[0][0] + fwhm = instrument.fwhm(inst, lam_line, as_bounded=True)[0][0] return fwhm - def amp_guess(row, fwhm): - w = row["wavelength"][0] - if not instrument.within_segment(w, inst): + def power_guess(row, fwhm): + # local integration for the lines + lam_line = row["wavelength"][0] + if not instrument.within_segment(lam_line, inst): return 0 factor = 1.5 - wmin = w - factor * fwhm - wmax = w + factor * fwhm - xz_window = (xz > wmin) & (xz < wmax) - xpoints = xz[xz_window] - ypoints = yz[xz_window] - if np.count_nonzero(xz_window) >= 2: + lam_min = lam_line - factor * fwhm + lam_max = lam_line + factor * fwhm + lam_window = (lam > lam_min) & (lam < lam_max) + xpoints = lam[lam_window] + ypoints = flux[lam_window] + if np.count_nonzero(lam_window) >= 2: # difference between flux in window and flux around it - power_guess = integrate.trapezoid(yz[xz_window], xz[xz_window]) + Fnu_dlambda = integrate.trapezoid(flux[lam_window], lam[lam_window]) # subtract continuum estimate, but make sure we don't go negative continuum = (ypoints[0] + ypoints[-1]) / 2 * (xpoints[-1] - xpoints[0]) - if continuum < power_guess: - power_guess -= continuum + if continuum < Fnu_dlambda: + Fnu_dlambda -= continuum else: - power_guess = 0 - return power_guess / fwhm + Fnu_dlambda = 0 + + # this is an unphysical power (Fnu * dlambda), but we + # convert to Fnu dnu = Fnu dnu/dlambda dlambda = Fnu c / + # lambda **2 dlambda + Fnu_dlambda *= units.intensity * units.wavelength + Fnu_dnu = Fnu_dlambda * constants.c / (lam_line * units.wavelength) ** 2 + return Fnu_dnu.to(units.intensity_power).value - # Same logic as in the old function: just use same amp for all - # dust features. - some_flux = 0.5 * np.median(yz) - loop_over_non_fixed("dust_feature", "power", lambda row: some_flux) + def drude_power_guess(row): + # multiply total power by some fraction to guess Drude power + fwhm = row["fwhm"][0] * units.wavelength + delta_w = spec.spectral_axis[-1] - spec.spectral_axis[0] + return (total_power * fwhm / delta_w).to(units.intensity_power).value + + loop_over_non_fixed("dust_feature", "power", drude_power_guess) if integrate_line_flux: - # calc line amplitude using instrumental fwhm and integral over data + # calc line power using instrumental fwhm and integral over data loop_over_non_fixed( - "line", "power", lambda row: amp_guess(row, line_fwhm_guess(row)) + "line", "power", lambda row: power_guess(row, line_fwhm_guess(row)) ) else: - loop_over_non_fixed("line", "power", lambda row: some_flux) + loop_over_non_fixed( + "line", "power", lambda row: median_flux * line_fwhm_guess(row) + ) - # set the fwhms in the features table + # Set the fwhms in the features table. Slightly different logic, + # as the fwhm for lines are masked by default. TODO: leave FWHM + # masked for lines, and instead have a sigma_v option. Any + # requirements to guess and fit the line width, should be + # encapsulated in sigma_v (the "broadening" of the line), as + # opposed to fwhm which is the normal instrumental width. if calc_line_fwhm: - loop_over_non_fixed("line", "fwhm", line_fwhm_guess, force=True) + for row_index in np.where(self.features["kind"] == "line")[0]: + row = self.features[row_index] + if row["fwhm"] is np.ma.masked: + self.features[row_index]["fwhm"] = ( + line_fwhm_guess(row), + np.nan, + np.nan, + ) + elif not bounded_is_fixed(row["fwhm"]): + self.features[row_index]["fwhm"]["val"] = line_fwhm_guess(row) @staticmethod def _convert_spec_data(spec, z): - """ - Turn astropy quantities stored in Spectrum1D into fittable - numbers. + """Convert Spectrum1D Quantities to fittable numbers. + + The unit of the input spectrum has to be a multiple of MJy / sr, + the internal intensity unit. The output of this function + consists of simple unitless arrays (the numbers in these arrays + are assumed to be consistent with the internal units). Also corrects for redshift. @@ -298,18 +334,23 @@ def _convert_spec_data(spec, z): ------- x, y, unc: wavelength in micron, flux, uncertainty - xz, yz, uncz: wavelength in micron, flux, uncertainty + lam, flux, unc: wavelength in micron, flux, uncertainty corrected for redshift + """ - x = spec.spectral_axis.to(u.micron).value - y = spec.flux.value - unc = spec.uncertainty.array + if not spec.flux.unit.is_equivalent(units.intensity): + raise PAHFITModelError( + "For now, PAHFIT only supports intensity units, i.e. convertible to MJy / sr." + ) + flux_obs = spec.flux.to(units.intensity).value + lam_obs = spec.spectral_axis.to(u.micron).value + unc_obs = (spec.uncertainty.array * spec.flux.unit).to(units.intensity).value # transform observed wavelength to "physical" wavelength - xz = x / (1 + z) # wavelength shorter - yz = y * (1 + z) # energy higher - uncz = unc * (1 + z) # uncertainty scales with flux - return x, y, unc, xz, yz, uncz + lam = lam_obs / (1 + z) # wavelength shorter + flux = flux_obs * (1 + z) # energy higher + unc = unc_obs * (1 + z) # uncertainty scales with flux + return lam_obs, flux_obs, unc_obs, lam, flux, unc def fit( self, @@ -370,7 +411,7 @@ def fit( # parse spectral data self.features.meta["user_unit"]["flux"] = spec.flux.unit inst, z = self._parse_instrument_and_redshift(spec, redshift) - x, _, _, xz, yz, uncz = self._convert_spec_data(spec, z) + x, _, _, lam, flux, unc = self._convert_spec_data(spec, z) # save these as part of the model (will be written to disk too) self.features.meta["redshift"] = inst @@ -379,66 +420,281 @@ def fit( # check if observed spectrum is compatible with instrument model instrument.check_range([min(x), max(x)], inst) - # weigths - w = 1.0 / uncz - - # construct model and perform fit - astropy_model = self._construct_astropy_model(inst, z, use_instrument_fwhm) - fit = LevMarLSQFitter(calc_uncertainties=True) - self.astropy_result = fit( - astropy_model, - xz, - yz, - weights=w, - maxiter=maxiter, - epsilon=1e-10, - acc=1e-10, - ) - self.fit_info = fit.fit_info + self._set_up_fitter(inst, z, lam=x, use_instrument_fwhm=use_instrument_fwhm) + self.fitter.fit(lam, flux, unc, maxiter=maxiter) + + # copy the fit results to the features table + self._ingest_fit_result_to_features() + if verbose: - print(fit.fit_info["message"]) + print(self.fitter.message) + + def _ingest_fit_result_to_features(self): + """Copy the results from the Fitter to the features table - self._parse_astropy_result(self.astropy_result) + This is a utility method, executed only at the end of fit(), + where Fitter.fit() has been applied. - def info(self): - """Print out the last fit results.""" - print(self.astropy_result) + """ + # iterate over the list stored in fitter, so we only get + # components that were set up by _set_up_fitter. Having an + # ENABLED/DISABLED flag for every feature would be a nice + # alternative (and clear for the user). + + self.features.meta["fitter_message"] = self.fitter.message + + for name in self.enabled_features: + for column, value in self.fitter.get_result(name).items(): + try: + i = np.where(self.features["name"] == name)[0] + # deal with fwhm usually being masked + if not bounded_is_missing(self.features[column][i]): + self.features[column]["val"][i] = value + else: + self.features[column][i] = (value, np.nan, np.nan) + except Exception as e: + print( + f"Could not assign to name {name} in features table. Some diagnostic output below" + ) + print(f"Index i is {i}") + print("Features table:", self.features) + raise e - def plot(self, spec=None, redshift=None): + def plot( + self, + spec=None, + redshift=None, + use_instrument_fwhm=False, + label_lines=False, + scalefac_resid=2, + **errorbar_kwargs, + ): """Plot model, and optionally compare to observational data. Parameters ---------- spec : Spectrum1D - Observational data. Does not have to be the same data that - was used for guessing or fitting. + Observational data. The units should be compatible with the + data that were used for the fit, but it does not have to be + the exact same spectrum. The spectrum will be converted to + internal units before plotting. redshift : float Redshift used to shift from the physical model, to the - observed model. + observed model. If None, it will be taken from spec.redshift - If None, will be taken from spec.redshift + use_instrument_fwhm : bool + For the lines, the default is to use the fwhm values + contained in the Features table. When set to True, the fwhm + will be determined by the instrument model instead. + + label_lines : bool + Add labels with the names of the lines, at the position of + each line. + + scalefac_resid : float + Factor multiplying the standard deviation of the residuals + to adjust plot limits. + + errorbar_kwargs : dict + Customize the data points plot by passing the given keyword + arguments to matplotlib.pyplot.errorbar. """ inst, z = self._parse_instrument_and_redshift(spec, redshift) - _, _, _, xz, yz, uncz = self._convert_spec_data(spec, z) - # Always use the current FWHM here (use_instrument_fwhm would - # overwrite the fitted value in the instrument overlap regions!) - astropy_model = self._construct_astropy_model( - inst, z, use_instrument_fwhm=False - ) - - enough_samples = max(1000, len(spec.wavelength)) + _, _, _, lam, flux, unc = self._convert_spec_data(spec, z) + enough_samples = max(10000, len(spec.wavelength)) + lam_mod = np.logspace(np.log10(min(lam)), np.log10(max(lam)), enough_samples) fig, axs = plt.subplots( ncols=1, nrows=2, - figsize=(15, 10), + figsize=(10, 10), gridspec_kw={"height_ratios": [3, 1]}, sharex=True, ) - PAHFITBase.plot(axs, xz, yz, uncz, astropy_model, model_samples=enough_samples) + + # spectrum and best fit model + ax = axs[0] + ax.set_yscale("linear") + ax.set_xscale("log") + ax.minorticks_on() + ax.tick_params( + axis="both", which="major", top="on", right="on", direction="in", length=10 + ) + ax.tick_params( + axis="both", which="minor", top="on", right="on", direction="in", length=5 + ) + + ext_model = None + has_att = "attenuation" in self.features["kind"] + has_abs = "absorption" in self.features["kind"] + if has_att: + row = self.features[self.features["kind"] == "attenuation"][0] + tau = row["tau"][0] + ext_model = S07_attenuation(tau_sil=tau)(lam_mod) + + if has_abs: + raise NotImplementedError( + "plotting absorption features not implemented yet" + ) + + if has_att or has_abs: + ax_att = ax.twinx() # axis for plotting the extinction curve + ax_att.tick_params(which="minor", direction="in", length=5) + ax_att.tick_params(which="major", direction="in", length=10) + ax_att.minorticks_on() + ax_att.plot(lam_mod, ext_model, "k--", alpha=0.5) + ax_att.set_ylabel("Attenuation") + ax_att.set_ylim(0, 1.1) + else: + ext_model = np.ones(len(lam_mod)) + + # Define legend lines + Leg_lines = [ + mpl.lines.Line2D([0], [0], color="k", linestyle="--", lw=2), + mpl.lines.Line2D([0], [0], color="#FE6100", lw=2), + mpl.lines.Line2D([0], [0], color="#648FFF", lw=2, alpha=0.5), + mpl.lines.Line2D([0], [0], color="#DC267F", lw=2, alpha=0.5), + mpl.lines.Line2D([0], [0], color="#785EF0", lw=2, alpha=1), + mpl.lines.Line2D([0], [0], color="#FFB000", lw=2, alpha=0.5), + ] + + # local utility + def tabulate_components(kind): + ss = {} + for name in self.features[self.features["kind"] == kind]["name"]: + ss[name] = self.tabulate( + inst, z, lam_mod, self.features["name"] == name + ) + return {name: s.flux.value for name, s in ss.items()} + + cont_y = np.zeros(len(lam_mod)) + if "dust_continuum" in self.features["kind"]: + # one plot for every component + for y in tabulate_components("dust_continuum").values(): + ax.plot(lam_mod, y * ext_model, "#FFB000", alpha=0.5) + # keep track of total continuum + cont_y += y + + if "starlight" in self.features["kind"]: + star_y = self.tabulate( + inst, z, lam_mod, self.features["kind"] == "starlight" + ).flux.value + ax.plot(lam_mod, star_y * ext_model, "#ffB000", alpha=0.5) + cont_y += star_y + + # total continuum + ax.plot(lam_mod, cont_y * ext_model, "#785EF0", alpha=1) + + # now plot the dust bands and lines + if "dust_feature" in self.features["kind"]: + for y in tabulate_components("dust_feature").values(): + ax.plot( + lam_mod, + (cont_y + y) * ext_model, + "#648FFF", + alpha=0.5, + ) + + if "line" in self.features["kind"]: + for name, y in tabulate_components("line").items(): + ax.plot( + lam_mod, + (cont_y + y) * ext_model, + "#DC267F", + alpha=0.5, + ) + if label_lines: + i = np.argmax(y) + # ignore out of range lines + if i > 0 and i < len(y) - 1: + w = lam_mod[i] + ax.text( + w, + y[i], + name, + va="center", + ha="center", + rotation="vertical", + bbox=dict(facecolor="white", alpha=0.75, pad=0), + ) + + ax.plot(lam_mod, self.tabulate(inst, z, lam_mod).flux.value, "#FE6100", alpha=1) + + # data + default_kwargs = dict( + fmt="o", + markeredgecolor="k", + markerfacecolor="none", + ecolor="k", + elinewidth=0.2, + capsize=0.5, + markersize=6, + ) + + ax.errorbar(lam, flux, yerr=unc, **(default_kwargs | errorbar_kwargs)) + + ax.set_ylim(0) + ax.set_ylabel(r"$\nu F_{\nu}$") + + ax.legend( + Leg_lines, + [ + "S07_attenuation", + "Spectrum Fit", + "Dust Features", + r"Atomic and $H_2$ Lines", + "Total Continuum Emissions", + "Continuum Components", + ], + prop={"size": 10}, + loc="best", + facecolor="white", + framealpha=1, + ncol=3, + ) + + # residuals = data in rest frame - (model evaluated at rest frame wavelengths) + res = flux - self.tabulate(inst, 0, lam).flux.value + std = np.nanstd(res) + ax = axs[1] + + ax.set_yscale("linear") + ax.set_xscale("log") + ax.tick_params( + axis="both", which="major", top="on", right="on", direction="in", length=10 + ) + ax.tick_params( + axis="both", which="minor", top="on", right="on", direction="in", length=5 + ) + ax.minorticks_on() + + # Custom X axis ticks + ax.xaxis.set_ticks( + [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 20, 25, 30, 40] + ) + + ax.axhline(0, linestyle="--", color="gray", zorder=0) + ax.plot( + lam, + res, + "ko", + fillstyle="none", + zorder=1, + markersize=errorbar_kwargs.get("markersize", None), + alpha=errorbar_kwargs.get("alpha", None), + linestyle="none", + ) + ax.set_ylim(-scalefac_resid * std, scalefac_resid * std) + ax.set_xlim(np.floor(np.amin(lam)), np.ceil(np.amax(lam))) + ax.set_xlabel(r"$\lambda$ [$\mu m$]") + ax.set_ylabel("Residuals [%]") + + # scalar x-axis marks + ax.xaxis.set_major_formatter(mpl.ticker.ScalarFormatter()) fig.subplots_adjust(hspace=0) + return fig def copy(self): """Copy the model. @@ -501,185 +757,232 @@ def tabulate( The flux model, evaluated at the given wavelengths, packaged as a Spectrum1D object. """ + z = 0 if redshift is None else redshift + + # decide which wavelength grid to use + if wavelengths is None: + ranges = instrument.wave_range(instrumentname) + if isinstance(ranges[0], float): + lam_min, lam_max = ranges + else: + # In case of multiple ranges (multiple segments), choose + # the min and max instead + lam_min = min(r[0] for r in ranges) + lam_max = max(r[1] for r in ranges) + + wfwhm = instrument.fwhm(instrumentname, lam_min, as_bounded=True)[0, 0] + lam = np.arange(lam_min, lam_max, wfwhm / 2) * u.micron + + elif isinstance(wavelengths, Spectrum1D): + lam = wavelengths.spectral_axis.to(u.micron) / (1 + z) + else: + # any other iterable will be accepted and converted to array + lam = np.asarray(wavelengths) * u.micron + # apply feature mask, make sub model, and set up functional if feature_mask is not None: - features_copy = self.features.copy() - features_to_use = features_copy[feature_mask] + features_to_use = self.features[feature_mask] else: features_to_use = self.features + + # if nothing is in range, return early with zeros + if len(features_to_use) == 0: + return Spectrum1D( + spectral_axis=lam, flux=np.zeros(lam.shape) * u.dimensionless_unscaled + ) + alt_model = Model(features_to_use) # Always use the current FWHM here (use_instrument_fwhm would # overwrite the value in the instrument overlap regions!) - flux_function = alt_model._construct_astropy_model( - instrumentname, redshift, use_instrument_fwhm=False - ) - # decide which wavelength grid to use - if wavelengths is None: - ranges = instrument.wave_range(instrumentname) - wmin = min(r[0] for r in ranges) - wmax = max(r[1] for r in ranges) - wfwhm = instrument.fwhm(instrumentname, wmin, as_bounded=True)[0, 0] - wav = np.arange(wmin, wmax, wfwhm / 2) * u.micron - elif isinstance(wavelengths, Spectrum1D): - wav = wavelengths.spectral_axis.to(u.micron) - else: - # any other iterable will be accepted and converted to array - wav = np.asarray(wavelengths) * u.micron + # need to wrap in try block to avoid bug: if all components are + # removed (because of wavelength range considerations), it won't work + try: + alt_model._set_up_fitter(instrumentname, z, use_instrument_fwhm=False) + fitter = alt_model.fitter + except PAHFITModelError: + return Spectrum1D( + spectral_axis=lam, flux=np.zeros(lam.shape) * u.dimensionless_unscaled + ) - # shift the "observed wavelength grid" to "physical wavelength grid" - wav /= 1 + redshift - flux_values = flux_function(wav.value) + flux_values = fitter.evaluate(lam.value) # apply unit stored in features table (comes from from last fit # or from loading previous result from disk) if "flux" not in self.features.meta["user_unit"]: flux_quantity = flux_values * u.dimensionless_unscaled else: - flux_quantity = flux_values * self.features.meta["user_unit"]["flux"] + user_unit = self.features.meta["user_unit"]["flux"] + flux_quantity = (flux_values * units.intensity).to(user_unit) - return Spectrum1D(spectral_axis=wav, flux=flux_quantity) + return Spectrum1D(spectral_axis=lam, flux=flux_quantity) - def _kludge_param_info(self, instrumentname, redshift, use_instrument_fwhm=True): - param_info = PAHFITBase.parse_table(self.features) - # edit line widths and drop lines out of range + def _excluded_features(self, instrumentname, redshift, lam_obs=None): + """Determine excluded features Based on instrument wavelength range. - # h2_info - param_info[2] = PAHFITBase.update_dictionary( - param_info[2], - instrumentname, - update_fwhms=use_instrument_fwhm, - redshift=redshift, - ) - # ion_info - param_info[3] = PAHFITBase.update_dictionary( - param_info[3], - instrumentname, - update_fwhms=use_instrument_fwhm, - redshift=redshift, - ) - # abs_info - param_info[4] = PAHFITBase.update_dictionary( - param_info[4], instrumentname, redshift + instrumentname : str + Qualified instrument name + + lam_obs : array + Optional observed wavelength grid. Exclude any lines and + dust features outside of this range. + + Returns + ------- + array of bool, same length as self.features, True where features + are far outside the wavelength range. + """ + lam_feature_obs = self.features["wavelength"]["val"] * (1 + redshift) + + # has wavelength and not within instrument range + is_outside = ~instrument.within_segment(lam_feature_obs, instrumentname) + + # also apply observed range if provided + if lam_obs is not None: + is_outside |= (lam_feature_obs < np.amin(lam_obs)) | ( + lam_feature_obs > np.amax(lam_obs) + ) + + # restriction on the kind of feature that can be excluded + excludable = ["line", "dust_feature", "absorption"] + is_excludable = np.logical_or.reduce( + [kind == self.features["kind"] for kind in excludable] ) - return param_info + return is_outside & is_excludable - def _construct_astropy_model( - self, instrumentname, redshift, use_instrument_fwhm=True + def _set_up_fitter( + self, instrumentname, redshift, lam=None, use_instrument_fwhm=True ): - """Convert the features table into a fittable model. + """Convert features table to Fitter instance, set self.fitter. - Some nuances in the behavior - - If a line has a fwhm set, it will be ignored, and replaced by - the calculated fwhm provided by the instrument model. - - If a line has been masked by _parse_astropy_result, and this - function is called again, those masks will be ignored, as the - data range might have changed. + For every row of the features table, calls a function of Fitter + API to register an appropriate component. Finalizes the Fitter + at the end (details of this step depend on the Fitter subclass). - TODO: Make sure the features outside of the data range are - removed. The instrument-based feature check is done in - _kludge_param_info(), but the observational data might only - cover a part of the instrument range. + Any unit conversions and model-specific things need to happen + BEFORE giving them to the fitters. + - The instrument-derived FWHM is determined here using the + instrument model (the Fitter does not need to know about this + detail). + - Features outside the appropriate wavelength range should not + be added to the Fitter: the "trimming" is done here, using the + given wavelength range lam (optional). - """ - param_info = self._kludge_param_info( - instrumentname, redshift, use_instrument_fwhm - ) - return PAHFITBase.model_from_param_info(param_info) + TODO: flags to indicate which features were excluded. - def _parse_astropy_result(self, astropy_model): - """Store the result of the astropy fit into the features table. + Parameters + ---------- - Every relevant value inside the astropy model, is written to the - right position in the features table. + instrumentname and redshift : needed to apply the instrument + model and to determine which feature to exclude - For the unresolved lines, the widths are calculated by the - instrument model, or fitted when these lines are in a spectral - overlap region. The calculated or fitted result is written to - the fwhm field of the table. When a new model is constructed - from the features table, this fwhm value will be ignored. + lam : array of observed wavelengths + Used to exclude features from the model based on the actual + observed data given. - For features that do not correspond to the data range, all - parameter values will be masked. Their numerical values remain - accessible by '.data' on the masked entity. This way, We still - keep their parameter values around (as opposed to removing the - rows entirely). When data with a larger range are passed for - another fitting call, those features can be unmasked if - necessary. + use_instrument_fwhm : bool + When set to False, the instrument model is not used and the + FWHM values are taken from the features table as-is. This is + the current workaround to fit the widths of lines, until the + "physical" and "instrumental" widths are conceptually + separated (see issue #293). """ - if len(self.features) < 2: - # Plotting and tabulating works fine, but the code below - # will not work with only one component. This can be - # addressed later, when the internal API is made agnostic of - # the fitting backend (astropy vs our own). - raise PAHFITModelError("Fit with fewer than 2 components not allowed!") - - # Some translation rules between astropy model components and - # feature table names and values. - - # these have the same value but different (fixed) names - param_name_equivalent = { - "temperature": "temperature", - "fwhm": "fwhm", - "x_0": "wavelength", - "mean": "wavelength", - "tau_sil": "tau", - } - - def param_conversion(features_kind, param_name, param_value): - # default conversion - if param_name in param_name_equivalent: - new_name = param_name_equivalent[param_name] - new_value = param_value - # certain types of features use tau instead of amplitude - elif param_name == "amplitude": - if features_kind in ["starlight", "dust_continuum", "absorption"]: - new_name = "tau" - else: - new_name = "power" - new_value = param_value - # convert stddev to fwhm - elif param_name == "stddev": - new_name = "fwhm" - new_value = param_value * 2.355 + # Fitting implementation can be changed by choosing another + # Fitter class. TODO: make this configurable. + self.fitter = APFitter() + + excluded = self._excluded_features(instrumentname, redshift, lam) + self.enabled_features = self.features["name"][~excluded] + + def cleaned(features_tuple3): + val = features_tuple3[0] + if bounded_is_fixed(features_tuple3): + return val else: - raise NotImplementedError( - f"no conversion rule for model parameter {param_name}" - ) - return new_name, new_value + vmin = -np.inf if np.isnan(features_tuple3[1]) else features_tuple3[1] + vmax = np.inf if np.isnan(features_tuple3[2]) else features_tuple3[2] + return np.array([val, vmin, vmax]) - # Go over all features. - for row in self.features: + for row in self.features[~excluded]: + kind = row["kind"] name = row["name"] - if name in astropy_model.submodel_names: - # undo any previous masking that might have occured - self.features.unmask_feature(name) - - # copy or translate, and store the parameters - component = astropy_model[name] - for param_name in component.param_names: - param_value = getattr(component, param_name).value - col_name, col_value = param_conversion( - row["kind"], param_name, param_value - ) - row[col_name][0] = col_value - # for the unresolved lines, indicate when the line fwhm was made non-fixed - if row["kind"] == "line" and col_name == "fwhm": - row["fwhm"].mask[1:] = component.fixed[param_name] + if kind == "starlight": + self.fitter.add_feature_starlight( + name, cleaned(row["temperature"]), cleaned(row["tau"]) + ) + + elif kind == "dust_continuum": + self.fitter.add_feature_dust_continuum( + name, cleaned(row["temperature"]), cleaned(row["tau"]) + ) + + elif kind == "line": + # be careful with lines that have masked FWHM values here + if use_instrument_fwhm or row["fwhm"] is np.ma.masked: + # One caveat here: redshift. We do the necessary + # adjustment as follows : 1. shift to observed wav; + # 2. evaluate fwhm at oberved wav; 3. shift back to + # rest frame wav (width in rest frame will be + # narrower than observed width) + lam_obs = row["wavelength"]["val"] * (1.0 + redshift) + # returned value is tuple (value, min, max). And + # min/max are already masked in case of fixed value + # (output of instrument.resolution is designed to be + # very similar to an entry in the features table) + calculated_fwhm = instrument.fwhm( + instrumentname, lam_obs, as_bounded=True + )[0] / (1.0 + redshift) + + # decide if scalar (fixed) or tuple (fitted fwhm + # between upper and lower fwhm limits, happens in + # case of overlapping instruments) + if calculated_fwhm[1] is np.ma.masked: + fwhm = calculated_fwhm[0] + else: + fwhm = calculated_fwhm + + else: + # if instrument model is not to be used, just take + # the value as is specified in the Features table + fwhm = cleaned(row["fwhm"]) + + self.fitter.add_feature_line( + name, cleaned(row["power"]), cleaned(row["wavelength"]), fwhm + ) + + elif kind == "dust_feature": + self.fitter.add_feature_dust_feature( + name, + cleaned(row["power"]), + cleaned(row["wavelength"]), + cleaned(row["fwhm"]), + ) + + elif kind == "attenuation": + self.fitter.add_feature_attenuation(name, cleaned(row["tau"])) + + elif kind == "absorption": + self.fitter.add_feature_absorption( + name, + cleaned(row["tau"]), + cleaned(row["wavelength"]), + cleaned(row["fwhm"]), + ) + else: - # signal that it was not fit by masking the feature - self.features.mask_feature(name) + raise PAHFITModelError( + f"Model components of kind {kind} are not implemented!" + ) + + self.fitter.finalize() @staticmethod def _parse_instrument_and_redshift(spec, redshift): - """Small utility to grab instrument and redshift from either - Spectrum1D metadata, or from arguments. - - """ + """Get instrument redshift from Spectrum1D metadata or arguments.""" # the rest of the implementation doesn't like Quantity... z = spec.redshift.value if redshift is None else redshift if z is None: diff --git a/pahfit/scripts/plot_pahfit.py b/pahfit/scripts/plot_pahfit.py index d5bb8532..a8765058 100755 --- a/pahfit/scripts/plot_pahfit.py +++ b/pahfit/scripts/plot_pahfit.py @@ -6,11 +6,8 @@ import matplotlib as mpl from pahfit.model import Model -from pahfit.base import PAHFITBase from pahfit.helpers import read_spectrum -from astropy import units as u - def initialize_parser(): """ @@ -78,14 +75,6 @@ def main(): def default_layout_plot(spec, model, scalefac_resid): - """ - Returns - ------- - fig : Figure object - - """ - - # plot result fontsize = 18 font = {"size": fontsize} mpl.rc("font", **font) @@ -96,29 +85,7 @@ def default_layout_plot(spec, model, scalefac_resid): mpl.rc("xtick.minor", size=3, width=1) mpl.rc("ytick.minor", size=3, width=1) - fig, axs = plt.subplots( - ncols=1, - nrows=2, - figsize=(15, 10), - gridspec_kw={"height_ratios": [3, 1]}, - sharex=True, - ) - - enough_samples = max(1000, len(spec.wavelength)) - - PAHFITBase.plot( - axs, - spec.wavelength.to(u.micron).value, - spec.flux, - spec.uncertainty.array, - model._construct_astropy_model( - spec.meta["instrument"], spec.redshift, use_instrument_fwhm=False - ), - model_samples=enough_samples, - scalefac_resid=scalefac_resid, - ) - - # use the whitespace better + fig = model.plot(spec, scalefac_resid=scalefac_resid) fig.subplots_adjust(hspace=0) return fig diff --git a/pahfit/tests/test_feature_parsing.py b/pahfit/tests/test_feature_parsing.py index ee0e39e5..d658a689 100644 --- a/pahfit/tests/test_feature_parsing.py +++ b/pahfit/tests/test_feature_parsing.py @@ -7,6 +7,7 @@ def test_feature_parsing(): """ + Goal ---- Test if the model is built successfully with certain features removed @@ -21,7 +22,8 @@ def test_feature_parsing(): Desired behavior ---------------- - The PAHFITBase instance is generated correctly, without crashing. + The Fitter instance underlying model is generated correctly, without + crashing. Functions that depend on specific model contents (lines, dust features, ...) can deal with those feature not being there. @@ -38,8 +40,8 @@ def test_feature_parsing(): def test_parsing(features_edit): m = Model(features_edit) - amodel = m._construct_astropy_model(instrumentname, 0) - m._parse_astropy_result(amodel) + m._set_up_fitter(instrumentname, 0) + m._ingest_fit_result_to_features() # Case 0: the whole table test_parsing(features) diff --git a/pahfit/tests/test_fitting_spitzer.py b/pahfit/tests/test_fitting_spitzer.py index d56975bf..56e4aa56 100644 --- a/pahfit/tests/test_fitting_spitzer.py +++ b/pahfit/tests/test_fitting_spitzer.py @@ -66,7 +66,7 @@ def test_fitting_m101(): try: np.testing.assert_allclose( - model.astropy_result.parameters, expvals, rtol=1e-6, atol=1e-6 + model.fitter.model, expvals, rtol=1e-6, atol=1e-6 ) except AssertionError as error: print(error) diff --git a/pahfit/tests/test_model_impl.py b/pahfit/tests/test_model_impl.py index fce3acc3..6007c175 100644 --- a/pahfit/tests/test_model_impl.py +++ b/pahfit/tests/test_model_impl.py @@ -10,9 +10,10 @@ def assert_features_table_equality(features1, features2): for string_col in ["name", "group", "kind", "model", "geometry"]: assert (features1[string_col] == features2[string_col]).all() for param_col in ["temperature", "tau", "wavelength", "power", "fwhm"]: - np.testing.assert_allclose( - features1[param_col], features2[param_col], rtol=1e-6, atol=1e-6 - ) + for k in ("val", "min", "max"): + np.testing.assert_allclose( + features1[param_col][k], features2[param_col][k], rtol=1e-6, atol=1e-6 + ) def default_spec_and_model_fit(fit=True): @@ -35,14 +36,17 @@ def test_feature_table_model_conversion(): # results were then written to model.features. If everything went # correct, reconstructing the model from model.features should # result in the exact same model. - fit_result = model.astropy_result - reconstructed_fit_result = model._construct_astropy_model( + + fit_result = model.fitter + model._set_up_fitter( instrumentname=spec.meta["instrument"], redshift=0, use_instrument_fwhm=False ) - for p in fit_result.param_names: - p1 = getattr(fit_result, p) - p2 = getattr(reconstructed_fit_result, p) - assert p1 == p2 + reconstructed_fit_result = model.fitter + for name in model.enabled_features: + par_dict1 = fit_result.get_result(name) + par_dict2 = reconstructed_fit_result.get_result(name) + for key in par_dict1: + assert par_dict1[key] == par_dict2[key] def test_model_edit(): @@ -57,19 +61,23 @@ def test_model_edit(): # edit the same parameter in the copy newT = 123 - model_to_edit.features.loc[feature][col][0] = newT + + i = np.where(model_to_edit.features["name"] == feature)[0] + model_to_edit.features[col]["val"][i] = newT # make sure the original value is still the same - assert model.features.loc[feature][col][0] == originalT + j = np.where(model.features["name"] == feature)[0] + assert model.features[col]["val"][j] == originalT # construct astropy model with dummy instrument - astropy_model_edit = model_to_edit._construct_astropy_model( + model_to_edit._set_up_fitter( instrumentname="spitzer.irs.*", redshift=0 ) + fitter_edit = model_to_edit.fitter # Make sure the change is reflected in this model. Very handy that # we can access the right component by the feature name! - assert astropy_model_edit[feature].temperature == newT + assert fitter_edit.get_result(feature)["temperature"] == newT def test_model_tabulate(): diff --git a/pahfit/units.py b/pahfit/units.py new file mode 100644 index 00000000..103e42f3 --- /dev/null +++ b/pahfit/units.py @@ -0,0 +1,17 @@ +import astropy.units as u +from astropy.units import CompositeUnit + + +# Working/default unitParameter default units: flux density/intensity/power +# These are PAHFITs default science packs parameter and output units +temperature = u.K +wavelength = u.um +flux_density = u.mJy +flux_power = CompositeUnit(1e-22, (u.W, u.m), (1, -2)) +solid_angle = u.sr +intensity = u.MJy / u.sr +intensity_power = CompositeUnit(1e-10, (u.W, u.m, u.sr), (1, -2, -1)) + +# Note: integrated power units of 1e-22 W/m^2 (from flux) corresponds +# to the unit 1e-10 W/m^2/sr (from intensity) if it occurs uniformly +# over a solid angle 0.21" on a side (about a small JWST IFU pixel)