From c3abd93c82ffe571151f3ca405909ed4f3a18ea1 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Mon, 13 May 2024 15:31:42 -0400 Subject: [PATCH 01/27] Add Fitter API and APFitter implementation --- pahfit/apfitter.py | 420 +++++++++++++++++++++++++++++++++++++++++++++ pahfit/fitter.py | 203 ++++++++++++++++++++++ 2 files changed, 623 insertions(+) create mode 100644 pahfit/apfitter.py create mode 100644 pahfit/fitter.py diff --git a/pahfit/apfitter.py b/pahfit/apfitter.py new file mode 100644 index 0000000..4a3d051 --- /dev/null +++ b/pahfit/apfitter.py @@ -0,0 +1,420 @@ +from pahfit.fitter import Fitter +from pahfit.errors import PAHFITModelError +from pahfit.component_models 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 register_() functions to start + setting up a model, then call finalize_model(). + + """ + self.clear() + + def clear(self): + """Reset model. + + After reset, register_() and finalize_model() can be used again. + + """ + self.additive_components = [] + self.multiplicative_components = [] + self.component_types = {} + self.model = None + self.message = None + + def components(self): + """Return list of component names. + + Only works after finalize_model(). + + """ + if hasattr(self.model, "submodel_names"): + return self.model.submodel_names + else: + # single-component edge case + return [self.model.name] + + def finalize_model(self): + """Sum the registered components into one CompoundModel. + + To be called after a series of "register_()" 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 _register_component(self, astropy_model_class, multiplicative=False, **kwargs): + """Register any component in a generic way. + + 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 register_() 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 register_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.component_types[name] = "starlight" + kwargs = self._astropy_model_kwargs( + name, ["temperature", "amplitude"], [temperature, tau] + ) + self._register_component(BlackBody1D, **kwargs) + + def register_dust_continuum(self, name, temperature, tau): + """Register a ModifiedBlackBody1D. + + Analogous. Temperature and tau are used as temperature and + amplitude + + """ + self.component_types[name] = "dust_continuum" + kwargs = self._astropy_model_kwargs( + name, ["temperature", "amplitude"], [temperature, tau] + ) + self._register_component(ModifiedBlackBody1D, **kwargs) + + def register_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.component_types[name] = "line" + + kwargs = self._astropy_model_kwargs( + name, + ["amplitude", "mean", "stddev"], + [power, wavelength, np.array([x for x in fwhm]) / 2.355], + ) + self._register_component(Gaussian1D, **kwargs) + + def register_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.component_types[name] = "dust_feature" + kwargs = self._astropy_model_kwargs( + name, ["amplitude", "x_0", "fwhm"], [power, wavelength, fwhm] + ) + self._register_component(Drude1D, **kwargs) + + def register_attenuation(self, name, tau): + """Register the S07 attenuation component. + + Analogous. Uses tau as tau_sil for S07_attenuation. Is + multiplicative. + + """ + self.component_types[name] = "attenuation" + kwargs = self._astropy_model_kwargs(name, ["tau_sil"], [tau]) + self._register_component(S07_attenuation, multiplicative=True, **kwargs) + + def register_absorption(self, name, tau, wavelength, fwhm): + """Register an absorbing Drude1D component. + + Analogous. Is multiplicative. + + """ + self.component_types[name] = "absorption" + kwargs = self._astropy_model_kwargs( + name, ["tau", "x_0", "fwhm"], [tau, wavelength, fwhm] + ) + self._register_component(att_Drude1D, multiplicative=True, **kwargs) + + def evaluate_model(self, xz): + """Evaluate internal astropy model with its current parameters. + + Parameters + ---------- + xz : array + Rest frame wavelengths in micron + + Returns + ------- + yz : array + Rest frame flux in internal units + """ + return self.model(xz) + + def fit(self, xz, yz, uncz, 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 (yz) 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 wavelenght unit). + + Parameters + ---------- + xz : array + Rest frame wavelengths in micron + + yz : array + Rest frame flux in internal units. + + uncz : array + Uncertainty on rest frame flux. Same units as yz. + + """ + # clean, because astropy does not like nan + w = 1 / uncz + + # make sure there are no zero uncertainties either + mask = np.isfinite(xz) & np.isfinite(yz) & np.isfinite(w) + + self.fit_info = [] + + fit = LevMarLSQFitter(calc_uncertainties=True) + astropy_result = fit( + self.model, + xz[mask], + yz[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 register_*() calls + made during setup. See also Fitter.components(). + + 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.component_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, value_tuples): + """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"] + + value_tuples : list of (array of size 3) + One for each param name, each in the format of [value, + min_bound, max_bound], i.e. in the format as stored in the + Features table. This means that [value, masked, masked] will + result in a fixed parameter. + + 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, value_tuple in zip(param_names, value_tuples): + kwargs[param_name] = value_tuple[0] + + # astropy does not like numpy bools, so we do this silly + # conversion. + is_fixed = np.isnan(value_tuple[1]) and np.isnan(value_tuple[2]) + kwargs["fixed"][param_name] = True if is_fixed else False + + # For the limits, use 0 if fixed, the raw values if + # variable, and None if variable but unbounded. + min_max = value_tuple[1], value_tuple[2] + limits = [0 if is_fixed else v for v in min_max] + kwargs["bounds"][param_name] = [None if np.isinf(x) else x for x in limits] + + return kwargs diff --git a/pahfit/fitter.py b/pahfit/fitter.py new file mode 100644 index 0000000..5c8b657 --- /dev/null +++ b/pahfit/fitter.py @@ -0,0 +1,203 @@ +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.). + + For the model setup, multiple functions are used, so a few notes are + provided here. There is one function per type of component supported + by PAHFIT, and the arguments of these functions will ask for + different "standard" PAHFIT numbers, i.e. those from the Features + table. These functions have the same signature between all Fitter + implementations, so that the Model class can use a single + implementation to set up the Fitter. The Model has access to the + Features table and the instrument model, and needs to set up Fitter + with the correct initial values, bounds, and "fixed" flags (e.g. + setting a fixed FWHM based on the instrument for the lines). After + all the components have been added, the finalize_model() function + can be called to finish setting up the internal astropy model. After + this has finished, fit() can be called to apply the model and the + astropy fitter to the data. + + """ + + @abstractmethod + def clear(self): + """Reset model. + + After reset, register_() and finalize_model() can be used again. + + """ + pass + + @abstractmethod + def components(self): + """Return list of features. + + Only works after finalize_model(). Will return the names passed + using the register functions. + + """ + pass + + @abstractmethod + def finalize_model(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 register_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 + Blackbody temperature, given as [value, lower_bound, + upper_bound]. The bounds assume the same system as the + features table: if they are masked, the parameter will be + fixed, while +inf or -inf mean unbounded. + + tau : array of size 3 + Analogously, used as power. + + """ + pass + + @abstractmethod + def register_dust_continuum(self, name, temperature, tau): + """Register a dust continuum feature.""" + pass + + @abstractmethod + def register_line(self, name, power, wavelength, fwhm): + """Register an emission line feature. + + Typically a Gaussian profile. + + """ + pass + + @abstractmethod + def register_dust_feature(self, name, power, wavelength, fwhm): + """Register a dust feature. + + Typically a Drude profile. + + """ + pass + + @abstractmethod + def register_attenuation(self, name, tau): + """Register the S07 attenuation component. + + Other types of attenuation might be possible in the future. Is + multiplicative. + + """ + pass + + @abstractmethod + def register_absorption(self, name, tau, wavelength, fwhm): + """Register an absorption feature. + + Typically a Drude profile. Is multiplicative. + + """ + pass + + @abstractmethod + def evaluate_model(self, xz): + """Evaluate the model at the given wavelengths. + + Parameters + ---------- + xz : array + Rest frame wavelengths in micron + + Returns + ------- + yz : array + Rest frame flux in internal units + + """ + pass + + @abstractmethod + def fit(self, xz, yz, uncz, 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 + ---------- + xz : array + Rest frame wavelengths in micron + + yz : array + Rest frame flux in internal units. + + uncz : array + Uncertainty on rest frame flux. Same units as yz. + + """ + 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 register_() calls + made during setup. See also Fitter.components(). + + 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'} + + """ + pass From a73dcb4b8a662659b2010c60f93d24e586013a70 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Mon, 13 May 2024 20:04:34 -0400 Subject: [PATCH 02/27] Delete base.py and use Fitter API in Model class implementation --- pahfit/base.py | 567 --------------------------- pahfit/model.py | 364 +++++++++-------- pahfit/tests/test_feature_parsing.py | 4 +- pahfit/tests/test_fitting_spitzer.py | 2 +- pahfit/tests/test_model_impl.py | 23 +- 5 files changed, 218 insertions(+), 742 deletions(-) delete mode 100644 pahfit/base.py diff --git a/pahfit/base.py b/pahfit/base.py deleted file mode 100644 index 0a55a88..0000000 --- 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/model.py b/pahfit/model.py index 1dd15a5..5b77204 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -1,19 +1,19 @@ 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.fitter import Fitter +from pahfit.apfitter import APFitter class Model: @@ -284,9 +284,19 @@ 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 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 np.ma.is_masked(row["fwhm"]): + 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): @@ -387,30 +397,54 @@ 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.fitter = self._set_up_fitter( + inst, z, x=x, use_instrument_fwhm=use_instrument_fwhm ) - self.fit_info = fit.fit_info + self.fitter.fit(xz, yz, uncz, maxiter=maxiter) + + # copy the fit results to the features table + self._ingest_fit_result_to_features(self.fitter) + if verbose: - print(fit.fit_info["message"]) + print(self.fitter.message) + + def _ingest_fit_result_to_features(self, fitter: Fitter): + """Copy the results from a 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 the Fitter is passed after Fitter.fit() has been applied. + Passing a different Fitter object could be useful for + testing. - 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 _construct_model. Having an + # ENABLED/DISABLED flag for every feature would be a nice + # alternative (and clear for the user). + + self.features.meta["fitter_message"] = fitter.message + + for name in fitter.components(): + for column, value in 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) + print(f'ingested {name} {column} as {value:6e}') + print('value in table ', self.features[column][i]) + if not self.features[column]["val"][i] == value: + print("assigment went wrong") + raise PAHFITModelError + 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, @@ -707,25 +741,18 @@ def tabulate( The flux model, evaluated at the given wavelengths, packaged as a Spectrum1D object. """ - # 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] - else: - features_to_use = self.features - 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 - ) + z = 0 if redshift is None else redshift # 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) + 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] wav = np.arange(wmin, wmax, wfwhm / 2) * u.micron elif isinstance(wavelengths, Spectrum1D): @@ -734,158 +761,171 @@ def tabulate( # any other iterable will be accepted and converted to array wav = np.asarray(wavelengths) * u.micron - # shift the "observed wavelength grid" to "physical wavelength grid" - wav /= 1 + redshift - flux_values = flux_function(wav.value) + # 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] + 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=wav, flux=np.zeros(wav.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!) + + # need to wrap in try block to avoid bug: if all components are + # removed (because of wavelength range considerations), it won't work + try: + fitter = alt_model._set_up_fitter( + instrumentname, z, use_instrument_fwhm=False + ) + except PAHFITModelError: + return Spectrum1D( + spectral_axis=wav, flux=np.zeros(wav.shape) * u.dimensionless_unscaled + ) + + # shift the "observed wavelength grid" to "physical wavelength grid + wav /= 1 + z + flux_values = fitter.evaluate_model(wav.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.internal_flux_unit(user_unit)).to( + user_unit + ) return Spectrum1D(spectral_axis=wav, 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, x=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 - return param_info + x : array + Optional observed wavelength grid. Exclude any lines and + dust features outside of this range. - def _construct_astropy_model( - self, instrumentname, redshift, use_instrument_fwhm=True - ): - """Convert the features table into a fittable model. + Returns + ------- + array of bool, same length as self.features, True where features + are far outside the wavelength range. + """ + observed_wavs = self.features["wavelength"]["val"] * (1 + redshift) - 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. + # has wavelength and not within instrument range + is_outside = ~instrument.within_segment(observed_wavs, instrumentname) - 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. + # also apply observed range if provided + if x is not None: + is_outside |= (observed_wavs < np.amin(x)) | (observed_wavs > np.amax(x)) - """ - param_info = self._kludge_param_info( - instrumentname, redshift, use_instrument_fwhm + # 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 PAHFITBase.model_from_param_info(param_info) - def _parse_astropy_result(self, astropy_model): - """Store the result of the astropy fit into the features table. + return is_outside & is_excludable - Every relevant value inside the astropy model, is written to the - right position in the features table. + def _set_up_fitter( + self, instrumentname, redshift, x=None, use_instrument_fwhm=True + ): + """Convert features table to Fitter instance. + + 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). - 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. + 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 xz (optional). - 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. + TODO: flags to indicate which features were excluded. + + Returns + ------- + Fitter """ - 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" + # Fitting implementation can be changed by choosing another + # Fitter class. TODO: make this configurable. + fitter = APFitter() + + excluded = self._excluded_features(instrumentname, redshift, x) + + def array3(features_tuple3): + return np.array([x for x in features_tuple3]) + + for row in self.features[~excluded]: + kind = row["kind"] + name = row["name"] + + if kind == "starlight": + fitter.register_starlight(name, row["temperature"], row["tau"]) + + elif kind == "dust_continuum": + fitter.register_dust_continuum(name, row["temperature"], row["tau"]) + + elif kind == "line": + if use_instrument_fwhm: + # one caveat here: redshift. Correct way to do it: + # 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) + w_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) + fwhm = instrument.fwhm(instrumentname, w_obs, as_bounded=True)[ + 0 + ] / (1.0 + redshift) + fwhm = np.ma.filled(fwhm, np.nan) 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 - else: - raise NotImplementedError( - f"no conversion rule for model parameter {param_name}" + fwhm = row["fwhm"] + fitter.register_line(name, row["power"], row["wavelength"], fwhm) + + elif kind == "dust_feature": + fitter.register_dust_feature( + name, row["power"], row["wavelength"], row["fwhm"] ) - return new_name, new_value - # Go over all features. - for row in self.features: - 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 + elif kind == "attenuation": + fitter.register_attenuation(name, row["tau"]) + + elif kind == "absorption": + fitter.register_absorption( + name, row["tau"], row["wavelength"], row["fwhm"] + ) - # 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] 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!" + ) + + fitter.finalize_model() + return fitter @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 ee0e39e..cc7c6b9 100644 --- a/pahfit/tests/test_feature_parsing.py +++ b/pahfit/tests/test_feature_parsing.py @@ -38,8 +38,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) + fitter = m._set_up_fitter(instrumentname, 0) + m._ingest_fit_result_to_features(fitter) # 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 d56975b..56e4aa5 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 fce3acc..c15a429 100644 --- a/pahfit/tests/test_model_impl.py +++ b/pahfit/tests/test_model_impl.py @@ -35,14 +35,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( + + # test only works for the astropy-based implementation at the moment. + fit_result = model.fitter + reconstructed_fit_result = 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 + for name in fit_result.components(): + 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(): @@ -63,13 +66,13 @@ def test_model_edit(): assert model.features.loc[feature][col][0] == originalT # construct astropy model with dummy instrument - astropy_model_edit = model_to_edit._construct_astropy_model( + fitter_edit = model_to_edit._construct_model( instrumentname="spitzer.irs.*", redshift=0 ) - # 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 + # Make sure the change is reflected in this astropy model. Very + # handy that we can access the right component by the feature name! + assert fitter_edit.get_result(feature)["temperature"] == newT def test_model_tabulate(): From 87a6e184955899119d3a8eacfb6c3d14ec6502d8 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Mon, 13 May 2024 20:06:07 -0400 Subject: [PATCH 03/27] Regular floats for Features _bounds_dtype to exactly match test --- pahfit/features/features.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 4bb56f7..31da313 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): From a5e370d061a708ef0b353cf107b5683befefd7ea Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Mon, 13 May 2024 20:29:38 -0400 Subject: [PATCH 04/27] Fixes based on test cases --- pahfit/model.py | 15 +++++---------- pahfit/tests/test_model_impl.py | 16 ++++++++++------ 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index 5b77204..a5a1a76 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -318,7 +318,9 @@ def _convert_spec_data(spec, z): """ if not spec.flux.unit.is_equivalent(units.intensity): - raise PAHFITModelError("For now, PAHFIT only supports intensity units, i.e. convertible to MJy / sr.") + 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 @@ -418,7 +420,7 @@ def _ingest_fit_result_to_features(self, fitter: Fitter): """ # iterate over the list stored in fitter, so we only get - # components that were set up by _construct_model. Having an + # 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). @@ -433,11 +435,6 @@ def _ingest_fit_result_to_features(self, fitter: Fitter): self.features[column]["val"][i] = value else: self.features[column][i] = (value, np.nan, np.nan) - print(f'ingested {name} {column} as {value:6e}') - print('value in table ', self.features[column][i]) - if not self.features[column]["val"][i] == value: - print("assigment went wrong") - raise PAHFITModelError except Exception as e: print( f"Could not assign to name {name} in features table. Some diagnostic output below" @@ -800,9 +797,7 @@ def tabulate( flux_quantity = flux_values * u.dimensionless_unscaled else: user_unit = self.features.meta["user_unit"]["flux"] - flux_quantity = (flux_values * units.internal_flux_unit(user_unit)).to( - user_unit - ) + flux_quantity = (flux_values * units.intensity).to(user_unit) return Spectrum1D(spectral_axis=wav, flux=flux_quantity) diff --git a/pahfit/tests/test_model_impl.py b/pahfit/tests/test_model_impl.py index c15a429..72a1aca 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): @@ -60,13 +61,16 @@ 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 - fitter_edit = model_to_edit._construct_model( + fitter_edit = model_to_edit._set_up_fitter( instrumentname="spitzer.irs.*", redshift=0 ) From 9ded5723395896eaaaae3e4f1805f4a1baa31238 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Mon, 13 May 2024 20:41:44 -0400 Subject: [PATCH 05/27] Remove obsolete helpers.calculate_compounds, more formatting --- pahfit/features/features.py | 2 +- pahfit/helpers.py | 114 +++--------------------------------- 2 files changed, 9 insertions(+), 107 deletions(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 31da313..0efed64 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", float), ("min", float), ("max", float)]) # 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/helpers.py b/pahfit/helpers.py index fc84ee5..f6bcabc 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 From c1a2665c58cb49677ec1ff2ab3b468bf8ecd0657 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Mon, 13 May 2024 20:43:30 -0400 Subject: [PATCH 06/27] Remove mentions of PAHFITBase --- docs/index.rst | 4 ---- pahfit/tests/test_feature_parsing.py | 4 +++- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/docs/index.rst b/docs/index.rst index 21352c0..905cf99 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -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 diff --git a/pahfit/tests/test_feature_parsing.py b/pahfit/tests/test_feature_parsing.py index cc7c6b9..9791e4d 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. From 96757a664fd47b54130458948f5b5b8ac604948b Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 14 May 2024 11:23:02 -0400 Subject: [PATCH 07/27] Clarify Fitter API operation description --- pahfit/fitter.py | 36 +++++++++++++++++++++--------------- 1 file changed, 21 insertions(+), 15 deletions(-) diff --git a/pahfit/fitter.py b/pahfit/fitter.py index 5c8b657..f009da5 100644 --- a/pahfit/fitter.py +++ b/pahfit/fitter.py @@ -15,7 +15,7 @@ class Fitter(ABC): 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.) + 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 @@ -23,20 +23,26 @@ class Fitter(ABC): 4. Access to the evaluation of the underlying model (again with no assumptions like in step 2.). - For the model setup, multiple functions are used, so a few notes are - provided here. There is one function per type of component supported - by PAHFIT, and the arguments of these functions will ask for - different "standard" PAHFIT numbers, i.e. those from the Features - table. These functions have the same signature between all Fitter - implementations, so that the Model class can use a single - implementation to set up the Fitter. The Model has access to the - Features table and the instrument model, and needs to set up Fitter - with the correct initial values, bounds, and "fixed" flags (e.g. - setting a fixed FWHM based on the instrument for the lines). After - all the components have been added, the finalize_model() function - can be called to finish setting up the internal astropy model. After - this has finished, fit() can be called to apply the model and the - astropy fitter to the data. + 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. + register_line(). Once all components have been added, the + finalize_model() 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(). """ From a9311cc0ef6a401e07f4bf07ba5475cc6aa34357 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Fri, 17 May 2024 16:07:28 -0400 Subject: [PATCH 08/27] Array vs scalar arguments to determine fitted/fixed in Fitter Clean up Features contents before giving them to Fitter. Nan is converted to np.inf or -np.inf. Fixed values (both bounds nan) are converted to a single scalar. The feature adding functions of Fitter will interpret scalars as fixed parameters, and 3-arrays as the initial value and bounds of parameters that will be fit. --- pahfit/apfitter.py | 34 +++++++++++++++++----------------- pahfit/fitter.py | 9 ++++----- pahfit/model.py | 44 ++++++++++++++++++++++++++++++++++---------- 3 files changed, 55 insertions(+), 32 deletions(-) diff --git a/pahfit/apfitter.py b/pahfit/apfitter.py index 4a3d051..37f1412 100644 --- a/pahfit/apfitter.py +++ b/pahfit/apfitter.py @@ -180,7 +180,7 @@ def register_line(self, name, power, wavelength, fwhm): kwargs = self._astropy_model_kwargs( name, ["amplitude", "mean", "stddev"], - [power, wavelength, np.array([x for x in fwhm]) / 2.355], + [power, wavelength, fwhm / 2.355], ) self._register_component(Gaussian1D, **kwargs) @@ -357,7 +357,7 @@ def get_result(self, component_name): raise PAHFITModelError(f"Unsupported component type: {c_type}") @staticmethod - def _astropy_model_kwargs(component_name, param_names, value_tuples): + 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 @@ -389,11 +389,10 @@ def _astropy_model_kwargs(component_name, param_names, value_tuples): Names of the parameters for the astropy model, e.g. ["dust_feature1", "dust_feature2"] - value_tuples : list of (array of size 3) + 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], i.e. in the format as stored in the - Features table. This means that [value, masked, masked] will - result in a fixed parameter. + min_bound, max_bound] for variable parameters, or a scalar + (single float) for fixed parameters. Returns ------- @@ -403,18 +402,19 @@ def _astropy_model_kwargs(component_name, param_names, value_tuples): # basic format of constructor parameters of astropy model kwargs = {"name": component_name, "bounds": {}, "fixed": {}} - for param_name, value_tuple in zip(param_names, value_tuples): - kwargs[param_name] = value_tuple[0] - - # astropy does not like numpy bools, so we do this silly - # conversion. - is_fixed = np.isnan(value_tuple[1]) and np.isnan(value_tuple[2]) - kwargs["fixed"][param_name] = True if is_fixed else False + 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, and None if variable but unbounded. - min_max = value_tuple[1], value_tuple[2] - limits = [0 if is_fixed else v for v in min_max] - kwargs["bounds"][param_name] = [None if np.isinf(x) else x for x in limits] + # 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/fitter.py b/pahfit/fitter.py index f009da5..877b09d 100644 --- a/pahfit/fitter.py +++ b/pahfit/fitter.py @@ -90,11 +90,10 @@ def register_starlight(self, name, temperature, tau): Unique name. Will be used to allow retrieval of the results after the fitting. - temperature : array of size 3 - Blackbody temperature, given as [value, lower_bound, - upper_bound]. The bounds assume the same system as the - features table: if they are masked, the parameter will be - fixed, while +inf or -inf mean unbounded. + 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. diff --git a/pahfit/model.py b/pahfit/model.py index a5a1a76..9e77e65 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -864,18 +864,28 @@ def _set_up_fitter( excluded = self._excluded_features(instrumentname, redshift, x) - def array3(features_tuple3): - return np.array([x for x in features_tuple3]) + def cleaned(features_tuple3): + val = features_tuple3[0] + if bounded_is_fixed(features_tuple3): + return val + else: + 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]) for row in self.features[~excluded]: kind = row["kind"] name = row["name"] if kind == "starlight": - fitter.register_starlight(name, row["temperature"], row["tau"]) + fitter.register_starlight( + name, cleaned(row["temperature"]), cleaned(row["tau"]) + ) elif kind == "dust_continuum": - fitter.register_dust_continuum(name, row["temperature"], row["tau"]) + fitter.register_dust_continuum( + name, cleaned(row["temperature"]), cleaned(row["tau"]) + ) elif kind == "line": if use_instrument_fwhm: @@ -892,22 +902,36 @@ def array3(features_tuple3): fwhm = instrument.fwhm(instrumentname, w_obs, as_bounded=True)[ 0 ] / (1.0 + redshift) - fwhm = np.ma.filled(fwhm, np.nan) + + # decide if scalar (fixed) or tuple (fitted fwhm + # between upper and lower fwhm limits, happens in + # case of overlapping instruments) + if np.ma.is_masked(fwhm): + fwhm = fwhm[0] else: - fwhm = row["fwhm"] - fitter.register_line(name, row["power"], row["wavelength"], fwhm) + fwhm = cleaned(row["fwhm"]) + + fitter.register_line( + name, cleaned(row["power"]), cleaned(row["wavelength"]), fwhm + ) elif kind == "dust_feature": fitter.register_dust_feature( - name, row["power"], row["wavelength"], row["fwhm"] + name, + cleaned(row["power"]), + cleaned(row["wavelength"]), + cleaned(row["fwhm"]), ) elif kind == "attenuation": - fitter.register_attenuation(name, row["tau"]) + fitter.register_attenuation(name, cleaned(row["tau"])) elif kind == "absorption": fitter.register_absorption( - name, row["tau"], row["wavelength"], row["fwhm"] + name, + cleaned(row["tau"]), + cleaned(row["wavelength"]), + cleaned(row["fwhm"]), ) else: From 02cd40bdaa172a3902b8e5d797473f094b32beb7 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Fri, 17 May 2024 16:15:29 -0400 Subject: [PATCH 09/27] Remove clear() from Fitter and APFitter --- pahfit/apfitter.py | 8 -------- pahfit/fitter.py | 10 ---------- 2 files changed, 18 deletions(-) diff --git a/pahfit/apfitter.py b/pahfit/apfitter.py index 37f1412..0ce006c 100644 --- a/pahfit/apfitter.py +++ b/pahfit/apfitter.py @@ -63,14 +63,6 @@ def __init__(self): After construction, use the register_() functions to start setting up a model, then call finalize_model(). - """ - self.clear() - - def clear(self): - """Reset model. - - After reset, register_() and finalize_model() can be used again. - """ self.additive_components = [] self.multiplicative_components = [] diff --git a/pahfit/fitter.py b/pahfit/fitter.py index 877b09d..3ec5052 100644 --- a/pahfit/fitter.py +++ b/pahfit/fitter.py @@ -45,16 +45,6 @@ class Fitter(ABC): passing the component name to get_result(). """ - - @abstractmethod - def clear(self): - """Reset model. - - After reset, register_() and finalize_model() can be used again. - - """ - pass - @abstractmethod def components(self): """Return list of features. From 45d14975eababec1221af88fa5dad2e74661a726 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Fri, 17 May 2024 16:38:34 -0400 Subject: [PATCH 10/27] Rename register_* to add_feature_*, clean up APFitter var names Also worked on more consistent naming of variables in APFitter implementation --- pahfit/apfitter.py | 66 ++++++++++++++++++++++++---------------------- pahfit/fitter.py | 23 ++++++++-------- pahfit/model.py | 14 +++++----- 3 files changed, 54 insertions(+), 49 deletions(-) diff --git a/pahfit/apfitter.py b/pahfit/apfitter.py index 0ce006c..0e5035e 100644 --- a/pahfit/apfitter.py +++ b/pahfit/apfitter.py @@ -60,20 +60,20 @@ class APFitter(Fitter): def __init__(self): """Construct a new fitter. - After construction, use the register_() functions to start - setting up a model, then call finalize_model(). + After construction, use the add_feature_() functions to start + setting up a model, then call finalize(). """ self.additive_components = [] self.multiplicative_components = [] - self.component_types = {} + self.feature_types = {} self.model = None self.message = None def components(self): """Return list of component names. - Only works after finalize_model(). + Only works after finalize(). """ if hasattr(self.model, "submodel_names"): @@ -82,10 +82,10 @@ def components(self): # single-component edge case return [self.model.name] - def finalize_model(self): + def finalize(self): """Sum the registered components into one CompoundModel. - To be called after a series of "register_()" calls, and before + To be called after a series of "add_feature_()" calls, and before using the other functionality (fitting and evaluating). """ @@ -99,8 +99,10 @@ def finalize_model(self): for c in self.multiplicative_components: self.model *= c - def _register_component(self, astropy_model_class, multiplicative=False, **kwargs): - """Register any component in a generic way. + 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 ---------- @@ -115,7 +117,7 @@ def _register_component(self, astropy_model_class, multiplicative=False, **kwarg kwargs : dict Arguments for the astropy model constructor, including a unique value for "name". Should be generated with - self._astropy_model_kwargs; the register_() functions show how + self._astropy_model_kwargs; the add_feature_() functions show how to do this for each type of feature. """ @@ -124,7 +126,7 @@ def _register_component(self, astropy_model_class, multiplicative=False, **kwarg else: self.additive_components.append(astropy_model_class(**kwargs)) - def register_starlight(self, name, temperature, tau): + def add_feature_starlight(self, name, temperature, tau): """Register a BlackBody1D. Parameters @@ -141,76 +143,76 @@ def register_starlight(self, name, temperature, tau): tau : analogous. Used as amplitude. """ - self.component_types[name] = "starlight" + self.feature_types[name] = "starlight" kwargs = self._astropy_model_kwargs( name, ["temperature", "amplitude"], [temperature, tau] ) - self._register_component(BlackBody1D, **kwargs) + self._add_component(BlackBody1D, **kwargs) - def register_dust_continuum(self, name, temperature, tau): + def add_feature_dust_continuum(self, name, temperature, tau): """Register a ModifiedBlackBody1D. Analogous. Temperature and tau are used as temperature and amplitude """ - self.component_types[name] = "dust_continuum" + self.feature_types[name] = "dust_continuum" kwargs = self._astropy_model_kwargs( name, ["temperature", "amplitude"], [temperature, tau] ) - self._register_component(ModifiedBlackBody1D, **kwargs) + self._add_component(ModifiedBlackBody1D, **kwargs) - def register_line(self, name, power, wavelength, fwhm): + 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.component_types[name] = "line" + self.feature_types[name] = "line" kwargs = self._astropy_model_kwargs( name, ["amplitude", "mean", "stddev"], - [power, wavelength, fwhm / 2.355], + [power, wavelength, fwhm / 2.355], ) - self._register_component(Gaussian1D, **kwargs) + self._add_component(Gaussian1D, **kwargs) - def register_dust_feature(self, name, power, wavelength, fwhm): + 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.component_types[name] = "dust_feature" + self.feature_types[name] = "dust_feature" kwargs = self._astropy_model_kwargs( name, ["amplitude", "x_0", "fwhm"], [power, wavelength, fwhm] ) - self._register_component(Drude1D, **kwargs) + self._add_component(Drude1D, **kwargs) - def register_attenuation(self, name, tau): + def add_feature_attenuation(self, name, tau): """Register the S07 attenuation component. Analogous. Uses tau as tau_sil for S07_attenuation. Is multiplicative. """ - self.component_types[name] = "attenuation" + self.feature_types[name] = "attenuation" kwargs = self._astropy_model_kwargs(name, ["tau_sil"], [tau]) - self._register_component(S07_attenuation, multiplicative=True, **kwargs) + self._add_component(S07_attenuation, multiplicative=True, **kwargs) - def register_absorption(self, name, tau, wavelength, fwhm): + def add_feature_absorption(self, name, tau, wavelength, fwhm): """Register an absorbing Drude1D component. Analogous. Is multiplicative. """ - self.component_types[name] = "absorption" + self.feature_types[name] = "absorption" kwargs = self._astropy_model_kwargs( name, ["tau", "x_0", "fwhm"], [tau, wavelength, fwhm] ) - self._register_component(att_Drude1D, multiplicative=True, **kwargs) + self._add_component(att_Drude1D, multiplicative=True, **kwargs) def evaluate_model(self, xz): """Evaluate internal astropy model with its current parameters. @@ -298,7 +300,7 @@ def get_result(self, component_name): Parameters ---------- component_name : str - One of the names provided to any of the register_*() calls + One of the names provided to any of the add_feature_*() calls made during setup. See also Fitter.components(). Returns @@ -319,7 +321,7 @@ def get_result(self, component_name): # CompoundModel but normal single-component model. component = self.model - c_type = self.component_types[component_name] + c_type = self.feature_types[component_name] if c_type == "starlight" or c_type == "dust_continuum": return { "temperature": component.temperature.value, @@ -407,6 +409,8 @@ def _astropy_model_kwargs(component_name, param_names, param_values): # 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)] + kwargs["bounds"][param_name] = [ + None if np.isinf(x) else x for x in (lo_bound, up_bound) + ] return kwargs diff --git a/pahfit/fitter.py b/pahfit/fitter.py index 3ec5052..f81b959 100644 --- a/pahfit/fitter.py +++ b/pahfit/fitter.py @@ -36,8 +36,8 @@ class Fitter(ABC): During the Fitter setup, the initial values, bounds, and "fixed" flags are passed using one function call for each component, e.g. - register_line(). Once all components have been added, the - finalize_model() function should be called; some subclasses (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 @@ -45,18 +45,19 @@ class Fitter(ABC): passing the component name to get_result(). """ + @abstractmethod def components(self): """Return list of features. - Only works after finalize_model(). Will return the names passed + Only works after finalize(). Will return the names passed using the register functions. """ pass @abstractmethod - def finalize_model(self): + def finalize(self): """Process the registered features and prepare for fitting. The register functions below allow adding individual features. @@ -68,7 +69,7 @@ def finalize_model(self): pass @abstractmethod - def register_starlight(self, name, temperature, tau): + def add_feature_starlight(self, name, temperature, tau): """Register a starlight feature. The exact representation depends on the implementation, but the @@ -92,12 +93,12 @@ def register_starlight(self, name, temperature, tau): pass @abstractmethod - def register_dust_continuum(self, name, temperature, tau): + def add_feature_dust_continuum(self, name, temperature, tau): """Register a dust continuum feature.""" pass @abstractmethod - def register_line(self, name, power, wavelength, fwhm): + def add_feature_line(self, name, power, wavelength, fwhm): """Register an emission line feature. Typically a Gaussian profile. @@ -106,7 +107,7 @@ def register_line(self, name, power, wavelength, fwhm): pass @abstractmethod - def register_dust_feature(self, name, power, wavelength, fwhm): + def add_feature_dust_feature(self, name, power, wavelength, fwhm): """Register a dust feature. Typically a Drude profile. @@ -115,7 +116,7 @@ def register_dust_feature(self, name, power, wavelength, fwhm): pass @abstractmethod - def register_attenuation(self, name, tau): + def add_feature_attenuation(self, name, tau): """Register the S07 attenuation component. Other types of attenuation might be possible in the future. Is @@ -125,7 +126,7 @@ def register_attenuation(self, name, tau): pass @abstractmethod - def register_absorption(self, name, tau, wavelength, fwhm): + def add_feature_absorption(self, name, tau, wavelength, fwhm): """Register an absorption feature. Typically a Drude profile. Is multiplicative. @@ -182,7 +183,7 @@ def get_result(self, feature_name): Parameters ---------- component_name : str - One of the names provided to any of the register_() calls + One of the names provided to any of the add_feature_() calls made during setup. See also Fitter.components(). Returns diff --git a/pahfit/model.py b/pahfit/model.py index 9e77e65..3a2493a 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -878,12 +878,12 @@ def cleaned(features_tuple3): name = row["name"] if kind == "starlight": - fitter.register_starlight( + fitter.add_feature_starlight( name, cleaned(row["temperature"]), cleaned(row["tau"]) ) elif kind == "dust_continuum": - fitter.register_dust_continuum( + fitter.add_feature_dust_continuum( name, cleaned(row["temperature"]), cleaned(row["tau"]) ) @@ -911,12 +911,12 @@ def cleaned(features_tuple3): else: fwhm = cleaned(row["fwhm"]) - fitter.register_line( + fitter.add_feature_line( name, cleaned(row["power"]), cleaned(row["wavelength"]), fwhm ) elif kind == "dust_feature": - fitter.register_dust_feature( + fitter.add_feature_dust_feature( name, cleaned(row["power"]), cleaned(row["wavelength"]), @@ -924,10 +924,10 @@ def cleaned(features_tuple3): ) elif kind == "attenuation": - fitter.register_attenuation(name, cleaned(row["tau"])) + fitter.add_feature_attenuation(name, cleaned(row["tau"])) elif kind == "absorption": - fitter.register_absorption( + fitter.add_feature_absorption( name, cleaned(row["tau"]), cleaned(row["wavelength"]), @@ -939,7 +939,7 @@ def cleaned(features_tuple3): f"Model components of kind {kind} are not implemented!" ) - fitter.finalize_model() + fitter.finalize() return fitter @staticmethod From add53abdb0c7d4e9b2937a39fcb3b348836ab93d Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Mon, 20 May 2024 15:57:13 -0400 Subject: [PATCH 11/27] Model and geometry keywords for add_feature_attenuation/absorption --- pahfit/apfitter.py | 9 +++++++-- pahfit/fitter.py | 6 +++--- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/pahfit/apfitter.py b/pahfit/apfitter.py index 0e5035e..22b1db0 100644 --- a/pahfit/apfitter.py +++ b/pahfit/apfitter.py @@ -191,18 +191,23 @@ def add_feature_dust_feature(self, name, power, wavelength, fwhm): ) self._add_component(Drude1D, **kwargs) - def add_feature_attenuation(self, name, tau): + 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): + def add_feature_absorption(self, name, tau, wavelength, fwhm, geometry='screen'): """Register an absorbing Drude1D component. Analogous. Is multiplicative. diff --git a/pahfit/fitter.py b/pahfit/fitter.py index f81b959..ad474a6 100644 --- a/pahfit/fitter.py +++ b/pahfit/fitter.py @@ -116,7 +116,7 @@ def add_feature_dust_feature(self, name, power, wavelength, fwhm): pass @abstractmethod - def add_feature_attenuation(self, name, tau): + 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 @@ -126,10 +126,10 @@ def add_feature_attenuation(self, name, tau): pass @abstractmethod - def add_feature_absorption(self, name, tau, wavelength, fwhm): + def add_feature_absorption(self, name, tau, wavelength, fwhm, geometry='screen'): """Register an absorption feature. - Typically a Drude profile. Is multiplicative. + Modeled by a Drude profile. Is multiplicative. """ pass From 7aed8516ee6285fffb3351e46b239ffb1819b957 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Mon, 20 May 2024 16:51:24 -0400 Subject: [PATCH 12/27] Replace xz, yz by lam, flux --- pahfit/apfitter.py | 32 +++++++-------- pahfit/fitter.py | 18 ++++----- pahfit/model.py | 98 ++++++++++++++++++++++++---------------------- 3 files changed, 76 insertions(+), 72 deletions(-) diff --git a/pahfit/apfitter.py b/pahfit/apfitter.py index 22b1db0..14a3468 100644 --- a/pahfit/apfitter.py +++ b/pahfit/apfitter.py @@ -219,22 +219,22 @@ def add_feature_absorption(self, name, tau, wavelength, fwhm, geometry='screen') ) self._add_component(att_Drude1D, multiplicative=True, **kwargs) - def evaluate_model(self, xz): + def evaluate(self, lam): """Evaluate internal astropy model with its current parameters. Parameters ---------- - xz : array + lam : array Rest frame wavelengths in micron Returns ------- - yz : array + flux : array Rest frame flux in internal units """ - return self.model(xz) + return self.model(lam) - def fit(self, xz, yz, uncz, maxiter=10000): + 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 @@ -249,37 +249,37 @@ def fit(self, xz, yz, uncz, maxiter=10000): Retrieval of uncertainties and fit details is yet to be implemented. - CAVEAT: flux unit (yz) is still ambiguous, since it can be flux - density or intensity, according to the options defined in + 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 wavelenght unit). + times wavelength unit). Parameters ---------- - xz : array + lam : array Rest frame wavelengths in micron - yz : array + flux : array Rest frame flux in internal units. - uncz : array - Uncertainty on rest frame flux. Same units as yz. + unc : array + Uncertainty on rest frame flux. Same units as flux. """ # clean, because astropy does not like nan - w = 1 / uncz + w = 1 / unc # make sure there are no zero uncertainties either - mask = np.isfinite(xz) & np.isfinite(yz) & np.isfinite(w) + mask = np.isfinite(lam) & np.isfinite(flux) & np.isfinite(w) self.fit_info = [] fit = LevMarLSQFitter(calc_uncertainties=True) astropy_result = fit( self.model, - xz[mask], - yz[mask], + lam[mask], + flux[mask], weights=w[mask], maxiter=maxiter, epsilon=1e-10, diff --git a/pahfit/fitter.py b/pahfit/fitter.py index ad474a6..df4f193 100644 --- a/pahfit/fitter.py +++ b/pahfit/fitter.py @@ -135,24 +135,24 @@ def add_feature_absorption(self, name, tau, wavelength, fwhm, geometry='screen') pass @abstractmethod - def evaluate_model(self, xz): - """Evaluate the model at the given wavelengths. + def evaluate(self, lam): + """Evaluate the fitting function at the given wavelengths. Parameters ---------- - xz : array + lam : array Rest frame wavelengths in micron Returns ------- - yz : array + flux : array Rest frame flux in internal units """ pass @abstractmethod - def fit(self, xz, yz, uncz, maxiter=1000): + 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 @@ -164,14 +164,14 @@ def fit(self, xz, yz, uncz, maxiter=1000): Parameters ---------- - xz : array + lam : array Rest frame wavelengths in micron - yz : array + flux : array Rest frame flux in internal units. - uncz : array - Uncertainty on rest frame flux. Same units as yz. + unc : array + Uncertainty on rest frame flux. Same units as flux. """ pass diff --git a/pahfit/model.py b/pahfit/model.py index 3a2493a..004c3e7 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -187,12 +187,12 @@ def guess( # 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): @@ -235,7 +235,7 @@ def dust_continuum_guess(row): else: w_ref = wmin - flux_ref = np.median(yz[(xz > w_ref - 0.2) & (xz < w_ref + 0.2)]) + flux_ref = np.median(flux[(lam > w_ref - 0.2) & (lam < w_ref + 0.2)]) amp_guess = flux_ref / bb(w_ref) return amp_guess / nbb @@ -257,12 +257,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 +273,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: @@ -313,7 +313,7 @@ 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 """ @@ -321,15 +321,15 @@ def _convert_spec_data(spec, z): 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 + 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, @@ -390,7 +390,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 @@ -402,7 +402,7 @@ def fit( self.fitter = self._set_up_fitter( inst, z, x=x, use_instrument_fwhm=use_instrument_fwhm ) - self.fitter.fit(xz, yz, uncz, maxiter=maxiter) + self.fitter.fit(lam, flux, unc, maxiter=maxiter) # copy the fit results to the features table self._ingest_fit_result_to_features(self.fitter) @@ -485,9 +485,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, @@ -515,7 +515,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( @@ -527,11 +527,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 = [ @@ -547,32 +547,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, @@ -581,7 +583,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, @@ -590,7 +592,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], @@ -601,7 +603,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( @@ -614,7 +616,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}$") @@ -637,7 +639,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] @@ -658,7 +660,7 @@ def tabulate_components(kind): ax.axhline(0, linestyle="--", color="gray", zorder=0) ax.plot( - xz, + lam, res, "ko", fillstyle="none", @@ -668,7 +670,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 [%]") @@ -801,13 +803,13 @@ def tabulate( return Spectrum1D(spectral_axis=wav, flux=flux_quantity) - def _excluded_features(self, instrumentname, redshift, x=None): + def _excluded_features(self, instrumentname, redshift, lam_obs=None): """Determine excluded features Based on instrument wavelength range. instrumentname : str Qualified instrument name - x : array + lam_obs : array Optional observed wavelength grid. Exclude any lines and dust features outside of this range. @@ -816,14 +818,16 @@ def _excluded_features(self, instrumentname, redshift, x=None): array of bool, same length as self.features, True where features are far outside the wavelength range. """ - observed_wavs = self.features["wavelength"]["val"] * (1 + redshift) + lam_feature_obs = self.features["wavelength"]["val"] * (1 + redshift) # has wavelength and not within instrument range - is_outside = ~instrument.within_segment(observed_wavs, instrumentname) + is_outside = ~instrument.within_segment(lam_feature_obs, instrumentname) # also apply observed range if provided - if x is not None: - is_outside |= (observed_wavs < np.amin(x)) | (observed_wavs > np.amax(x)) + 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"] @@ -849,7 +853,7 @@ def _set_up_fitter( detail). - Features outside the appropriate wavelength range should not be added to the Fitter: the "trimming" is done here, using the - given wavelength range xz (optional). + given wavelength range lam (optional). TODO: flags to indicate which features were excluded. From edcbbf68a168adb69746272ab12ef7e954504850 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Mon, 20 May 2024 17:19:56 -0400 Subject: [PATCH 13/27] Docstring and formatting fixes --- pahfit/fitter.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pahfit/fitter.py b/pahfit/fitter.py index df4f193..e3eda0b 100644 --- a/pahfit/fitter.py +++ b/pahfit/fitter.py @@ -116,7 +116,7 @@ def add_feature_dust_feature(self, name, power, wavelength, fwhm): pass @abstractmethod - def add_feature_attenuation(self, name, tau, model='S07', geometry='screen'): + 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 @@ -126,7 +126,7 @@ def add_feature_attenuation(self, name, tau, model='S07', geometry='screen'): pass @abstractmethod - def add_feature_absorption(self, name, tau, wavelength, fwhm, geometry='screen'): + def add_feature_absorption(self, name, tau, wavelength, fwhm, geometry="screen"): """Register an absorption feature. Modeled by a Drude profile. Is multiplicative. @@ -193,7 +193,7 @@ def get_result(self, feature_name): 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'} + e.g. {'name': 'line0', 'power': value, 'fwhm': value, 'wavelength': value} """ pass From 4057012da193d0dcfa044d4149cf20dd8f4f7d39 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 21 May 2024 15:22:59 -0400 Subject: [PATCH 14/27] Add TODO in guess to replace fwhm guessing by sigma_v --- pahfit/model.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pahfit/model.py b/pahfit/model.py index 004c3e7..60a059d 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -285,7 +285,11 @@ def amp_guess(row, fwhm): loop_over_non_fixed("line", "power", lambda row: some_flux) # Set the fwhms in the features table. Slightly different logic, - # as the fwhm for lines are masked by default + # 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: for row_index in np.where(self.features["kind"] == "line")[0]: row = self.features[row_index] From 8f4d2bc78fdfab99856185493cb26f72188bd681 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 21 May 2024 17:38:08 -0400 Subject: [PATCH 15/27] Work with self.fitter directly instead fitter argument/return --- pahfit/model.py | 43 +++++++++++----------------- pahfit/tests/test_feature_parsing.py | 4 +-- 2 files changed, 19 insertions(+), 28 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index 60a059d..9df5084 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -403,24 +403,20 @@ def fit( # check if observed spectrum is compatible with instrument model instrument.check_range([min(x), max(x)], inst) - self.fitter = self._set_up_fitter( - inst, z, x=x, use_instrument_fwhm=use_instrument_fwhm - ) + self._set_up_fitter(inst, z, x=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(self.fitter) + self._ingest_fit_result_to_features() if verbose: print(self.fitter.message) - def _ingest_fit_result_to_features(self, fitter: Fitter): - """Copy the results from a Fitter to the features table + def _ingest_fit_result_to_features(self): + """Copy the results from the Fitter to the features table This is a utility method, executed only at the end of fit(), - where the Fitter is passed after Fitter.fit() has been applied. - Passing a different Fitter object could be useful for - testing. + where Fitter.fit() has been applied. """ # iterate over the list stored in fitter, so we only get @@ -428,10 +424,10 @@ def _ingest_fit_result_to_features(self, fitter: Fitter): # ENABLED/DISABLED flag for every feature would be a nice # alternative (and clear for the user). - self.features.meta["fitter_message"] = fitter.message + self.features.meta["fitter_message"] = self.fitter.message - for name in fitter.components(): - for column, value in fitter.get_result(name).items(): + for name in self.fitter.components(): + 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 @@ -844,7 +840,7 @@ def _excluded_features(self, instrumentname, redshift, lam_obs=None): def _set_up_fitter( self, instrumentname, redshift, x=None, use_instrument_fwhm=True ): - """Convert features table to Fitter instance. + """Convert features table to Fitter instance, set self.fitter. For every row of the features table, calls a function of Fitter API to register an appropriate component. Finalizes the Fitter @@ -861,14 +857,10 @@ def _set_up_fitter( TODO: flags to indicate which features were excluded. - Returns - ------- - Fitter - """ # Fitting implementation can be changed by choosing another # Fitter class. TODO: make this configurable. - fitter = APFitter() + self.fitter = APFitter() excluded = self._excluded_features(instrumentname, redshift, x) @@ -886,12 +878,12 @@ def cleaned(features_tuple3): name = row["name"] if kind == "starlight": - fitter.add_feature_starlight( + self.fitter.add_feature_starlight( name, cleaned(row["temperature"]), cleaned(row["tau"]) ) elif kind == "dust_continuum": - fitter.add_feature_dust_continuum( + self.fitter.add_feature_dust_continuum( name, cleaned(row["temperature"]), cleaned(row["tau"]) ) @@ -919,12 +911,12 @@ def cleaned(features_tuple3): else: fwhm = cleaned(row["fwhm"]) - fitter.add_feature_line( + self.fitter.add_feature_line( name, cleaned(row["power"]), cleaned(row["wavelength"]), fwhm ) elif kind == "dust_feature": - fitter.add_feature_dust_feature( + self.fitter.add_feature_dust_feature( name, cleaned(row["power"]), cleaned(row["wavelength"]), @@ -932,10 +924,10 @@ def cleaned(features_tuple3): ) elif kind == "attenuation": - fitter.add_feature_attenuation(name, cleaned(row["tau"])) + self.fitter.add_feature_attenuation(name, cleaned(row["tau"])) elif kind == "absorption": - fitter.add_feature_absorption( + self.fitter.add_feature_absorption( name, cleaned(row["tau"]), cleaned(row["wavelength"]), @@ -947,8 +939,7 @@ def cleaned(features_tuple3): f"Model components of kind {kind} are not implemented!" ) - fitter.finalize() - return fitter + self.fitter.finalize() @staticmethod def _parse_instrument_and_redshift(spec, redshift): diff --git a/pahfit/tests/test_feature_parsing.py b/pahfit/tests/test_feature_parsing.py index 9791e4d..d658a68 100644 --- a/pahfit/tests/test_feature_parsing.py +++ b/pahfit/tests/test_feature_parsing.py @@ -40,8 +40,8 @@ def test_feature_parsing(): def test_parsing(features_edit): m = Model(features_edit) - fitter = m._set_up_fitter(instrumentname, 0) - m._ingest_fit_result_to_features(fitter) + m._set_up_fitter(instrumentname, 0) + m._ingest_fit_result_to_features() # Case 0: the whole table test_parsing(features) From ac2146fb95480a7c7690944b738cd33f7c18595c Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 21 May 2024 17:52:26 -0400 Subject: [PATCH 16/27] Remove Fitter.components(), instead keep 'enabled' list in Model --- pahfit/apfitter.py | 16 ++-------------- pahfit/fitter.py | 12 +----------- pahfit/model.py | 11 +++++------ pahfit/tests/test_model_impl.py | 13 +++++++------ 4 files changed, 15 insertions(+), 37 deletions(-) diff --git a/pahfit/apfitter.py b/pahfit/apfitter.py index 14a3468..9f11f33 100644 --- a/pahfit/apfitter.py +++ b/pahfit/apfitter.py @@ -70,18 +70,6 @@ def __init__(self): self.model = None self.message = None - def components(self): - """Return list of component names. - - Only works after finalize(). - - """ - if hasattr(self.model, "submodel_names"): - return self.model.submodel_names - else: - # single-component edge case - return [self.model.name] - def finalize(self): """Sum the registered components into one CompoundModel. @@ -207,7 +195,7 @@ def add_feature_attenuation(self, name, tau, model="S07", geometry="screen"): 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'): + def add_feature_absorption(self, name, tau, wavelength, fwhm, geometry="screen"): """Register an absorbing Drude1D component. Analogous. Is multiplicative. @@ -306,7 +294,7 @@ def get_result(self, component_name): ---------- component_name : str One of the names provided to any of the add_feature_*() calls - made during setup. See also Fitter.components(). + made during setup. Returns ------- diff --git a/pahfit/fitter.py b/pahfit/fitter.py index e3eda0b..a75f2bd 100644 --- a/pahfit/fitter.py +++ b/pahfit/fitter.py @@ -46,16 +46,6 @@ class Fitter(ABC): """ - @abstractmethod - def components(self): - """Return list of features. - - Only works after finalize(). Will return the names passed - using the register functions. - - """ - pass - @abstractmethod def finalize(self): """Process the registered features and prepare for fitting. @@ -184,7 +174,7 @@ def get_result(self, feature_name): ---------- component_name : str One of the names provided to any of the add_feature_() calls - made during setup. See also Fitter.components(). + made during setup. Returns ------- diff --git a/pahfit/model.py b/pahfit/model.py index 9df5084..bbe46d2 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -12,7 +12,6 @@ from pahfit import instrument from pahfit.errors import PAHFITModelError from pahfit.component_models import BlackBody1D, S07_attenuation -from pahfit.fitter import Fitter from pahfit.apfitter import APFitter @@ -426,7 +425,7 @@ def _ingest_fit_result_to_features(self): self.features.meta["fitter_message"] = self.fitter.message - for name in self.fitter.components(): + 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] @@ -781,9 +780,8 @@ def tabulate( # need to wrap in try block to avoid bug: if all components are # removed (because of wavelength range considerations), it won't work try: - fitter = alt_model._set_up_fitter( - instrumentname, z, use_instrument_fwhm=False - ) + alt_model._set_up_fitter(instrumentname, z, use_instrument_fwhm=False) + fitter = alt_model.fitter except PAHFITModelError: return Spectrum1D( spectral_axis=wav, flux=np.zeros(wav.shape) * u.dimensionless_unscaled @@ -791,7 +789,7 @@ def tabulate( # shift the "observed wavelength grid" to "physical wavelength grid wav /= 1 + z - flux_values = fitter.evaluate_model(wav.value) + flux_values = fitter.evaluate(wav.value) # apply unit stored in features table (comes from from last fit # or from loading previous result from disk) @@ -863,6 +861,7 @@ def _set_up_fitter( self.fitter = APFitter() excluded = self._excluded_features(instrumentname, redshift, x) + self.enabled_features = self.features["name"][~excluded] def cleaned(features_tuple3): val = features_tuple3[0] diff --git a/pahfit/tests/test_model_impl.py b/pahfit/tests/test_model_impl.py index 72a1aca..6007c17 100644 --- a/pahfit/tests/test_model_impl.py +++ b/pahfit/tests/test_model_impl.py @@ -37,12 +37,12 @@ def test_feature_table_model_conversion(): # correct, reconstructing the model from model.features should # result in the exact same model. - # test only works for the astropy-based implementation at the moment. fit_result = model.fitter - reconstructed_fit_result = model._set_up_fitter( + model._set_up_fitter( instrumentname=spec.meta["instrument"], redshift=0, use_instrument_fwhm=False ) - for name in fit_result.components(): + 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: @@ -70,12 +70,13 @@ def test_model_edit(): assert model.features[col]["val"][j] == originalT # construct astropy model with dummy instrument - fitter_edit = model_to_edit._set_up_fitter( + 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 astropy model. Very - # handy that we can access the right component by the feature name! + # Make sure the change is reflected in this model. Very handy that + # we can access the right component by the feature name! assert fitter_edit.get_result(feature)["temperature"] == newT From 2942c08ddbed353151648eb16277c61fcd6a11f5 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 25 Jun 2024 09:59:10 -0400 Subject: [PATCH 17/27] Don't copy when applying feature mask in tabulate --- pahfit/model.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index bbe46d2..034c9d3 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -761,8 +761,7 @@ def tabulate( # 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 From 30b0c2eb8320c809e3f5b7e6fc4d59a8619f9daa Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 25 Jun 2024 10:14:38 -0400 Subject: [PATCH 18/27] Better documentation and variable names for _set_up_fitter --- pahfit/model.py | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index 034c9d3..212e223 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -402,7 +402,7 @@ def fit( # check if observed spectrum is compatible with instrument model instrument.check_range([min(x), max(x)], inst) - self._set_up_fitter(inst, z, x=x, use_instrument_fwhm=use_instrument_fwhm) + 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 @@ -834,8 +834,8 @@ def _excluded_features(self, instrumentname, redshift, lam_obs=None): return is_outside & is_excludable - def _set_up_fitter( - self, instrumentname, redshift, x=None, use_instrument_fwhm=True + def _set_up_fitter( + self, instrumentname, redshift, lam=None, use_instrument_fwhm=True ): """Convert features table to Fitter instance, set self.fitter. @@ -854,12 +854,29 @@ def _set_up_fitter( TODO: flags to indicate which features were excluded. + Parameters + ---------- + + instrumentname and redshift : needed to apply the instrument + model and to determine which feature to exclude + + lam : array of observed wavelengths + Used to exclude features from the model based on the actual + observed data given. + + 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). + """ # Fitting implementation can be changed by choosing another # Fitter class. TODO: make this configurable. self.fitter = APFitter() - excluded = self._excluded_features(instrumentname, redshift, x) + excluded = self._excluded_features(instrumentname, redshift, lam) self.enabled_features = self.features["name"][~excluded] def cleaned(features_tuple3): From e2410772a819aaf26efded3d84a3aa4fb9b40e0d Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 25 Jun 2024 10:24:15 -0400 Subject: [PATCH 19/27] Replace wavelengths named w by lam --- pahfit/model.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index 212e223..53c8518 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -228,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(flux[(lam > w_ref - 0.2) & (lam < 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) @@ -909,12 +909,12 @@ def cleaned(features_tuple3): # oberved wav; 3. shift back to rest frame wav # (width in rest frame will be narrower than # observed width) - w_obs = row["wavelength"]["val"] * (1.0 + redshift) + 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) - fwhm = instrument.fwhm(instrumentname, w_obs, as_bounded=True)[ + fwhm = instrument.fwhm(instrumentname, lam_obs, as_bounded=True)[ 0 ] / (1.0 + redshift) From 040798d6ed99a4e4db800887bbcc31933d2ab7cc Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 25 Jun 2024 10:48:41 -0400 Subject: [PATCH 20/27] Fix instrument and redshift being swapped in Features metadata --- pahfit/model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index 53c8518..e16ed8d 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -181,8 +181,8 @@ 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 From 032ed4b2e85db0e96934570135256d5bcdc56c10 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 25 Jun 2024 11:05:50 -0400 Subject: [PATCH 21/27] Fix style issue --- pahfit/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pahfit/model.py b/pahfit/model.py index e16ed8d..5258735 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -834,7 +834,7 @@ def _excluded_features(self, instrumentname, redshift, lam_obs=None): return is_outside & is_excludable - def _set_up_fitter( + def _set_up_fitter( self, instrumentname, redshift, lam=None, use_instrument_fwhm=True ): """Convert features table to Fitter instance, set self.fitter. From f53c8db69c5e0b7a15609cc2c24c53b5b1b54eaa Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 25 Jun 2024 13:57:12 -0400 Subject: [PATCH 22/27] Put fitter modules in pahfit.fitters. Rename apfitter to ap_fitter. --- docs/index.rst | 2 +- pahfit/feature_strengths.py | 2 +- pahfit/fitters/__init__.py | 0 pahfit/{component_models.py => fitters/ap_components.py} | 0 pahfit/{apfitter.py => fitters/ap_fitter.py} | 4 ++-- pahfit/{ => fitters}/fitter.py | 0 pahfit/model.py | 4 ++-- 7 files changed, 6 insertions(+), 6 deletions(-) create mode 100644 pahfit/fitters/__init__.py rename pahfit/{component_models.py => fitters/ap_components.py} (100%) rename pahfit/{apfitter.py => fitters/ap_fitter.py} (99%) rename pahfit/{ => fitters}/fitter.py (100%) diff --git a/docs/index.rst b/docs/index.rst index 905cf99..97952a0 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -94,4 +94,4 @@ Reference API Component models not provided by astropy.models. -.. automodapi:: pahfit.component_models +.. automodapi:: pahfit.fitters.ap_components diff --git a/pahfit/feature_strengths.py b/pahfit/feature_strengths.py index f1c6722..40eba46 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/fitters/__init__.py b/pahfit/fitters/__init__.py new file mode 100644 index 0000000..e69de29 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/apfitter.py b/pahfit/fitters/ap_fitter.py similarity index 99% rename from pahfit/apfitter.py rename to pahfit/fitters/ap_fitter.py index 9f11f33..92760bf 100644 --- a/pahfit/apfitter.py +++ b/pahfit/fitters/ap_fitter.py @@ -1,6 +1,6 @@ -from pahfit.fitter import Fitter +from pahfit.fitters.fitter import Fitter from pahfit.errors import PAHFITModelError -from pahfit.component_models import ( +from pahfit.fitters.ap_components import ( BlackBody1D, ModifiedBlackBody1D, S07_attenuation, diff --git a/pahfit/fitter.py b/pahfit/fitters/fitter.py similarity index 100% rename from pahfit/fitter.py rename to pahfit/fitters/fitter.py diff --git a/pahfit/model.py b/pahfit/model.py index 5258735..e6d4c4e 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -11,8 +11,8 @@ from pahfit.features import Features from pahfit import instrument from pahfit.errors import PAHFITModelError -from pahfit.component_models import BlackBody1D, S07_attenuation -from pahfit.apfitter import APFitter +from pahfit.fitters.ap_components import BlackBody1D, S07_attenuation +from pahfit.fitters.ap_fitter import APFitter class Model: From a990b47e95334822fedb1b76e8bba10eaf235fc0 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Wed, 3 Jul 2024 14:43:40 -0400 Subject: [PATCH 23/27] Fix astropy contribution guide link --- docs/index.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/index.rst b/docs/index.rst index 97952a0..7073faf 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -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 `_ From 2dc3c801f3e3957cb9342ba3dd98598ebc319806 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Wed, 3 Jul 2024 15:33:23 -0400 Subject: [PATCH 24/27] Replace links to PAHFIT classic --- CHANGES.rst | 2 +- README.rst | 2 +- docs/index.rst | 2 +- docs/version_description/version_description.rst | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 7d75f0a..c94f28e 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 dc92abe..1a02367 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 7073faf..e607236 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 `_ diff --git a/docs/version_description/version_description.rst b/docs/version_description/version_description.rst index 609592d..882f12a 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 ------------ From 5d22489f95e1346a07f9c833cb6c8757c53cc194 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Wed, 3 Jul 2024 16:28:26 -0400 Subject: [PATCH 25/27] Rephrase comment explaining redshift adjustment of fwhm --- pahfit/model.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index e6d4c4e..e2ef45e 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -903,20 +903,23 @@ def cleaned(features_tuple3): ) elif kind == "line": - if use_instrument_fwhm: - # one caveat here: redshift. Correct way to do it: - # 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) + # be careful with lines that have masked FWHM values here + fwhm_masked = row['fwhm'].ndim == 0 + if use_instrument_fwhm or fwhm_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) - fwhm = instrument.fwhm(instrumentname, lam_obs, as_bounded=True)[ - 0 - ] / (1.0 + redshift) + 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 From da9ed249803028088555750c532a26dd840a3ea2 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Mon, 8 Jul 2024 19:26:53 -0400 Subject: [PATCH 26/27] Replace more cases of "wav" by "lam", fix redshift bug in plot --- pahfit/model.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index e2ef45e..2293580 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -752,12 +752,12 @@ def tabulate( 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 + lam = np.arange(wmin, wmax, wfwhm / 2) * u.micron elif isinstance(wavelengths, Spectrum1D): - wav = wavelengths.spectral_axis.to(u.micron) + lam = wavelengths.spectral_axis.to(u.micron) / (1 + z) else: # any other iterable will be accepted and converted to array - wav = np.asarray(wavelengths) * u.micron + lam = np.asarray(wavelengths) * u.micron # apply feature mask, make sub model, and set up functional if feature_mask is not None: @@ -768,7 +768,7 @@ def tabulate( # if nothing is in range, return early with zeros if len(features_to_use) == 0: return Spectrum1D( - spectral_axis=wav, flux=np.zeros(wav.shape) * u.dimensionless_unscaled + spectral_axis=lam, flux=np.zeros(lam.shape) * u.dimensionless_unscaled ) alt_model = Model(features_to_use) @@ -783,12 +783,10 @@ def tabulate( fitter = alt_model.fitter except PAHFITModelError: return Spectrum1D( - spectral_axis=wav, flux=np.zeros(wav.shape) * u.dimensionless_unscaled + spectral_axis=lam, flux=np.zeros(lam.shape) * u.dimensionless_unscaled ) - # shift the "observed wavelength grid" to "physical wavelength grid - wav /= 1 + z - flux_values = fitter.evaluate(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) @@ -798,7 +796,7 @@ def tabulate( 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 _excluded_features(self, instrumentname, redshift, lam_obs=None): """Determine excluded features Based on instrument wavelength range. From 4e16008d1052541a78096ab809fe357485fea8ef Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 9 Jul 2024 16:11:11 -0400 Subject: [PATCH 27/27] Use "is np.ma.masked" when interpreting fwhm column --- pahfit/model.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index 2293580..83380b0 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -292,7 +292,7 @@ def amp_guess(row, fwhm): if calc_line_fwhm: for row_index in np.where(self.features["kind"] == "line")[0]: row = self.features[row_index] - if np.ma.is_masked(row["fwhm"]): + if row["fwhm"] is np.ma.masked: self.features[row_index]["fwhm"] = ( line_fwhm_guess(row), np.nan, @@ -902,8 +902,7 @@ def cleaned(features_tuple3): elif kind == "line": # be careful with lines that have masked FWHM values here - fwhm_masked = row['fwhm'].ndim == 0 - if use_instrument_fwhm or fwhm_masked: + 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 @@ -922,9 +921,14 @@ def cleaned(features_tuple3): # decide if scalar (fixed) or tuple (fitted fwhm # between upper and lower fwhm limits, happens in # case of overlapping instruments) - if np.ma.is_masked(fwhm): - fwhm = fwhm[0] + 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(