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 0a55a881..00000000 --- a/pahfit/base.py +++ /dev/null @@ -1,567 +0,0 @@ -import numpy as np - -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 - - 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/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 4bb56f7c..0efed64b 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -137,7 +137,7 @@ class Features(Table): _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', 'kind', 'geometry', 'model')) # str attributes (no bounds) - _bounds_dtype = np.dtype([("val", "f4"), ("min", "f4"), ("max", "f4")]) # bounded param type + _bounds_dtype = np.dtype([("val", float), ("min", float), ("max", float)]) # bounded param type @classmethod def read(cls, file, *args, **kwargs): diff --git a/pahfit/fitters/__init__.py b/pahfit/fitters/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/pahfit/component_models.py b/pahfit/fitters/ap_components.py similarity index 100% rename from pahfit/component_models.py rename to pahfit/fitters/ap_components.py diff --git a/pahfit/fitters/ap_fitter.py b/pahfit/fitters/ap_fitter.py new file mode 100644 index 00000000..92760bf0 --- /dev/null +++ b/pahfit/fitters/ap_fitter.py @@ -0,0 +1,409 @@ +from pahfit.fitters.fitter import Fitter +from pahfit.errors import PAHFITModelError +from pahfit.fitters.ap_components import ( + BlackBody1D, + ModifiedBlackBody1D, + S07_attenuation, + att_Drude1D, + Drude1D, +) +from astropy.modeling.functional_models import Gaussian1D +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, + ["amplitude", "mean", "stddev"], + [power, wavelength, fwhm / 2.355], + ) + self._add_component(Gaussian1D, **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, ["amplitude", "x_0", "fwhm"], [power, wavelength, fwhm] + ) + self._add_component(Drude1D, **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). + + NOTE: for now, the return units for "power" are (flux unit) x + (micron). Still ambiguous, because the flux unit could be flux + density or intensity. + + 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. {'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.amplitude.value, + "wavelength": component.mean.value, + "fwhm": component.stddev.value * 2.355, + } + elif c_type == "dust_feature": + return { + "power": component.amplitude.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 fc84ee51..f6bcabc2 100644 --- a/pahfit/helpers.py +++ b/pahfit/helpers.py @@ -1,16 +1,13 @@ import os from importlib import resources -from pahfit.component_models import BlackBody1D, S07_attenuation from pahfit import units from specutils import Spectrum1D -from astropy.modeling.physical_models import Drude1D -from astropy.modeling.functional_models import Gaussian1D 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): @@ -90,108 +87,13 @@ def read_spectrum(specfile, format=None): # 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 + 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)) + 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 - -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 diff --git a/pahfit/model.py b/pahfit/model.py index 1dd15a57..83380b07 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -1,19 +1,18 @@ from specutils import Spectrum1D from astropy import units as u 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 import units -from pahfit.features.util import bounded_is_fixed +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, S07_attenuation +from pahfit.fitters.ap_components import BlackBody1D, S07_attenuation +from pahfit.fitters.ap_fitter import APFitter class Model: @@ -182,17 +181,17 @@ 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) + wmin = min(lam) + wmax = max(lam) # 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): @@ -229,14 +228,14 @@ def dust_continuum_guess(row): fmax_lam = 2898.0 / temp bb = BlackBody1D(1, temp) if fmax_lam >= wmin and fmax_lam <= wmax: - w_ref = fmax_lam + lam_ref = fmax_lam elif fmax_lam > wmax: - w_ref = wmax + lam_ref = wmax else: - w_ref = wmin + lam_ref = wmin - 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) @@ -257,12 +256,12 @@ def amp_guess(row, fwhm): 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_window = (lam > wmin) & (lam < wmax) + 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]) + power_guess = 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: @@ -273,7 +272,7 @@ def amp_guess(row, fwhm): # Same logic as in the old function: just use same amp for all # dust features. - some_flux = 0.5 * np.median(yz) + some_flux = 0.5 * np.median(flux) loop_over_non_fixed("dust_feature", "power", lambda row: some_flux) if integrate_line_flux: @@ -284,9 +283,23 @@ def amp_guess(row, fwhm): else: loop_over_non_fixed("line", "power", lambda row: some_flux) - # 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): @@ -303,21 +316,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 """ if not spec.flux.unit.is_equivalent(units.intensity): - raise PAHFITModelError("For now, PAHFIT only supports intensity units, i.e. convertible to MJy / sr.") - y = spec.flux.to(units.intensity).value - x = spec.spectral_axis.to(u.micron).value - unc = (spec.uncertainty.array * spec.flux.unit).to(units.intensity).value + 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, @@ -378,7 +393,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 @@ -387,30 +402,45 @@ 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) - self._parse_astropy_result(self.astropy_result) + def _ingest_fit_result_to_features(self): + """Copy the results from the Fitter to the features table - def info(self): - """Print out the last fit results.""" - print(self.astropy_result) + This is a utility method, executed only at the end of fit(), + where Fitter.fit() has been applied. + + """ + # 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, @@ -454,9 +484,9 @@ def plot( """ inst, z = self._parse_instrument_and_redshift(spec, redshift) - _, _, _, xz, yz, uncz = self._convert_spec_data(spec, z) + _, _, _, lam, flux, unc = self._convert_spec_data(spec, z) enough_samples = max(10000, len(spec.wavelength)) - x_mod = np.logspace(np.log10(min(xz)), np.log10(max(xz)), enough_samples) + lam_mod = np.logspace(np.log10(min(lam)), np.log10(max(lam)), enough_samples) fig, axs = plt.subplots( ncols=1, @@ -484,7 +514,7 @@ def plot( if has_att: row = self.features[self.features["kind"] == "attenuation"][0] tau = row["tau"][0] - ext_model = S07_attenuation(tau_sil=tau)(x_mod) + ext_model = S07_attenuation(tau_sil=tau)(lam_mod) if has_abs: raise NotImplementedError( @@ -496,11 +526,11 @@ def plot( 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(x_mod, ext_model, "k--", alpha=0.5) + 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(x_mod)) + ext_model = np.ones(len(lam_mod)) # Define legend lines Leg_lines = [ @@ -516,32 +546,34 @@ def plot( def tabulate_components(kind): ss = {} for name in self.features[self.features["kind"] == kind]["name"]: - ss[name] = self.tabulate(inst, z, x_mod, self.features["name"] == 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(x_mod)) + 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(x_mod, y * ext_model, "#FFB000", alpha=0.5) + 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, x_mod, self.features["kind"] == "starlight" + inst, z, lam_mod, self.features["kind"] == "starlight" ).flux.value - ax.plot(x_mod, star_y * ext_model, "#ffB000", alpha=0.5) + ax.plot(lam_mod, star_y * ext_model, "#ffB000", alpha=0.5) cont_y += star_y # total continuum - ax.plot(x_mod, cont_y * ext_model, "#785EF0", alpha=1) + 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( - x_mod, + lam_mod, (cont_y + y) * ext_model, "#648FFF", alpha=0.5, @@ -550,7 +582,7 @@ def tabulate_components(kind): if "line" in self.features["kind"]: for name, y in tabulate_components("line").items(): ax.plot( - x_mod, + lam_mod, (cont_y + y) * ext_model, "#DC267F", alpha=0.5, @@ -559,7 +591,7 @@ def tabulate_components(kind): i = np.argmax(y) # ignore out of range lines if i > 0 and i < len(y) - 1: - w = x_mod[i] + w = lam_mod[i] ax.text( w, y[i], @@ -570,7 +602,7 @@ def tabulate_components(kind): bbox=dict(facecolor="white", alpha=0.75, pad=0), ) - ax.plot(x_mod, self.tabulate(inst, z, x_mod).flux.value, "#FE6100", alpha=1) + ax.plot(lam_mod, self.tabulate(inst, z, lam_mod).flux.value, "#FE6100", alpha=1) # data default_kwargs = dict( @@ -583,7 +615,7 @@ def tabulate_components(kind): markersize=6, ) - ax.errorbar(xz, yz, yerr=uncz, **(default_kwargs | errorbar_kwargs)) + ax.errorbar(lam, flux, yerr=unc, **(default_kwargs | errorbar_kwargs)) ax.set_ylim(0) ax.set_ylabel(r"$\nu F_{\nu}$") @@ -606,7 +638,7 @@ def tabulate_components(kind): ) # residuals = data in rest frame - (model evaluated at rest frame wavelengths) - res = yz - self.tabulate(inst, 0, xz).flux.value + res = flux - self.tabulate(inst, 0, lam).flux.value std = np.nanstd(res) ax = axs[1] @@ -627,7 +659,7 @@ def tabulate_components(kind): ax.axhline(0, linestyle="--", color="gray", zorder=0) ax.plot( - xz, + lam, res, "ko", fillstyle="none", @@ -637,7 +669,7 @@ def tabulate_components(kind): linestyle="none", ) ax.set_ylim(-scalefac_resid * std, scalefac_resid * std) - ax.set_xlim(np.floor(np.amin(xz)), np.ceil(np.amax(xz))) + ax.set_xlim(np.floor(np.amin(lam)), np.ceil(np.amax(lam))) ax.set_xlabel(r"$\lambda$ [$\mu m$]") ax.set_ylabel("Residuals [%]") @@ -707,185 +739,231 @@ 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): + wmin, wmax = ranges + else: + # In case of multiple ranges (multiple segments), choose + # the min and max instead + 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] + lam = np.arange(wmin, wmax, 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/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():