From 8bc62a068152e9f7b15cce8c98152dd315448404 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Thu, 18 Apr 2024 19:39:01 -0400 Subject: [PATCH 01/68] Use importlib.resources instead of pkg_resources (from setuptools) --- pahfit/features/features.py | 5 ++--- pahfit/helpers.py | 12 +++++------- pahfit/instrument.py | 5 ++--- 3 files changed, 9 insertions(+), 13 deletions(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index da5550e5..67302ef1 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -22,7 +22,7 @@ from astropy.table import vstack, Table, TableAttribute from astropy.io.misc.yaml import yaml import astropy.units as u -from pkg_resources import resource_filename +from importlib import resources from pahfit.errors import PAHFITFeatureError from pahfit.features.features_format import BoundedMaskedColumn, BoundedParTableFormatter @@ -158,8 +158,7 @@ def _read_scipack(cls, file): feat_tables = dict() if not os.path.isfile(file): - pack_path = resource_filename("pahfit", "packs/science") - file = os.path.join(pack_path, file) + file = resources.files("pahfit") / "packs/science" / file try: with open(file) as fd: scipack = yaml.load(fd, Loader=UniqueKeyLoader) diff --git a/pahfit/helpers.py b/pahfit/helpers.py index cc42836e..65c46065 100644 --- a/pahfit/helpers.py +++ b/pahfit/helpers.py @@ -1,5 +1,5 @@ import os -import pkg_resources +from importlib import resources from specutils import Spectrum1D @@ -32,14 +32,13 @@ def find_packfile(packfile): if os.path.isfile(packfile): packfile_found = packfile else: - pack_path = pkg_resources.resource_filename("pahfit", "packs/science/") - test_packfile = "{}/{}".format(pack_path, packfile) + test_packfile = resources.files("pahfit") / "packs/science" / packfile if os.path.isfile(test_packfile): packfile_found = test_packfile else: raise ValueError("Input packfile {} not found".format(packfile)) - return packfile_found + return str(packfile_found) def read_spectrum(specfile, format=None): @@ -63,8 +62,7 @@ def read_spectrum(specfile, format=None): """ # resolve filename if not os.path.isfile(specfile): - pack_path = pkg_resources.resource_filename("pahfit", "data/") - test_specfile = "{}/{}".format(pack_path, specfile) + test_specfile = resources.files("pahfit") / "data" / specfile if os.path.isfile(test_specfile): specfile = test_specfile else: @@ -75,7 +73,7 @@ def read_spectrum(specfile, format=None): tformat = None # process user-specified or filename extension based format if format is None: - suffix = specfile.split(".")[-1].lower() + suffix = specfile.name.split(".")[-1].lower() if suffix == "ecsv": tformat = "ECSV" elif suffix == "ipac": diff --git a/pahfit/instrument.py b/pahfit/instrument.py index 3828d9c7..c55bdb8a 100644 --- a/pahfit/instrument.py +++ b/pahfit/instrument.py @@ -11,11 +11,10 @@ import os from pathlib import Path -import glob import numpy as np from numpy.polynomial import Polynomial from astropy.io.misc import yaml -from pkg_resources import resource_filename +from importlib import resources from pahfit.errors import PAHFITPackError, PAHFITWarning from warnings import warn @@ -25,7 +24,7 @@ def read_instrument_packs(): """Read all instrument packs into the 'packs' variable.""" global packs - for pack in glob.glob(resource_filename("pahfit", "packs/instrument/*.yaml")): + for pack in (resources.files("pahfit") / "packs/instrument").glob("*.yaml"): try: with open(pack) as fd: p = yaml.load(fd) From d53422fecbccfed622f9d3059c5628d96914b0ae Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Wed, 8 May 2024 21:40:49 -0400 Subject: [PATCH 02/68] Simplify feature bounds using np.nan instead of masking Partially masked arrays are not well supported, and required elaborate workaround for pretty printing. Here we use a structured array type, and represent "fixed" bounds with np.nan instead of masking. I.e. non is interpreted as the same as negative or positive infinity, unless both bound are nan, which indicates a fixed parameter. --- pahfit/features/features.py | 13 +++++---- pahfit/features/features_format.py | 47 +++++++++--------------------- pahfit/features/util.py | 13 ++++----- 3 files changed, 27 insertions(+), 46 deletions(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index da5550e5..1d6d5019 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -82,7 +82,7 @@ def value_bounds(val, bounds): if val is None: val = np.ma.masked if not bounds: - return (val,) + 2 * (np.ma.masked,) # Fixed + return (val,) + 2 * (np.nan,) # Fixed ret = [val] for i, b in enumerate(bounds): if isinstance(b, str): @@ -123,7 +123,8 @@ 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', 'geometry', 'model')) # String attributes (no bounds) - + _bounds_dtype = np.dtype([("val", "f4"), ("min", "f4"), ("max", "f4")]) + @classmethod def read(cls, file, *args, **kwargs): """Read a table from file. @@ -308,7 +309,7 @@ def _construct_table(cls, inp: dict): tables = [] for (kind, features) in inp.items(): kind_params = cls._kind_params[kind] # All params for this kind - rows = [] + rows, dtypes = [], [] for (name, params) in features.items(): for missing in kind_params - params.keys(): if missing in cls._no_bounds: @@ -316,8 +317,8 @@ def _construct_table(cls, inp: dict): else: params[missing] = value_bounds(0.0, bounds=(0.0, None)) rows.append(dict(name=name, **params)) - table_columns = rows[0].keys() - t = cls(rows, names=table_columns) + dtypes.append(None if name in cls._no_bounds else cls._bounds_dtypes) + t = cls(rows, names=rows[0].keys(), dtype=dtypes) for p in cls._kind_params[kind]: if p not in cls._no_bounds: t[p].info.format = "0.4g" # Nice format (customized by Formatter) @@ -352,7 +353,7 @@ def mask_feature(self, name, mask_value=True): pass else: # mask only the value, not the bounds - row[col_name].mask[0] = mask_value + row[col_name].mask['val'] = mask_value def unmask_feature(self, name): """Remove the mask for all parameters of a feature.""" diff --git a/pahfit/features/features_format.py b/pahfit/features/features_format.py index 78369727..9526b6ed 100644 --- a/pahfit/features/features_format.py +++ b/pahfit/features/features_format.py @@ -1,53 +1,34 @@ -import numpy.ma as ma -from astropy.table import MaskedColumn +import numpy as np from astropy.table.pprint import TableFormatter - # * Special table formatting for bounded (val, min, max) values def fmt_func(fmt): - def _fmt(v): - if ma.is_masked(v[0]): - return " " - if ma.is_masked(v[1]): - return f"{v[0]:{fmt}} (Fixed)" - return f"{v[0]:{fmt}} ({v[1]:{fmt}}, {v[2]:{fmt}})" - + def _fmt(x): + ret = f"{x['val']:{fmt}}" + if np.isnan(x['min']) and np.isnan(x['max']): + return ret + " (fixed)" + else: + mn = ("-∞" if np.isnan(x['min']) or x['min'] == -np.inf + else f"{x['min']:{fmt}}") + mx = ("∞" if np.isnan(x['max']) or x['max'] == np.inf + else f"{x['max']:{fmt}}") + return f"{ret} ({mn}, {mx})" return _fmt -class BoundedMaskedColumn(MaskedColumn): - """Masked column which can be toggled to group rows into one item - for formatting. To be set as Table's `MaskedColumn'. - """ - - _omit_shape = False - - @property - def shape(self): - sh = super().shape - return sh[0:-1] if self._omit_shape and len(sh) > 1 else sh - - def is_fixed(self): - return ma.getmask(self)[:, 1:].all(1) - - class BoundedParTableFormatter(TableFormatter): """Format bounded parameters. Bounded parameters are 3-field structured arrays, with fields - 'var', 'min', and 'max'. To be set as Table's `TableFormatter'. + 'val', 'min', and 'max'. To be set as Table's `TableFormatter'. """ - def _pformat_table(self, table, *args, **kwargs): bpcols = [] try: - colsh = [(col, col.shape) for col in table.columns.values()] - BoundedMaskedColumn._omit_shape = True - for col, sh in colsh: - if len(sh) == 2 and sh[1] == 3: + for col in table.columns.values(): + if len(col.dtype) == 3 # bounded! bpcols.append((col, col.info.format)) col.info.format = fmt_func(col.info.format or "g") return super()._pformat_table(table, *args, **kwargs) finally: - BoundedMaskedColumn._omit_shape = False for col, fmt in bpcols: col.info.format = fmt diff --git a/pahfit/features/util.py b/pahfit/features/util.py index 68513aa0..92c18020 100644 --- a/pahfit/features/util.py +++ b/pahfit/features/util.py @@ -1,29 +1,28 @@ """pahfit.util General pahfit.features utility functions.""" import numpy as np -import numpy.ma as ma def bounded_is_missing(val): """Return a mask array indicating which of the bounded values are missing. A missing bounded value has a masked value.""" - return ma.getmask(val)[..., 0] + return getattr(a['val'], 'mask', None) or np.zeros_like(a['val'], dtype=bool) def bounded_is_fixed(val): """Return a mask array indicating which of the bounded values are fixed. A fixed bounded value has masked bounds.""" - return ma.getmask(val)[..., -2:].all(-1) + return np.isnan(val['min']) & np.isnan(val['max']) def bounded_min(val): """Return the minimum of each bounded value passed. Either the lower bound, or, if no such bound is set, the value itself.""" - lower = val[..., 1] - return np.where(lower, lower, val[..., 0]) + lower = val['min'] + return np.where(lower, lower, val['val']) def bounded_max(val): """Return the maximum of each bounded value passed. Either the upper bound, or, if no such bound is set, the value itself.""" - upper = val[..., 2] - return np.where(upper, upper, val[..., 0]) + upper = val['max'] + return np.where(upper, upper, val['val']) From 3017dda4828604815029ade8a6d2b37413db5ff5 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 9 May 2024 09:14:45 -0400 Subject: [PATCH 03/68] Add missing colon --- pahfit/features/features_format.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pahfit/features/features_format.py b/pahfit/features/features_format.py index 9526b6ed..6bb8eec1 100644 --- a/pahfit/features/features_format.py +++ b/pahfit/features/features_format.py @@ -25,7 +25,7 @@ def _pformat_table(self, table, *args, **kwargs): bpcols = [] try: for col in table.columns.values(): - if len(col.dtype) == 3 # bounded! + if len(col.dtype) == 3: # bounded! bpcols.append((col, col.info.format)) col.info.format = fmt_func(col.info.format or "g") return super()._pformat_table(table, *args, **kwargs) From ec979326e3fff065e6745bd71e87bf98a78b7848 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 9 May 2024 10:20:15 -0400 Subject: [PATCH 04/68] Remove BoundedMaskedColumn A normal masked column using structured numpy arrays and nan for both bounds when fixed now suffices. --- pahfit/features/features.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 1d6d5019..6b11a6a7 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -24,7 +24,7 @@ import astropy.units as u from pkg_resources import resource_filename from pahfit.errors import PAHFITFeatureError -from pahfit.features.features_format import BoundedMaskedColumn, BoundedParTableFormatter +from pahfit.features.features_format import BoundedParTableFormatter class UniqueKeyLoader(yaml.SafeLoader): @@ -109,7 +109,6 @@ class Features(Table): """ TableFormatter = BoundedParTableFormatter - MaskedColumn = BoundedMaskedColumn param_covar = TableAttribute(default=[]) _kind_params = {'starlight': {'temperature', 'tau'}, From afc0bc9e66b3ba7195fe0a0210dadb4707df0c60 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 9 May 2024 10:59:26 -0400 Subject: [PATCH 05/68] Fix bounded_is_missing arg usage --- pahfit/features/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pahfit/features/util.py b/pahfit/features/util.py index 92c18020..d9493f94 100644 --- a/pahfit/features/util.py +++ b/pahfit/features/util.py @@ -5,7 +5,7 @@ def bounded_is_missing(val): """Return a mask array indicating which of the bounded values are missing. A missing bounded value has a masked value.""" - return getattr(a['val'], 'mask', None) or np.zeros_like(a['val'], dtype=bool) + return getattr(val['val'], 'mask', None) or np.zeros_like(val['val'], dtype=bool) def bounded_is_fixed(val): From ef094afbde5619f6435ddbbb4cbe86b22db9e1b7 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 9 May 2024 11:08:43 -0400 Subject: [PATCH 06/68] features: doc improvements --- pahfit/features/features.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 6b11a6a7..c9e71b28 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -69,8 +69,8 @@ def value_bounds(val, bounds): Returns: ------- - The value, if unbounded, or a 3 element tuple (value, min, max). - Any missing bound is replaced with the numpy `masked' value. + A 3 element tuple (value, min, max). + Any missing bound is replaced with the numpy.nan value. Raises: ------- @@ -82,7 +82,7 @@ def value_bounds(val, bounds): if val is None: val = np.ma.masked if not bounds: - return (val,) + 2 * (np.nan,) # Fixed + return (val,) + 2 * (np.nan,) # (val,nan,nan) indicates fixed ret = [val] for i, b in enumerate(bounds): if isinstance(b, str): From 4fefb0adec781ab75c076886bbd443534dd49b06 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 9 May 2024 12:56:39 -0400 Subject: [PATCH 07/68] Fix dtype var mention --- pahfit/features/features.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index c9e71b28..84067725 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -123,7 +123,7 @@ class Features(Table): _param_attrs = set(('value', 'bounds')) # Each parameter can have these attributes _no_bounds = set(('name', 'group', 'geometry', 'model')) # String attributes (no bounds) _bounds_dtype = np.dtype([("val", "f4"), ("min", "f4"), ("max", "f4")]) - + @classmethod def read(cls, file, *args, **kwargs): """Read a table from file. @@ -316,7 +316,7 @@ def _construct_table(cls, inp: dict): else: params[missing] = value_bounds(0.0, bounds=(0.0, None)) rows.append(dict(name=name, **params)) - dtypes.append(None if name in cls._no_bounds else cls._bounds_dtypes) + dtypes.append(None if name in cls._no_bounds else cls._bounds_dtype) t = cls(rows, names=rows[0].keys(), dtype=dtypes) for p in cls._kind_params[kind]: if p not in cls._no_bounds: From 3c9fb7a1ae32be563383a0cebc5ecf91f22adb64 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 9 May 2024 13:29:39 -0400 Subject: [PATCH 08/68] _construct_table: fix dtypes formulation --- pahfit/features/features.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 84067725..65b06d76 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -316,7 +316,8 @@ def _construct_table(cls, inp: dict): else: params[missing] = value_bounds(0.0, bounds=(0.0, None)) rows.append(dict(name=name, **params)) - dtypes.append(None if name in cls._no_bounds else cls._bounds_dtype) + dtypes = [None if x in cls._no_bounds else cls._bounds_dtype for x in kind_params] + dtypes.insert(0, None) # For the name t = cls(rows, names=rows[0].keys(), dtype=dtypes) for p in cls._kind_params[kind]: if p not in cls._no_bounds: From 2729f89f1e2141f0d5d8314930824ce48cb8fa6c Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 9 May 2024 13:30:01 -0400 Subject: [PATCH 09/68] Make linter happier --- pahfit/features/features_format.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pahfit/features/features_format.py b/pahfit/features/features_format.py index 6bb8eec1..084bf199 100644 --- a/pahfit/features/features_format.py +++ b/pahfit/features/features_format.py @@ -1,6 +1,7 @@ import numpy as np from astropy.table.pprint import TableFormatter + # * Special table formatting for bounded (val, min, max) values def fmt_func(fmt): def _fmt(x): @@ -25,7 +26,7 @@ def _pformat_table(self, table, *args, **kwargs): bpcols = [] try: for col in table.columns.values(): - if len(col.dtype) == 3: # bounded! + if len(col.dtype) == 3: # bounded! bpcols.append((col, col.info.format)) col.info.format = fmt_func(col.info.format or "g") return super()._pformat_table(table, *args, **kwargs) From ce49a70026c7a013c71ae18ca4840c21476c7aeb Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 9 May 2024 13:32:18 -0400 Subject: [PATCH 10/68] _construct_table: hard-code str for column types for no-bounds cols --- pahfit/features/features.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 65b06d76..c453c679 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -316,8 +316,8 @@ def _construct_table(cls, inp: dict): else: params[missing] = value_bounds(0.0, bounds=(0.0, None)) rows.append(dict(name=name, **params)) - dtypes = [None if x in cls._no_bounds else cls._bounds_dtype for x in kind_params] - dtypes.insert(0, None) # For the name + dtypes = [str if x in cls._no_bounds else cls._bounds_dtype for x in kind_params] + dtypes.insert(0, str) # For the name t = cls(rows, names=rows[0].keys(), dtype=dtypes) for p in cls._kind_params[kind]: if p not in cls._no_bounds: From 1cd8404cd7ffaba85a40d4525625cd48286cc492 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 9 May 2024 13:55:04 -0400 Subject: [PATCH 11/68] Fix dtype and format for structured array Features cells --- pahfit/features/features.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index c453c679..021fcfe3 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -121,7 +121,7 @@ class Features(Table): _units = {'temperature': u.K, 'wavelength': u.um, 'fwhm': u.um} _group_attrs = set(('bounds', 'features', 'kind')) # group-level attributes _param_attrs = set(('value', 'bounds')) # Each parameter can have these attributes - _no_bounds = set(('name', 'group', 'geometry', 'model')) # String attributes (no bounds) + _no_bounds = set(('name', 'group', 'kind', 'geometry', 'model')) # str attributes (no bounds) _bounds_dtype = np.dtype([("val", "f4"), ("min", "f4"), ("max", "f4")]) @classmethod @@ -316,12 +316,9 @@ def _construct_table(cls, inp: dict): else: params[missing] = value_bounds(0.0, bounds=(0.0, None)) rows.append(dict(name=name, **params)) - dtypes = [str if x in cls._no_bounds else cls._bounds_dtype for x in kind_params] - dtypes.insert(0, str) # For the name - t = cls(rows, names=rows[0].keys(), dtype=dtypes) - for p in cls._kind_params[kind]: - if p not in cls._no_bounds: - t[p].info.format = "0.4g" # Nice format (customized by Formatter) + param_names = rows[0].keys() + dtypes = [str if x in cls._no_bounds else cls._bounds_dtype for x in param_names] + t = cls(rows, names=param_names, dtype=dtypes) tables.append(t) tables = vstack(tables) for cn, col in tables.columns.items(): From cf180eef62210fa58edf5a719bba2f7880ef6091 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 9 May 2024 14:07:18 -0400 Subject: [PATCH 12/68] Suppress mention of [val, min, max] in column headers --- pahfit/features/features_format.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pahfit/features/features_format.py b/pahfit/features/features_format.py index 084bf199..44d3cc1d 100644 --- a/pahfit/features/features_format.py +++ b/pahfit/features/features_format.py @@ -33,3 +33,7 @@ def _pformat_table(self, table, *args, **kwargs): finally: for col, fmt in bpcols: col.info.format = fmt + + def _name_and_structure(self, name, *args): + "Simplified column name: no val, min, max needed." + return name From be7cca1d212b5cd29b7a0d8120123cf9df50f540 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 9 May 2024 14:54:38 -0400 Subject: [PATCH 13/68] Features: Strip dtype from default __repr__ --- pahfit/features/features.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 021fcfe3..eca9732e 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -355,3 +355,7 @@ def mask_feature(self, name, mask_value=True): def unmask_feature(self, name): """Remove the mask for all parameters of a feature.""" self.mask_feature(name, mask_value=False) + + def _base_repr_(self, *args, **kwargs): + """Omit dtype on self-print.""" + return super()._base_repr_(*args, ** kwargs | dict(show_dtype=False)) From 19e9c25b02fab79bd0a1418b8ac99dc0cc47dfa7 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 9 May 2024 14:55:04 -0400 Subject: [PATCH 14/68] Allow meta['pahfit_format'] on table and columns --- pahfit/features/features_format.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/pahfit/features/features_format.py b/pahfit/features/features_format.py index 44d3cc1d..2871e936 100644 --- a/pahfit/features/features_format.py +++ b/pahfit/features/features_format.py @@ -3,7 +3,11 @@ # * Special table formatting for bounded (val, min, max) values -def fmt_func(fmt): +def fmt_func(fmt: str): + """Format bounded variables specially.""" + if fmt.startswith('%'): + fmt = fmt[1:] + def _fmt(x): ret = f"{x['val']:{fmt}}" if np.isnan(x['min']) and np.isnan(x['max']): @@ -24,11 +28,13 @@ class BoundedParTableFormatter(TableFormatter): """ def _pformat_table(self, table, *args, **kwargs): bpcols = [] + tlfmt = table.meta.get('pahfit_format') try: for col in table.columns.values(): if len(col.dtype) == 3: # bounded! bpcols.append((col, col.info.format)) - col.info.format = fmt_func(col.info.format or "g") + fmt = col.meta.get('pahfit_format') or tlfmt or "g" + col.info.format = fmt_func(fmt) return super()._pformat_table(table, *args, **kwargs) finally: for col, fmt in bpcols: From 646ffd1d2a1de7994db9404f8e8c26915cef5025 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 27 Oct 2022 10:07:08 -0400 Subject: [PATCH 15/68] Mask unused value for linter --- 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 67302ef1..3ad97576 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -30,7 +30,7 @@ class UniqueKeyLoader(yaml.SafeLoader): def construct_mapping(self, node, deep=False): mapping = set() - for key_node, value_node in node.value: + for key_node, _ in node.value: key = self.construct_object(key_node, deep=deep) if key in mapping: raise PAHFITFeatureError(f"Duplicate {key!r} key found in YAML.") From 6bd5e1db07f2093ef82af101f92bda605de9811c Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Fri, 17 Mar 2023 15:49:18 -0400 Subject: [PATCH 16/68] Expose KIND_PARAMS and PARAM_UNITS as module level constants --- pahfit/features/features.py | 63 +++++++++++++++++++++++-------------- pahfit/units.py | 19 +++++++++++ 2 files changed, 59 insertions(+), 23 deletions(-) create mode 100644 pahfit/units.py diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 3ad97576..48629133 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -25,6 +25,20 @@ from importlib import resources from pahfit.errors import PAHFITFeatureError from pahfit.features.features_format import BoundedMaskedColumn, BoundedParTableFormatter +from pahfit.units import UNITS + +# Feature kinds and associated parameters +KIND_PARAMS = {'starlight': {'temperature', 'tau'}, + 'dust_continuum': {'temperature', 'tau'}, + 'line': {'wavelength', 'power'}, # 'fwhm', Instrument Pack detail! + 'dust_feature': {'wavelength', 'fwhm', 'power'}, + 'attenuation': {'model', 'tau', 'geometry'}, + 'absorption': {'wavelength', 'fwhm', 'tau', 'geometry'}} + +# Parameter default units: flux density/intensity/power (other units determined on fit) +PARAM_UNITS = {'temperature': UNITS.temperature, + 'wavelength': UNITS.wavelength, + 'fwhm': UNITS.wavelength} class UniqueKeyLoader(yaml.SafeLoader): @@ -102,26 +116,24 @@ def value_bounds(val, bounds): class Features(Table): - """A class for holding PAHFIT features and their associated - parameter information. Note that each parameter has an associated - `kind', and that each kind has an associated set of allowable - parameters (see _kind_params, below). + """A class for holding a table of PAHFIT features and associated + parameter information. + + Note that each parameter has an associated `kind', and that each + kind has an associated set of allowable parameters (see + `KIND_PARAMS`). + + See Also + -------- + `~astropy.table.Table`: The parent table class. """ TableFormatter = BoundedParTableFormatter MaskedColumn = BoundedMaskedColumn param_covar = TableAttribute(default=[]) - _kind_params = {'starlight': {'temperature', 'tau'}, - 'dust_continuum': {'temperature', 'tau'}, - 'line': {'wavelength', 'power'}, # 'fwhm', Instrument Pack detail! - 'dust_feature': {'wavelength', 'fwhm', 'power'}, - 'attenuation': {'model', 'tau', 'geometry'}, - 'absorption': {'wavelength', 'fwhm', 'tau', 'geometry'}} - - _units = {'temperature': u.K, 'wavelength': u.um, 'fwhm': u.um} - _group_attrs = set(('bounds', 'features', 'kind')) # group-level attributes - _param_attrs = set(('value', 'bounds')) # Each parameter can have these attributes + _param_attrs = set(('value', 'bounds', 'tied')) # params can have these attributes + _group_attrs = set(('bounds', 'features', 'kind', 'tied')) # group-level attributes _no_bounds = set(('name', 'group', 'geometry', 'model')) # String attributes (no bounds) @classmethod @@ -178,7 +190,7 @@ def _read_scipack(cls, file): raise PAHFITFeatureError(f"No kind found for {name}\n\t{file}") try: - valid_params = cls._kind_params[kind] + valid_params = KIND_PARAMS[kind] except KeyError: raise PAHFITFeatureError(f"Unknown kind {kind} for {name}\n\t{file}") unknown_params = [x for x in keys @@ -263,7 +275,7 @@ def _add_feature(cls, kind: str, t: dict, name: str, *, t[kind][name]['group'] = group t[kind][name]['kind'] = kind for (param, val) in pars.items(): - if param not in cls._kind_params[kind]: + if param not in KIND_PARAMS[kind]: continue if isinstance(val, dict): # A param attribute dictionary unknown_attrs = [x for x in val.keys() if x not in cls._param_attrs] @@ -306,10 +318,12 @@ def _construct_table(cls, inp: dict): """ tables = [] for (kind, features) in inp.items(): - kind_params = cls._kind_params[kind] # All params for this kind + if kind == "_ratios": + continue + kp = KIND_PARAMS[kind] # All params for this kind rows = [] for (name, params) in features.items(): - for missing in kind_params - params.keys(): + for missing in kp - params.keys(): if missing in cls._no_bounds: params[missing] = 0.0 else: @@ -317,15 +331,18 @@ def _construct_table(cls, inp: dict): rows.append(dict(name=name, **params)) table_columns = rows[0].keys() t = cls(rows, names=table_columns) - for p in cls._kind_params[kind]: + for p in KIND_PARAMS[kind]: if p not in cls._no_bounds: t[p].info.format = "0.4g" # Nice format (customized by Formatter) tables.append(t) tables = vstack(tables) for cn, col in tables.columns.items(): - if cn in cls._units: - col.unit = cls._units[cn] - tables.add_index('name') + if cn in PARAM_UNITS: + col.unit = PARAM_UNITS[cn] + cls._index_table(tables) + + if '_ratios' in inp: + tables.meta['_ratios'] = inp['_ratios'] return tables def mask_feature(self, name, mask_value=True): @@ -344,7 +361,7 @@ def mask_feature(self, name, mask_value=True): """ row = self.loc[name] - relevant_params = self._kind_params[row['kind']] + relevant_params = KIND_PARAMS[row['kind']] for col_name in relevant_params: if col_name in self._no_bounds: # these are all strings, so can't mask diff --git a/pahfit/units.py b/pahfit/units.py new file mode 100644 index 00000000..c95549b7 --- /dev/null +++ b/pahfit/units.py @@ -0,0 +1,19 @@ +import astropy.units as u +from astropy.units import CompositeUnit + +from enum import Enum + +# Working/default unitParameter default units: flux density/intensity/power +# These are PAHFITs default science packs parameter and output units +class UNITS(Enum): + temperature = u.K + wavelength = u.um + fwhm = u.um + flux_density = u.mJy + flux_power = CompositeUnit(1e-22, (u.W, u.m), (1, -2)) + intensity = u.MJy / u.sr + intensity_power = CompositeUnit(1e-10, (u.W, u.m, u.sr), (1, -2, -1)) + +# integrated power units of 1e-22 W/m^2 (from flux) corresponds to the +# unit 1e-10 W/m^2/sr (from intensity) if it occurs uniformly over a +# solid angle 0.21" on a side (about a small JWST IFU pixel) From 3d0aae5e34441100bc0ed746274dfaea5cd97d04 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 15 Dec 2022 12:50:15 -0500 Subject: [PATCH 17/68] UNITS: remove fwhm, add solid_angle --- pahfit/units.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pahfit/units.py b/pahfit/units.py index c95549b7..8893f244 100644 --- a/pahfit/units.py +++ b/pahfit/units.py @@ -1,19 +1,19 @@ +from enum import Enum import astropy.units as u from astropy.units import CompositeUnit -from enum import Enum # Working/default unitParameter default units: flux density/intensity/power # These are PAHFITs default science packs parameter and output units class UNITS(Enum): temperature = u.K wavelength = u.um - fwhm = u.um flux_density = u.mJy flux_power = CompositeUnit(1e-22, (u.W, u.m), (1, -2)) + solid_angle = u.sr intensity = u.MJy / u.sr intensity_power = CompositeUnit(1e-10, (u.W, u.m, u.sr), (1, -2, -1)) -# integrated power units of 1e-22 W/m^2 (from flux) corresponds to the -# unit 1e-10 W/m^2/sr (from intensity) if it occurs uniformly over a -# solid angle 0.21" on a side (about a small JWST IFU pixel) +# Note: integrated power units of 1e-22 W/m^2 (from flux) corresponds +# to the unit 1e-10 W/m^2/sr (from intensity) if it occurs uniformly +# over a solid angle 0.21" on a side (about a small JWST IFU pixel) From a2740971834489682e132e9940beddc12c00adb3 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Thu, 15 Dec 2022 13:04:03 -0500 Subject: [PATCH 18/68] Bug: Enum members need .value to extract the unit value --- pahfit/features/features.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 48629133..68560ebf 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -36,9 +36,9 @@ 'absorption': {'wavelength', 'fwhm', 'tau', 'geometry'}} # Parameter default units: flux density/intensity/power (other units determined on fit) -PARAM_UNITS = {'temperature': UNITS.temperature, - 'wavelength': UNITS.wavelength, - 'fwhm': UNITS.wavelength} +PARAM_UNITS = {'temperature': UNITS.temperature.value, + 'wavelength': UNITS.wavelength.value, + 'fwhm': UNITS.wavelength.value} class UniqueKeyLoader(yaml.SafeLoader): From 548ded1ed3d6e23f83457f5442b5515959322521 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Fri, 16 Dec 2022 08:11:33 -0500 Subject: [PATCH 19/68] units: switch to simple module constants --- pahfit/features/features.py | 8 ++++---- pahfit/units.py | 16 +++++++--------- 2 files changed, 11 insertions(+), 13 deletions(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 68560ebf..7af056fa 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -25,7 +25,7 @@ from importlib import resources from pahfit.errors import PAHFITFeatureError from pahfit.features.features_format import BoundedMaskedColumn, BoundedParTableFormatter -from pahfit.units import UNITS +import pahfit.units # Feature kinds and associated parameters KIND_PARAMS = {'starlight': {'temperature', 'tau'}, @@ -36,9 +36,9 @@ 'absorption': {'wavelength', 'fwhm', 'tau', 'geometry'}} # Parameter default units: flux density/intensity/power (other units determined on fit) -PARAM_UNITS = {'temperature': UNITS.temperature.value, - 'wavelength': UNITS.wavelength.value, - 'fwhm': UNITS.wavelength.value} +PARAM_UNITS = {'temperature': pahfit.units.temperature, + 'wavelength': pahfit.units.wavelength, + 'fwhm': pahfit.units.wavelength} class UniqueKeyLoader(yaml.SafeLoader): diff --git a/pahfit/units.py b/pahfit/units.py index 8893f244..103e42f3 100644 --- a/pahfit/units.py +++ b/pahfit/units.py @@ -1,18 +1,16 @@ -from enum import Enum import astropy.units as u from astropy.units import CompositeUnit # Working/default unitParameter default units: flux density/intensity/power # These are PAHFITs default science packs parameter and output units -class UNITS(Enum): - temperature = u.K - wavelength = u.um - flux_density = u.mJy - flux_power = CompositeUnit(1e-22, (u.W, u.m), (1, -2)) - solid_angle = u.sr - intensity = u.MJy / u.sr - intensity_power = CompositeUnit(1e-10, (u.W, u.m, u.sr), (1, -2, -1)) +temperature = u.K +wavelength = u.um +flux_density = u.mJy +flux_power = CompositeUnit(1e-22, (u.W, u.m), (1, -2)) +solid_angle = u.sr +intensity = u.MJy / u.sr +intensity_power = CompositeUnit(1e-10, (u.W, u.m, u.sr), (1, -2, -1)) # Note: integrated power units of 1e-22 W/m^2 (from flux) corresponds # to the unit 1e-10 W/m^2/sr (from intensity) if it occurs uniformly From 8d637c434717131f666d4bfa335eb06093b94e63 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Fri, 17 Mar 2023 16:11:27 -0400 Subject: [PATCH 20/68] Remove 'tied' param/group type (for another PR) --- pahfit/features/features.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 7af056fa..2a8a34cb 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -132,8 +132,8 @@ class Features(Table): MaskedColumn = BoundedMaskedColumn param_covar = TableAttribute(default=[]) - _param_attrs = set(('value', 'bounds', 'tied')) # params can have these attributes - _group_attrs = set(('bounds', 'features', 'kind', 'tied')) # group-level attributes + _param_attrs = set(('value', 'bounds')) # params can have these attributes + _group_attrs = set(('bounds', 'features', 'kind')) # group-level attributes _no_bounds = set(('name', 'group', 'geometry', 'model')) # String attributes (no bounds) @classmethod From 1157f89f389bc3ae5cb7df131c524ee3c93a5f6f Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Fri, 17 Mar 2023 16:18:47 -0400 Subject: [PATCH 21/68] _index_table: re-add --- pahfit/features/features.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 2a8a34cb..7704bfd8 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -345,6 +345,11 @@ def _construct_table(cls, inp: dict): tables.meta['_ratios'] = inp['_ratios'] return tables + @staticmethod + def _index_table(tbl): + for indx in ('name', 'group'): + tbl.add_index(indx) + def mask_feature(self, name, mask_value=True): """Mask all the parameters of a feature. From b4a5ba40503b2c279a15d790a4185bbe23e6cec0 Mon Sep 17 00:00:00 2001 From: JD Smith <93749+jdtsmith@users.noreply.github.com> Date: Fri, 17 Mar 2023 16:26:46 -0400 Subject: [PATCH 22/68] linter happiness --- pahfit/features/features.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 7704bfd8..85adaba0 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -21,7 +21,6 @@ import numpy as np from astropy.table import vstack, Table, TableAttribute from astropy.io.misc.yaml import yaml -import astropy.units as u from importlib import resources from pahfit.errors import PAHFITFeatureError from pahfit.features.features_format import BoundedMaskedColumn, BoundedParTableFormatter From 21c61f064d05990453a9592edd2320269e7627f0 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Wed, 8 May 2024 14:18:50 -0400 Subject: [PATCH 23/68] Set power unit to intensity_power amplitude for now --- pahfit/features/features.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pahfit/features/features.py b/pahfit/features/features.py index 85adaba0..db7d3007 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -35,9 +35,13 @@ 'absorption': {'wavelength', 'fwhm', 'tau', 'geometry'}} # Parameter default units: flux density/intensity/power (other units determined on fit) +# Note: power is actually amplitude for now. Change this to +# intensity_power when the power features are implemented. + PARAM_UNITS = {'temperature': pahfit.units.temperature, 'wavelength': pahfit.units.wavelength, - 'fwhm': pahfit.units.wavelength} + 'fwhm': pahfit.units.wavelength, + 'power': pahfit.units.intensity} class UniqueKeyLoader(yaml.SafeLoader): From e7e5c70a0a3f8fc6043ee6d86f7813623dcbbf99 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Thu, 9 May 2024 11:27:39 -0400 Subject: [PATCH 24/68] Convert input spectrum to internal units, error if not ~ MJy sr-1 --- pahfit/model.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index af4382fb..df0a6c51 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -12,6 +12,7 @@ from pahfit import instrument from pahfit.errors import PAHFITModelError from pahfit.component_models import BlackBody1D +from pahfit import units class Model: @@ -288,9 +289,12 @@ def amp_guess(row, fwhm): @staticmethod def _convert_spec_data(spec, z): - """ - Turn astropy quantities stored in Spectrum1D into fittable - numbers. + """Convert Spectrum1D Quantities to fittable numbers. + + The unit of the input spectrum has to be a multiple of MJy / sr, + the internal intensity unit. The output of this function + consists of simple unitless arrays (the numbers in these arrays + are assumed to be consistent with the internal units). Also corrects for redshift. @@ -300,10 +304,13 @@ def _convert_spec_data(spec, z): xz, yz, uncz: wavelength in micron, flux, uncertainty corrected for redshift + """ + if not spec.flux.unit.is_equivalent(units.intensity): + raise PAHFITModelError("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 - y = spec.flux.value - unc = spec.uncertainty.array + unc = (spec.uncertainty.array * spec.flux.unit).to(units.intensity).value # transform observed wavelength to "physical" wavelength xz = x / (1 + z) # wavelength shorter From 1a3360bd8bedf2487f6909e5c9b04bf0c2ac90ca Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Thu, 9 May 2024 17:19:29 -0400 Subject: [PATCH 25/68] State that MJy / sr unit restriction is temporary --- pahfit/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pahfit/model.py b/pahfit/model.py index df0a6c51..69d22f3b 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -307,7 +307,7 @@ def _convert_spec_data(spec, z): """ if not spec.flux.unit.is_equivalent(units.intensity): - raise PAHFITModelError("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 From 4573a35132a1bbdee4d339cc725c735bc23f57d4 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Thu, 9 May 2024 17:59:09 -0400 Subject: [PATCH 26/68] Convert example spectra to MJy / sr (temporary fix) --- pahfit/helpers.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/pahfit/helpers.py b/pahfit/helpers.py index 65c46065..fc84ee51 100644 --- a/pahfit/helpers.py +++ b/pahfit/helpers.py @@ -1,11 +1,14 @@ import os from importlib import resources -from specutils import Spectrum1D - 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"] @@ -81,8 +84,18 @@ def read_spectrum(specfile, format=None): else: tformat = format - return Spectrum1D.read(specfile, format=tformat) + s = Spectrum1D.read(specfile, format=tformat) + + # Convert to intensity units by assuming an arbitrary solid angle + # for now. To be removed when dual unit support (intensity and flux + # density) is supported. + if s.flux.unit.is_equivalent(units.flux_density): + solid_angle = (3 * u.arcsec)**2 + alt_flux = (s.flux / solid_angle).to(units.intensity) + alt_unc_array = (s.uncertainty.array * s.flux.unit / solid_angle).to(units.intensity) + s = Spectrum1D(alt_flux, s.spectral_axis, uncertainty=StdDevUncertainty(alt_unc_array)) + return s def calculate_compounds(obsdata, pmodel): """ From 5ce6067528fc87d0111821bc2fd7b5ba0a9824a1 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 7 May 2024 14:56:40 -0400 Subject: [PATCH 27/68] Replace Model.plot() implementation by method that uses tabulate() --- pahfit/model.py | 222 +++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 208 insertions(+), 14 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index 69d22f3b..007bdbed 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -2,6 +2,7 @@ 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 @@ -11,9 +12,7 @@ from pahfit.base import PAHFITBase from pahfit import instrument from pahfit.errors import PAHFITModelError -from pahfit.component_models import BlackBody1D -from pahfit import units - +from pahfit.component_models import BlackBody1D, S07_attenuation class Model: """This class acts as the main API for PAHFIT. @@ -411,7 +410,14 @@ def info(self): """Print out the last fit results.""" print(self.astropy_result) - def plot(self, spec=None, redshift=None): + def plot( + self, + spec=None, + redshift=None, + use_instrument_fwhm=False, + label_lines=False, + **errorbar_kwargs, + ): """Plot model, and optionally compare to observational data. Parameters @@ -422,30 +428,218 @@ def plot(self, spec=None, redshift=None): redshift : float Redshift used to shift from the physical model, to the - observed model. + observed model. If None, it will be taken from spec.redshift - If None, will be taken from spec.redshift + use_instrument_fwhm : bool + For the lines, the default is to use the fwhm values + contained in the features table. When set to True, the fwhm + will be determined by the instrument model instead. + + label_lines : bool + Add labels with the names of the lines, at the position of each line. + + errorbar_kwargs : dict + Customize the data points plot by passing the given keyword + arguments to matplotlib.pyplot.errorbar. """ inst, z = self._parse_instrument_and_redshift(spec, redshift) _, _, _, xz, yz, uncz = self._convert_spec_data(spec, z) - # Always use the current FWHM here (use_instrument_fwhm would - # overwrite the fitted value in the instrument overlap regions!) - astropy_model = self._construct_astropy_model( - inst, z, use_instrument_fwhm=False + # total model + model = self._construct_astropy_model( + inst, z, use_instrument_fwhm=use_instrument_fwhm ) - - enough_samples = max(1000, len(spec.wavelength)) + enough_samples = max(10000, len(spec.wavelength)) + x_mod = np.logspace(np.log10(min(xz)), np.log10(max(xz)), enough_samples) fig, axs = plt.subplots( ncols=1, nrows=2, - figsize=(15, 10), + figsize=(10, 10), gridspec_kw={"height_ratios": [3, 1]}, sharex=True, ) - PAHFITBase.plot(axs, xz, yz, uncz, astropy_model, model_samples=enough_samples) + + # spectrum and best fit model + ax = axs[0] + ax.set_yscale("linear") + ax.set_xscale("log") + ax.minorticks_on() + ax.tick_params( + axis="both", which="major", top="on", right="on", direction="in", length=10 + ) + ax.tick_params( + axis="both", which="minor", top="on", right="on", direction="in", length=5 + ) + + ext_model = None + has_att = "attenuation" in self.features["kind"] + has_abs = "absorption" in self.features["kind"] + if has_att: + row = self.features[self.features["kind"] == "attenuation"][0] + tau = row["tau"][0] + ext_model = S07_attenuation(tau_sil=tau)(x_mod) + + if has_abs: + raise NotImplementedError( + "plotting absorption features not implemented yet" + ) + + if has_att or has_abs: + ax_att = ax.twinx() # axis for plotting the extinction curve + ax_att.tick_params(which="minor", direction="in", length=5) + ax_att.tick_params(which="major", direction="in", length=10) + ax_att.minorticks_on() + ax_att.plot(x_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)) + + # Define legend lines + Leg_lines = [ + mpl.lines.Line2D([0], [0], color="k", linestyle="--", lw=2), + mpl.lines.Line2D([0], [0], color="#FE6100", lw=2), + mpl.lines.Line2D([0], [0], color="#648FFF", lw=2, alpha=0.5), + mpl.lines.Line2D([0], [0], color="#DC267F", lw=2, alpha=0.5), + mpl.lines.Line2D([0], [0], color="#785EF0", lw=2, alpha=1), + mpl.lines.Line2D([0], [0], color="#FFB000", lw=2, alpha=0.5), + ] + + # local utility + def tabulate_components(kind): + ss = {} + for name in self.features[self.features["kind"] == kind]["name"]: + ss[name] = self.tabulate(inst, z, x_mod, self.features["name"] == name) + return {name: s.flux.value for name, s in ss.items()} + + cont_y = np.zeros(len(x_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) + # 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" + ).flux.value + ax.plot(x_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) + + # 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, + (cont_y + y) * ext_model, + "#648FFF", + alpha=0.5, + ) + + if "line" in self.features["kind"]: + for name, y in tabulate_components("line").items(): + ax.plot( + x_mod, + (cont_y + y) * ext_model, + "#DC267F", + alpha=0.5, + ) + if label_lines: + i = np.argmax(y) + # ignore out of range lines + if i > 0 and i < len(y) - 1: + w = x_mod[i] + ax.text( + w, + y[i], + name, + va="center", + ha="center", + rotation="vertical", + bbox=dict(facecolor="white", alpha=0.75, pad=0), + ) + + ax.plot(x_mod, self.tabulate(inst, z, x_mod).flux.value, "#FE6100", alpha=1) + + # data + default_kwargs = dict( + fmt="o", + markeredgecolor="k", + markerfacecolor="none", + ecolor="k", + elinewidth=0.2, + capsize=0.5, + markersize=6, + ) + + ax.errorbar(xz, yz, yerr=uncz, **(default_kwargs | errorbar_kwargs)) + + ax.set_ylim(0) + ax.set_ylabel(r"$\nu F_{\nu}$") + + ax.legend( + Leg_lines, + [ + "S07_attenuation", + "Spectrum Fit", + "Dust Features", + r"Atomic and $H_2$ Lines", + "Total Continuum Emissions", + "Continuum Components", + ], + prop={"size": 10}, + loc="best", + facecolor="white", + framealpha=1, + ncol=3, + ) + + # residuals, lower sub-figure + res = yz - model(xz) + std = np.nanstd(res) + ax = axs[1] + + ax.set_yscale("linear") + ax.set_xscale("log") + ax.tick_params( + axis="both", which="major", top="on", right="on", direction="in", length=10 + ) + ax.tick_params( + axis="both", which="minor", top="on", right="on", direction="in", length=5 + ) + ax.minorticks_on() + + # Custom X axis ticks + ax.xaxis.set_ticks( + [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 20, 25, 30, 40] + ) + + ax.axhline(0, linestyle="--", color="gray", zorder=0) + ax.plot( + xz, + res, + "ko-", + fillstyle="none", + zorder=1, + markersize=errorbar_kwargs.get("markersize", None), + alpha=errorbar_kwargs.get("alpha", None), + linestyle="none", + ) + scalefac_resid = 2 + ax.set_ylim(-scalefac_resid * std, scalefac_resid * std) + ax.set_xlim(np.floor(np.amin(xz)), np.ceil(np.amax(xz))) + ax.set_xlabel(r"$\lambda$ [$\mu m$]") + ax.set_ylabel("Residuals [%]") + + # scalar x-axis marks + ax.xaxis.set_major_formatter(mpl.ticker.ScalarFormatter()) fig.subplots_adjust(hspace=0) + return fig def copy(self): """Copy the model. From 1bbe11f4068b1faf679f19f35c42283be9f1ae98 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 7 May 2024 15:15:51 -0400 Subject: [PATCH 28/68] scalefac_resid argument for plot() and use it in plot_pahfit script --- pahfit/model.py | 9 +++++++-- pahfit/scripts/plot_pahfit.py | 35 +---------------------------------- 2 files changed, 8 insertions(+), 36 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index 007bdbed..5419597e 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -416,6 +416,7 @@ def plot( redshift=None, use_instrument_fwhm=False, label_lines=False, + scalefac_resid=2, **errorbar_kwargs, ): """Plot model, and optionally compare to observational data. @@ -436,7 +437,12 @@ def plot( will be determined by the instrument model instead. label_lines : bool - Add labels with the names of the lines, at the position of each line. + Add labels with the names of the lines, at the position of + each line. + + scalefac_resid : float + Factor multiplying the standard deviation of the residuals + to adjust plot limits. errorbar_kwargs : dict Customize the data points plot by passing the given keyword @@ -630,7 +636,6 @@ def tabulate_components(kind): alpha=errorbar_kwargs.get("alpha", None), linestyle="none", ) - scalefac_resid = 2 ax.set_ylim(-scalefac_resid * std, scalefac_resid * std) ax.set_xlim(np.floor(np.amin(xz)), np.ceil(np.amax(xz))) ax.set_xlabel(r"$\lambda$ [$\mu m$]") diff --git a/pahfit/scripts/plot_pahfit.py b/pahfit/scripts/plot_pahfit.py index d5bb8532..a8765058 100755 --- a/pahfit/scripts/plot_pahfit.py +++ b/pahfit/scripts/plot_pahfit.py @@ -6,11 +6,8 @@ import matplotlib as mpl from pahfit.model import Model -from pahfit.base import PAHFITBase from pahfit.helpers import read_spectrum -from astropy import units as u - def initialize_parser(): """ @@ -78,14 +75,6 @@ def main(): def default_layout_plot(spec, model, scalefac_resid): - """ - Returns - ------- - fig : Figure object - - """ - - # plot result fontsize = 18 font = {"size": fontsize} mpl.rc("font", **font) @@ -96,29 +85,7 @@ def default_layout_plot(spec, model, scalefac_resid): mpl.rc("xtick.minor", size=3, width=1) mpl.rc("ytick.minor", size=3, width=1) - fig, axs = plt.subplots( - ncols=1, - nrows=2, - figsize=(15, 10), - gridspec_kw={"height_ratios": [3, 1]}, - sharex=True, - ) - - enough_samples = max(1000, len(spec.wavelength)) - - PAHFITBase.plot( - axs, - spec.wavelength.to(u.micron).value, - spec.flux, - spec.uncertainty.array, - model._construct_astropy_model( - spec.meta["instrument"], spec.redshift, use_instrument_fwhm=False - ), - model_samples=enough_samples, - scalefac_resid=scalefac_resid, - ) - - # use the whitespace better + fig = model.plot(spec, scalefac_resid=scalefac_resid) fig.subplots_adjust(hspace=0) return fig From 99dda283c8b295ed853598bc44461a53f203af7d Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 7 May 2024 15:40:30 -0400 Subject: [PATCH 29/68] Remove old plot from PAHFITBase --- pahfit/base.py | 191 ------------------------------------------------- 1 file changed, 191 deletions(-) diff --git a/pahfit/base.py b/pahfit/base.py index a0ff47e9..0a55a881 100644 --- a/pahfit/base.py +++ b/pahfit/base.py @@ -1,5 +1,4 @@ import numpy as np -import matplotlib as mpl from pahfit.instrument import within_segment, fwhm from pahfit.errors import PAHFITModelError @@ -265,196 +264,6 @@ def model_from_param_info(param_info): return model - @staticmethod - def plot(axs, x, y, yerr, model, model_samples=1000, scalefac_resid=2): - """ - Plot model using axis object. - - Parameters - ---------- - axs : matplotlib.axis objects - where to put the plot - x : floats - wavelength points - y : floats - observed spectrum - yerr: floats - observed spectrum uncertainties - model : PAHFITBase model (astropy modeling CompoundModel) - model giving all the components and parameters - model_samples : int - Total number of wavelength points to allocate to the model display - scalefac_resid : float - Factor multiplying the standard deviation of the residuals to adjust plot limits - """ - # remove units if they are present - if hasattr(x, "value"): - x = x.value - if hasattr(y, "value"): - y = y.value - if hasattr(yerr, "value"): - yerr = yerr.value - - # Fine x samples for model fit - x_mod = np.logspace(np.log10(min(x)), np.log10(max(x)), model_samples) - - # spectrum and best fit model - ax = axs[0] - ax.set_yscale("linear") - ax.set_xscale("log") - ax.minorticks_on() - ax.tick_params(axis="both", - which="major", - top="on", - right="on", - direction="in", - length=10) - ax.tick_params(axis="both", - which="minor", - top="on", - right="on", - direction="in", - length=5) - - ax_att = ax.twinx() # axis for plotting the extinction curve - ax_att.tick_params(which="minor", direction="in", length=5) - ax_att.tick_params(which="major", direction="in", length=10) - ax_att.minorticks_on() - - # get the extinction model (probably a better way to do this) - ext_model = None - for cmodel in model: - if isinstance(cmodel, S07_attenuation): - ext_model = cmodel(x_mod) - - # get additional extinction components that can be - # characterized by functional forms (Drude profile in this case) - for cmodel in model: - if isinstance(cmodel, att_Drude1D): - if ext_model is not None: - ext_model *= cmodel(x_mod) - else: - ext_model = cmodel(x_mod) - ax_att.plot(x_mod, ext_model, "k--", alpha=0.5) - ax_att.set_ylabel("Attenuation") - ax_att.set_ylim(0, 1.1) - - # Define legend lines - Leg_lines = [ - mpl.lines.Line2D([0], [0], color="k", linestyle="--", lw=2), - mpl.lines.Line2D([0], [0], color="#FE6100", lw=2), - mpl.lines.Line2D([0], [0], color="#648FFF", lw=2, alpha=0.5), - mpl.lines.Line2D([0], [0], color="#DC267F", lw=2, alpha=0.5), - mpl.lines.Line2D([0], [0], color="#785EF0", lw=2, alpha=1), - mpl.lines.Line2D([0], [0], color="#FFB000", lw=2, alpha=0.5), - ] - - # create the continum compound model (base for plotting lines) - cont_components = [] - - for cmodel in model: - if isinstance(cmodel, BlackBody1D): - cont_components.append(cmodel) - # plot as we go - ax.plot(x_mod, - cmodel(x_mod) * ext_model / x_mod, - "#FFB000", - alpha=0.5) - cont_model = cont_components[0] - for cmodel in cont_components[1:]: - cont_model += cmodel - cont_y = cont_model(x_mod) - - # now plot the dust bands and lines - for cmodel in model: - if isinstance(cmodel, Gaussian1D): - ax.plot( - x_mod, - (cont_y + cmodel(x_mod)) * ext_model / x_mod, - "#DC267F", - alpha=0.5, - ) - if isinstance(cmodel, Drude1D): - ax.plot( - x_mod, - (cont_y + cmodel(x_mod)) * ext_model / x_mod, - "#648FFF", - alpha=0.5, - ) - - ax.plot(x_mod, cont_y * ext_model / x_mod, "#785EF0", alpha=1) - - ax.plot(x_mod, model(x_mod) / x_mod, "#FE6100", alpha=1) - ax.errorbar( - x, - y / x, - yerr=yerr / x, - fmt="o", - markeredgecolor="k", - markerfacecolor="none", - ecolor="k", - elinewidth=0.2, - capsize=0.5, - markersize=6, - ) - - ax.set_ylim(0) - ax.set_ylabel(r"$\nu F_{\nu}$") - - ax.legend( - Leg_lines, - [ - "S07_attenuation", - "Spectrum Fit", - "Dust Features", - r"Atomic and $H_2$ Lines", - "Total Continuum Emissions", - "Continuum Components", - ], - prop={"size": 10}, - loc="best", - facecolor="white", - framealpha=1, - ncol=3, - ) - - # residuals, lower sub-figure - res = (y - model(x)) / x - std = np.std(res) - ax = axs[1] - - ax.set_yscale("linear") - ax.set_xscale("log") - ax.tick_params(axis="both", - which="major", - top="on", - right="on", - direction="in", - length=10) - ax.tick_params(axis="both", - which="minor", - top="on", - right="on", - direction="in", - length=5) - ax.minorticks_on() - - # Custom X axis ticks - ax.xaxis.set_ticks( - [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 20, 25, 30, 40]) - - ax.axhline(0, linestyle="--", color="gray", zorder=0) - ax.plot(x, res, "ko-", fillstyle="none", zorder=1) - ax.set_ylim(-scalefac_resid * std, scalefac_resid * std) - ax.set_xlim(np.floor(np.amin(x)), np.ceil(np.amax(x))) - ax.set_xlabel(r"$\lambda$ [$\mu m$]") - ax.set_ylabel("Residuals [%]") - - # scalar x-axis marks - ax.xaxis.set_minor_formatter(mpl.ticker.ScalarFormatter()) - ax.xaxis.set_major_formatter(mpl.ticker.ScalarFormatter()) - - @staticmethod def update_dictionary(feature_dict, instrumentname, update_fwhms=False, redshift=0): """ Update parameter dictionary based on the instrument name. From 6a85e38f80cb02a1df66deb602a7b354bf56f871 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Mon, 13 May 2024 13:17:35 -0400 Subject: [PATCH 30/68] Fix accidentally removed pahfit.units import --- pahfit/model.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pahfit/model.py b/pahfit/model.py index 5419597e..851c845b 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -7,6 +7,7 @@ import numpy as np from scipy import interpolate, integrate +from pahfit import units from pahfit.features.util import bounded_is_fixed from pahfit.features import Features from pahfit.base import PAHFITBase From 50f232c6a776f5652087e2b04968ab71500473d5 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Mon, 13 May 2024 13:18:42 -0400 Subject: [PATCH 31/68] Tweaks to Model.plot() based on pull request #281 review --- pahfit/model.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/pahfit/model.py b/pahfit/model.py index 851c845b..1dd15a57 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -15,6 +15,7 @@ from pahfit.errors import PAHFITModelError from pahfit.component_models import BlackBody1D, S07_attenuation + class Model: """This class acts as the main API for PAHFIT. @@ -425,8 +426,10 @@ def plot( Parameters ---------- spec : Spectrum1D - Observational data. Does not have to be the same data that - was used for guessing or fitting. + Observational data. The units should be compatible with the + data that were used for the fit, but it does not have to be + the exact same spectrum. The spectrum will be converted to + internal units before plotting. redshift : float Redshift used to shift from the physical model, to the @@ -434,7 +437,7 @@ def plot( use_instrument_fwhm : bool For the lines, the default is to use the fwhm values - contained in the features table. When set to True, the fwhm + contained in the Features table. When set to True, the fwhm will be determined by the instrument model instead. label_lines : bool @@ -452,10 +455,6 @@ def plot( """ inst, z = self._parse_instrument_and_redshift(spec, redshift) _, _, _, xz, yz, uncz = self._convert_spec_data(spec, z) - # total model - model = self._construct_astropy_model( - inst, z, use_instrument_fwhm=use_instrument_fwhm - ) enough_samples = max(10000, len(spec.wavelength)) x_mod = np.logspace(np.log10(min(xz)), np.log10(max(xz)), enough_samples) @@ -606,8 +605,8 @@ def tabulate_components(kind): ncol=3, ) - # residuals, lower sub-figure - res = yz - model(xz) + # residuals = data in rest frame - (model evaluated at rest frame wavelengths) + res = yz - self.tabulate(inst, 0, xz).flux.value std = np.nanstd(res) ax = axs[1] @@ -630,7 +629,7 @@ def tabulate_components(kind): ax.plot( xz, res, - "ko-", + "ko", fillstyle="none", zorder=1, markersize=errorbar_kwargs.get("markersize", None), 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 32/68] 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 00000000..4a3d0510 --- /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 00000000..5c8b657e --- /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 33/68] 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 0a55a881..00000000 --- a/pahfit/base.py +++ /dev/null @@ -1,567 +0,0 @@ -import numpy as np - -from pahfit.instrument import within_segment, fwhm -from pahfit.errors import PAHFITModelError -from pahfit.component_models import ( - BlackBody1D, - ModifiedBlackBody1D, - S07_attenuation, - att_Drude1D, -) -from astropy.modeling.physical_models import Drude1D -from astropy.modeling.functional_models import Gaussian1D - -__all__ = ["PAHFITBase"] - - -def _ingest_limits(min_vals, max_vals): - """ - Ingest the limits read from yaml file and generate the appropriate - internal format (list of tuples). Needed to handle the case - where a limit is not desired as numpy arrays cannot have elements - of None, instead a value of nan is used. - - Limits that are not set are designated as 'nan' in files and - these are changed to the python None to be compatible with - the astropy.modeling convention. - - Parameters - ---------- - min_vals, - max_vals : numpy.array (masked arrays) - min/max values of the limits for a parameter - nan designates no limit - - Returns - ------- - plimits : list of tuples - tuples give the min/max limits for a parameter - """ - plimits = [] - mask_min = min_vals.mask - data_min = min_vals.data - mask_max = max_vals.mask - data_max = max_vals.data - - mask_min_ind = np.where(np.logical_not(mask_min))[0] - mask_max_ind = np.where(np.logical_not(mask_max))[0] - - min_vals = np.zeros(len(mask_min)) - min_vals[mask_min_ind] = data_min[mask_min_ind] - - max_vals = np.zeros(len(mask_max)) - max_vals[mask_max_ind] = data_max[mask_max_ind] - - plimits = [] - for cmin, cmax in zip(min_vals, max_vals): - if np.isinf(cmin): - cmin = None - if np.isinf(cmax): - cmax = None - plimits.append((cmin, cmax)) - - return plimits - - -def _ingest_fixed(fixed_vals): - """ - Ingest the fixed value read from a file and generate the appropriate - internal format (list of booleans). Since this information is indirectly - hidden in the parameter of a feature, this function is needed to - extract that. - - Parameters - ---------- - min_vals : numpy.array (masked array) - fixed designations - - Returns - ------- - pfixed : list (boolean) - True/False designation for parameters - """ - mask_false_ind = np.where(np.logical_not(fixed_vals.mask))[0] - fixed_vals = ["True"] * len(fixed_vals.mask) - for i in range(0, len(mask_false_ind)): - ind = mask_false_ind[i] - fixed_vals[ind] = "False" - - pfixed = [] - for cfixed in fixed_vals: - if cfixed == "True": - cfixed = True - if cfixed == "False": - cfixed = False - pfixed.append(cfixed) - - return pfixed - - -class PAHFITBase: - """ - Old implementation. Some functions are still used by the new Model - class. The unused functionality has been removed. - - Construct that is still used for now - - param_info: tuple of dicts called (bb_info, df_info, h2_info, ion_info, abs_info, att_info) - The dictionaries contain info for each type of component. Each - component of the dictonaries is a vector. - bb_info - - dict with {name, temps, temps_limits, temps_fixed, - amps, amps_limits, amps_fixed}, each a vector - dust_features, h2_features, ion_features - - dict with {name amps, amps_limits, amps_fixed, - x_0, x_0_limits, x_0_fixed, fwhms, fwhms_limits, fwhm_fixed}. - """ - @staticmethod - def model_from_param_info(param_info): - # setup the model - bb_info = param_info[0] - dust_features = param_info[1] - h2_features = param_info[2] - ion_features = param_info[3] - abs_info = param_info[4] - att_info = param_info[5] - - model = None - if bb_info is not None: - bbs = [] - for k in range(len(bb_info["names"])): - BBClass = ModifiedBlackBody1D if bb_info["modified"][ - k] else BlackBody1D - bbs.append( - BBClass( - name=bb_info["names"][k], - temperature=bb_info["temps"][k], - amplitude=bb_info["amps"][k], - bounds={ - "temperature": bb_info["temps_limits"][k], - "amplitude": bb_info["amps_limits"][k], - }, - fixed={ - "temperature": bb_info["temps_fixed"][k], - "amplitude": bb_info["amps_fixed"][k], - }, - ) - ) - model = sum(bbs[1:], bbs[0]) - - if dust_features is not None: - df = [] - for k in range(len(dust_features["names"])): - df.append( - Drude1D( - name=dust_features["names"][k], - amplitude=dust_features["amps"][k], - x_0=dust_features["x_0"][k], - fwhm=dust_features["fwhms"][k], - bounds={ - "amplitude": dust_features["amps_limits"][k], - "x_0": dust_features["x_0_limits"][k], - "fwhm": dust_features["fwhms_limits"][k], - }, - fixed={ - "amplitude": dust_features["amps_fixed"][k], - "x_0": dust_features["x_0_fixed"][k], - "fwhm": dust_features["fwhms_fixed"][k], - }, - )) - - df = sum(df[1:], df[0]) - if model: - model += df - else: - model = df - - if h2_features is not None: - h2 = [] - for k in range(len(h2_features["names"])): - h2.append( - Gaussian1D( - name=h2_features["names"][k], - amplitude=h2_features["amps"][k], - mean=h2_features["x_0"][k], - stddev=h2_features["fwhms"][k] / 2.355, - bounds={ - "amplitude": - h2_features["amps_limits"][k], - "mean": - h2_features["x_0_limits"][k], - "stddev": ( - h2_features["fwhms"][k] * 0.9 / 2.355, - h2_features["fwhms"][k] * 1.1 / 2.355, - ), - }, - fixed={ - "amplitude": h2_features["amps_fixed"][k], - "mean": h2_features["x_0_fixed"][k], - "stddev": h2_features["fwhms_fixed"][k], - }, - )) - h2 = sum(h2[1:], h2[0]) - if model: - model += h2 - else: - model = h2 - - if ion_features is not None: - ions = [] - for k in range(len(ion_features["names"])): - ions.append( - Gaussian1D( - name=ion_features["names"][k], - amplitude=ion_features["amps"][k], - mean=ion_features["x_0"][k], - stddev=ion_features["fwhms"][k] / 2.355, - bounds={ - "amplitude": - ion_features["amps_limits"][k], - "mean": - ion_features["x_0_limits"][k], - "stddev": ( - ion_features["fwhms"][k] * 0.9 / 2.355, - ion_features["fwhms"][k] * 1.1 / 2.355, - ), - }, - fixed={ - "amplitude": ion_features["amps_fixed"][k], - "mean": ion_features["x_0_fixed"][k], - "stddev": ion_features["fwhms_fixed"][k], - }, - )) - ions = sum(ions[1:], ions[0]) - if model: - model += ions - else: - model = ions - - # add additional att components to the model if necessary - if not model: - raise PAHFITModelError("No model components found") - - if abs_info is not None: - for k in range(len(abs_info["names"])): - model *= att_Drude1D( - name=abs_info["names"][k], - tau=abs_info["amps"][k], - x_0=abs_info["x_0"][k], - fwhm=abs_info["fwhms"][k], - bounds={ - "tau": abs_info["amps_limits"][k], - "fwhm": abs_info["fwhms_limits"][k], - }, - fixed={"x_0": abs_info["x_0_fixed"][k]}, - ) - - if att_info is not None: - model *= S07_attenuation( - name=att_info["name"], - tau_sil=att_info["tau_sil"], - bounds={"tau_sil": att_info["tau_sil_limits"]}, - fixed={"tau_sil": att_info["tau_sil_fixed"]}, - ) - - return model - - def update_dictionary(feature_dict, instrumentname, update_fwhms=False, redshift=0): - """ - Update parameter dictionary based on the instrument name. - Based on the instrument name, this function removes the - features outside of the wavelength range and - updates the FWHMs of the lines. - - - Parameters - ---------- - feature_dict : dictionary - Dictionary created by reading in a science pack. - - instrumentname : string - Name of the instrument with which the input spectrum - is observed. - - update_fwhms = Boolean - True for h2_info and ion_info - False for df_info - - Returns - ------- - updated feature_dict - """ - if feature_dict is None: - return None - - # convert from physical feature, to observed wavelength - def redshifted_waves(): - return feature_dict["x_0"] * (1 + redshift) - - ind = np.nonzero(within_segment(redshifted_waves(), instrumentname))[0] - - # select the valid entries in these arrays - array_keys = ("x_0", "amps", "fwhms", "names") - new_values_1 = {key: feature_dict[key][ind] for key in array_keys} - - # these are lists instead - list_keys = ( - "amps_fixed", - "fwhms_fixed", - "x_0_fixed", - "x_0_limits", - "amps_limits", - "fwhms_limits", - ) - new_values_2 = { - key: [feature_dict[key][i] for i in ind] for key in list_keys - } - - feature_dict.update(new_values_1) - feature_dict.update(new_values_2) - - if len(feature_dict['names']) == 0: - # if we removed all the things, be careful here. Setting to - # None should make the model construction function behave - # normally. - feature_dict = None - return feature_dict - - if update_fwhms: - # observe the lines at the redshifted wavelength - fwhm_min_max = fwhm(instrumentname, redshifted_waves(), as_bounded=True) - # shift the observed fwhm back to the rest frame (where the - # observed data will be moved, and its features will become - # narrower) - fwhm_min_max /= (1 + redshift) - # For astropy a numpy.bool does not work for the 'fixed' - # parameter. It needs to be a regular bool. Doing tolist() - # instead of using the array mask directly solves this. - feature_dict.update( - { - "fwhms": fwhm_min_max[:, 0], - # masked means there is no min/max, i.e. they need to be fixed - "fwhms_fixed": fwhm_min_max[:, 1].mask.tolist(), - "fwhms_limits": fwhm_min_max[:, 1:].tolist(), - } - ) - - return feature_dict - - @staticmethod - def parse_table(pack_table): - """ - Load the model parameters from a Table - - Parameters - ---------- - pack_table : Table - Table created by reading in a science pack. - - Returns - ------- - readout : tuple - Tuple containing dictionaries of all components from the - input Table. Can be used to create PAHFITBase instance using - param_info argument. Dictionary in tuple is None if no - components of that type were specified. - """ - # Getting indices for the different components - bb_ind = np.where((pack_table["kind"] == "starlight") - | (pack_table["kind"] == "dust_continuum"))[0] - df_ind = np.where(pack_table["kind"] == "dust_feature")[0] - ga_ind = np.where(pack_table["kind"] == "line")[0] - at_ind = np.where(pack_table["kind"] == "attenuation")[0] - ab_ind = np.where(pack_table["kind"] == "absorption")[0] - - # now split the gas emission lines between H2 and ions - names = [str(i) for i in pack_table["name"][ga_ind]] - if len(names) > 0: - # this has trouble with empty list - h2_temp = np.char.find(names, "H2") >= 0 - ion_temp = np.char.find(names, "H2") == -1 - h2_ind = ga_ind[h2_temp] - ion_ind = ga_ind[ion_temp] - else: - h2_ind = [] - ion_ind = [] - # the rest works fine with empty list - - # Creating the blackbody dict - bb_info = None - if len(bb_ind) > 0: - bb_info = { - "names": - np.array(pack_table["name"][bb_ind].data), - "temps": - np.array(pack_table["temperature"][:, 0][bb_ind].data), - "temps_limits": - _ingest_limits( - pack_table["temperature"][:, 1][bb_ind], - pack_table["temperature"][:, 2][bb_ind], - ), - "temps_fixed": - _ingest_fixed(pack_table["temperature"][:, 1][bb_ind]), - "amps": - np.array(pack_table["tau"][:, 0][bb_ind].data), - "amps_limits": - _ingest_limits( - pack_table["tau"][:, 1][bb_ind], - pack_table["tau"][:, 2][bb_ind], - ), - "amps_fixed": - _ingest_fixed(pack_table["tau"][:, 1][bb_ind]), - "modified": - np.array(pack_table["kind"][bb_ind] == "dust_continuum"), - } - - # Creating the dust_features dict - df_info = None - if len(df_ind) > 0: - df_info = { - "names": - np.array(pack_table["name"][df_ind].data), - "x_0": - np.array(pack_table["wavelength"][:, 0][df_ind].data), - "x_0_limits": - _ingest_limits( - pack_table["wavelength"][:, 1][df_ind], - pack_table["wavelength"][:, 2][df_ind], - ), - "x_0_fixed": - _ingest_fixed(pack_table["wavelength"][:, 1][df_ind]), - "amps": - np.array(pack_table["power"][:, 0][df_ind].data), - "amps_limits": - _ingest_limits( - pack_table["power"][:, 1][df_ind], - pack_table["power"][:, 2][df_ind], - ), - "amps_fixed": - _ingest_fixed(pack_table["power"][:, 1][df_ind]), - "fwhms": - np.array(pack_table["fwhm"][:, 0][df_ind].data), - "fwhms_limits": - _ingest_limits( - pack_table["fwhm"][:, 1][df_ind], - pack_table["fwhm"][:, 2][df_ind], - ), - "fwhms_fixed": - _ingest_fixed(pack_table["fwhm"][:, 1][df_ind]), - } - - # Creating the H2 dict - h2_info = None - if len(h2_ind) > 0: - h2_info = { - "names": - np.array(pack_table["name"][h2_ind].data), - "x_0": - np.array(pack_table["wavelength"][:, 0][h2_ind].data), - "x_0_limits": - _ingest_limits( - pack_table["wavelength"][:, 1][h2_ind], - pack_table["wavelength"][:, 2][h2_ind], - ), - "x_0_fixed": - _ingest_fixed(pack_table["wavelength"][:, 1][h2_ind]), - "amps": - np.array(pack_table["power"][:, 0][h2_ind].data), - "amps_limits": - _ingest_limits( - pack_table["power"][:, 1][h2_ind], - pack_table["power"][:, 2][h2_ind], - ), - "amps_fixed": - _ingest_fixed(pack_table["power"][:, 1][h2_ind]), - "fwhms": - np.array(pack_table["fwhm"][:, 0][h2_ind].data), - "fwhms_limits": - _ingest_limits( - pack_table["fwhm"][:, 1][h2_ind], - pack_table["fwhm"][:, 2][h2_ind], - ), - "fwhms_fixed": - _ingest_fixed(pack_table["fwhm"][:, 1][h2_ind]), - } - - # Creating the ion dict - ion_info = None - if len(ion_ind) > 0: - ion_info = { - "names": - np.array(pack_table["name"][ion_ind].data), - "x_0": - np.array(pack_table["wavelength"][:, 0][ion_ind].data), - "x_0_limits": - _ingest_limits( - pack_table["wavelength"][:, 1][ion_ind], - pack_table["wavelength"][:, 2][ion_ind], - ), - "x_0_fixed": - _ingest_fixed(pack_table["wavelength"][:, 1][ion_ind]), - "amps": - np.array(pack_table["power"][:, 0][ion_ind].data), - "amps_limits": - _ingest_limits( - pack_table["power"][:, 1][ion_ind], - pack_table["power"][:, 2][ion_ind], - ), - "amps_fixed": - _ingest_fixed(pack_table["power"][:, 1][ion_ind]), - "fwhms": - np.array(pack_table["fwhm"][:, 0][ion_ind].data), - "fwhms_limits": - _ingest_limits( - pack_table["fwhm"][:, 1][ion_ind].data, - pack_table["fwhm"][:, 2][ion_ind].data, - ), - "fwhms_fixed": - _ingest_fixed(pack_table["fwhm"][:, 1][ion_ind].data), - } - - # Create the attenuation dict (could be absorption drudes - # and S07 model) - abs_info = None - if len(ab_ind) > 0: - abs_info = { - "names": - np.array(pack_table["name"][at_ind].data), - "x_0": - np.array(pack_table["wavelength"][:, 0][at_ind].data), - "x_0_limits": - _ingest_limits( - pack_table["wavelength"][:, 1][at_ind], - pack_table["wavelength"][:, 2][at_ind], - ), - "x_0_fixed": - _ingest_fixed(pack_table["wavelength"][:, 1][at_ind]), - "amps": - np.array(pack_table["tau"][:, 0][at_ind].data), - "amps_limits": - _ingest_limits( - pack_table["tau"][:, 0][at_ind], - pack_table["tau"][:, 1][at_ind], - ), - "amps_fixed": - _ingest_fixed(pack_table["tau"][:, 1][at_ind]), - "fwhms": - np.array(pack_table["fwhm"][:, 0][at_ind].data), - "fwhms_limits": - _ingest_limits( - pack_table["fwhm"][:, 1][at_ind], - pack_table["fwhm"][:, 2][at_ind], - ), - "fwhms_fixed": - _ingest_fixed(pack_table["fwhm"][:, 1][at_ind]), - } - - att_info = None - if len(at_ind) > 1: - raise NotImplementedError("More than one attenuation component not supported") - elif len(at_ind) == 1: - i = at_ind[0] - att_info = {"name": pack_table["name"][i], - "tau_sil": pack_table["tau"][i][0], - "tau_sil_limits": pack_table["tau"][i][1:], - "tau_sil_fixed": True if pack_table["tau"][i].mask[1] else False} - - return [bb_info, df_info, h2_info, ion_info, abs_info, att_info] diff --git a/pahfit/model.py b/pahfit/model.py index 1dd15a57..5b772046 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 ee0e39e5..cc7c6b9d 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 d56975bf..56e4aa56 100644 --- a/pahfit/tests/test_fitting_spitzer.py +++ b/pahfit/tests/test_fitting_spitzer.py @@ -66,7 +66,7 @@ def test_fitting_m101(): try: np.testing.assert_allclose( - model.astropy_result.parameters, expvals, rtol=1e-6, atol=1e-6 + model.fitter.model, expvals, rtol=1e-6, atol=1e-6 ) except AssertionError as error: print(error) diff --git a/pahfit/tests/test_model_impl.py b/pahfit/tests/test_model_impl.py index fce3acc3..c15a429a 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 34/68] 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 4bb56f7c..31da3131 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 35/68] 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 5b772046..a5a1a76f 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 c15a429a..72a1aca8 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 36/68] 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 31da3131..0efed64b 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -137,7 +137,7 @@ class Features(Table): _group_attrs = set(('bounds', 'features', 'kind')) # group-level attributes _param_attrs = set(('value', 'bounds')) # Each parameter can have these attributes _no_bounds = set(('name', 'group', 'kind', 'geometry', 'model')) # str attributes (no bounds) - _bounds_dtype = np.dtype([("val", 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 fc84ee51..f6bcabc2 100644 --- a/pahfit/helpers.py +++ b/pahfit/helpers.py @@ -1,16 +1,13 @@ import os from importlib import resources -from pahfit.component_models import BlackBody1D, S07_attenuation from pahfit import units from specutils import Spectrum1D -from astropy.modeling.physical_models import Drude1D -from astropy.modeling.functional_models import Gaussian1D from astropy import units as u from astropy.nddata import StdDevUncertainty -__all__ = ["read_spectrum", "calculate_compounds"] +__all__ = ["find_packfile", "read_spectrum"] def find_packfile(packfile): @@ -90,108 +87,13 @@ def read_spectrum(specfile, format=None): # for now. To be removed when dual unit support (intensity and flux # density) is supported. if s.flux.unit.is_equivalent(units.flux_density): - solid_angle = (3 * u.arcsec)**2 + solid_angle = (3 * u.arcsec) ** 2 alt_flux = (s.flux / solid_angle).to(units.intensity) - alt_unc_array = (s.uncertainty.array * s.flux.unit / solid_angle).to(units.intensity) - s = Spectrum1D(alt_flux, s.spectral_axis, uncertainty=StdDevUncertainty(alt_unc_array)) + alt_unc_array = (s.uncertainty.array * s.flux.unit / solid_angle).to( + units.intensity + ) + s = Spectrum1D( + alt_flux, s.spectral_axis, uncertainty=StdDevUncertainty(alt_unc_array) + ) return s - -def calculate_compounds(obsdata, pmodel): - """ - Determine model compounds for total continuum, stellar continuum, - total dust continuum, combined dust features, - combined atomic and H2 lines, combined H2 lines, - combined atomic lines, and extinction model - - Parameters - ---------- - obsdata : dict - observed data where x = wavelength, y = SED, and unc = uncertainties - - pmodel : PAHFITBase model - model giving all the components and parameters - - Returns - ------- - compounds : dict - x = wavelength in microns; - tot_cont = total continuum; - stellar_cont = stellar continuum; - dust_cont = total dust continuum; - dust_features = combined dust features; - tot_lines = combined atomic and H2 emission lines; - h2_lines = combined H2 lines; - atomic_lines = combined atomic lines; - extinction_model = extinction model - """ - - # get wavelength array - x = obsdata["x"].value - - # calculate total dust continuum and total continuum (including stellar continuum) - # v2.0: first BlackBody1D is stellar continuum - cont_components = [] - - for cmodel in pmodel.model: - if isinstance(cmodel, BlackBody1D): - cont_components.append(cmodel) - stellar_cont_model = cont_components[0] - dust_cont_model = cont_components[1] - for cmodel in cont_components[2:]: - dust_cont_model += cmodel - totcont = dust_cont_model(x) + stellar_cont_model(x) - - # calculate total dust features - dust_features = [] - - for cmodel in pmodel.model: - if isinstance(cmodel, Drude1D): - dust_features.append(cmodel) - dust_features_model = dust_features[0] - for cmodel in dust_features[1:]: - dust_features_model += cmodel - - # calculate H2 spectrum - h2_features = [] - - for cmodel in pmodel.model: - if isinstance(cmodel, Gaussian1D): - if cmodel.name[0:2] == "H2": - h2_features.append(cmodel) - h2_features_model = h2_features[0] - for cmodel in h2_features[1:]: - h2_features_model += cmodel - - # calculate atomic line spectrum - atomic_features = [] - - for cmodel in pmodel.model: - if isinstance(cmodel, Gaussian1D): - if cmodel.name[0:2] != "H2": - atomic_features.append(cmodel) - atomic_features_model = atomic_features[0] - for cmodel in atomic_features[1:]: - atomic_features_model += cmodel - - # all atomic and H2 lines - totlines = h2_features_model(x) + atomic_features_model(x) - - # get extinction model - for cmodel in pmodel.model: - if isinstance(cmodel, S07_attenuation): - ext_model = cmodel(x) - - # save compounds in dictionary - compounds = {} - compounds["x"] = x - compounds["tot_cont"] = totcont - compounds["stellar_cont"] = stellar_cont_model(x) - compounds["dust_cont"] = dust_cont_model(x) - compounds["dust_features"] = dust_features_model(x) - compounds["tot_lines"] = totlines - compounds["h2_lines"] = h2_features_model(x) - compounds["atomic_lines"] = atomic_features_model(x) - compounds["extinction_model"] = ext_model - - return compounds 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 37/68] 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 21352c0d..905cf99c 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 cc7c6b9d..9791e4db 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 38/68] 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 5c8b657e..f009da52 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 39/68] 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 4a3d0510..37f14125 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 f009da52..877b09d1 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 a5a1a76f..9e77e656 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 40/68] 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 37f14125..0ce006c0 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 877b09d1..3ec5052a 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 41/68] 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 0ce006c0..0e5035e1 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 3ec5052a..f81b9590 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 9e77e656..3a2493a8 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 42/68] 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 0e5035e1..22b1db02 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 f81b9590..ad474a6e 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 43/68] 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 22b1db02..14a34680 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 ad474a6e..df4f193a 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 3a2493a8..004c3e7e 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 44/68] 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 df4f193a..e3eda0bc 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 45/68] 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 004c3e7e..60a059d9 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 46/68] 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 60a059d9..9df50848 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 9791e4db..d658a689 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 47/68] 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 14a34680..9f11f332 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 e3eda0bc..a75f2bde 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 9df50848..bbe46d2f 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 72a1aca8..6007c175 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 48/68] 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 bbe46d2f..034c9d35 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 49/68] 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 034c9d35..212e223d 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 50/68] 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 212e223d..53c85181 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 51/68] 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 53c85181..e16ed8d9 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 52/68] 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 e16ed8d9..52587354 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 53/68] 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 905cf99c..97952a00 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 f1c67228..40eba460 100644 --- a/pahfit/feature_strengths.py +++ b/pahfit/feature_strengths.py @@ -1,5 +1,5 @@ from astropy.modeling.functional_models import Gaussian1D -from pahfit.component_models import BlackBody1D +from pahfit.fitters.ap_components import BlackBody1D import numpy as np diff --git a/pahfit/fitters/__init__.py b/pahfit/fitters/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/pahfit/component_models.py b/pahfit/fitters/ap_components.py similarity index 100% rename from pahfit/component_models.py rename to pahfit/fitters/ap_components.py diff --git a/pahfit/apfitter.py b/pahfit/fitters/ap_fitter.py similarity index 99% rename from pahfit/apfitter.py rename to pahfit/fitters/ap_fitter.py index 9f11f332..92760bf0 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 52587354..e6d4c4ef 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 54/68] 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 97952a00..7073faff 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 55/68] 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 7d75f0a8..c94f28eb 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -12,4 +12,4 @@ ================== - IDL versions -- http://tir.astro.utoledo.edu/jdsmith/research/pahfit.php +- https://github.com/PAHFIT/pahfit_classic diff --git a/README.rst b/README.rst index dc92abed..1a02367c 100644 --- a/README.rst +++ b/README.rst @@ -2,7 +2,7 @@ PAHFIT ====== -``PAHFIT`` is a decomposition model and tool for astronomical infrared spectra, focusing on dust and gas emission features from the interstellar medium (see `Smith, J.D.T., Draine B.T., et al., 2007, ApJ, 656, 770 `_). +``PAHFIT`` is a decomposition model and tool for astronomical infrared spectra, focusing on dust and gas emission features from the interstellar medium (see `Smith, J.D.T., Draine B.T., et al., 2007, ApJ, 656, 770 `_). This package provides an updated python implementation of ``PAHFIT``. While the original versions of ``PAHFIT`` (``v1.x``) were written in IDL and focused mainly on Spitzer/IRS spectroscopic observations, the newer python-based versions (``>=v2.0``) will expand instrument coverage to other existing (e.g., AKARI) and planned (e.g., JWST) facilities, and will offer a more flexible modeling framework suitable for modeling a wider range of astrophysical sources. diff --git a/docs/index.rst b/docs/index.rst index 7073faff..e6072367 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -14,7 +14,7 @@ and a more flexible modeling framework suitable for modeling a wider range of astrophysical sources. For details for the IDL version of PAHFIT see -`Smith, J.D.T., Draine B.T., et al., 2007, ApJ, 656, 770 `_. +`Smith, J.D.T., Draine B.T., et al., 2007, ApJ, 656, 770 `_. This package is potentially an `astropy affiliated package `_ diff --git a/docs/version_description/version_description.rst b/docs/version_description/version_description.rst index 609592da..882f12a5 100644 --- a/docs/version_description/version_description.rst +++ b/docs/version_description/version_description.rst @@ -5,7 +5,7 @@ Version Description v1.0 - v1.4 ------------ -The original IDL version of PAHFIT can be found in: `http://tir.astro.utoledo.edu/jdsmith/research/pahfit.php `_. For more details and background of ``PAHFIT``, see the `background `_ page. +The original IDL version of PAHFIT can be found in: `https://github.com/PAHFIT/pahfit_classic `_. For more details and background of ``PAHFIT``, see the `background `_ page. v2.0 ------------ 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 56/68] 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 e6d4c4ef..e2ef45e1 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 57/68] 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 e2ef45e1..22935809 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 58/68] 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 22935809..83380b07 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( From f80c4605cd5b2827d33bb0c4a9bff44e147ddd77 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 25 Jun 2024 18:51:00 -0400 Subject: [PATCH 59/68] Add PowerDrude1D and PowerGaussian1D classes --- pahfit/fitters/ap_components.py | 145 ++++++++++++++++++++++++++++++++ 1 file changed, 145 insertions(+) diff --git a/pahfit/fitters/ap_components.py b/pahfit/fitters/ap_components.py index db719492..7604a61a 100644 --- a/pahfit/fitters/ap_components.py +++ b/pahfit/fitters/ap_components.py @@ -3,6 +3,8 @@ from astropy.modeling.physical_models import Drude1D from astropy.modeling import Fittable1DModel from astropy.modeling import Parameter +from astropy import constants +from pahfit import units __all__ = ["BlackBody1D", "ModifiedBlackBody1D", "S07_attenuation", "att_Drude1D"] @@ -134,3 +136,146 @@ def evaluate(x, tau, x_0, fwhm): profile = Drude1D(amplitude=1.0, fwhm=fwhm, x_0=x_0) tau_x = tau * profile(x) return (1.0 - np.exp(-1.0 * tau_x)) / tau_x + + +class PowerDrude1D(Fittable1DModel): + """Drude profile with amplitude determined by power. + + Special attention is needed for the units. An implementation for + this is 'unitful' because 'power' is defined as integral over + frequency, while the profile is formulated as a function of + wavelength. If we assume the output has unit(flux), then without + additional conversions the power unit will be unit(flux) * frequency + = unit(flux) * unit(c) / unit(lam). + + Example: if x_0 and fwhm are in micron, and the flux is in MJy / sr, + then the unit of the fitted power will be MJy sr-1 * unit(c) / + micron, which differs by a constant factor from MJy sr-1 Hz, + depending on the chosen unit of c. + + For efficiency and to prevent ambiguity, we assume that all units + are the internal pahfit units defined in pahfit.units, and + precalculate a conversion factor. + + TODO: We need to check if the flux is 'intensity' or 'flux_density', + and assume the power parameter has 'intensity_power' or + 'flux_density_power' units respectively. For now, only intensity is + supported. + + """ + + power = Parameter(min=0.0) + x_0 = Parameter(min=0.0) + fwhm = Parameter(default=1, min=0.0) + + # constant factors in the equation to convert power to amplitude of + # the profile. + intensity_amplitude_factor = ( + (2 * units.intensity_power * units.wavelength / (constants.c * np.pi)) + .to(units.intensity) + .value + ) + + @staticmethod + def evaluate(x, power, x_0, fwhm): + """Smith, et al. (2007) dust features model. Calculation is for + a Drude profile (equation in section 4.1.4). + + The intensity profile as a function of + wavelength is + + Inu(lambda) = (b * g**2) / ((lambda / x0 - x0 / lambda)**2 + g**2) + + With + b = amplitude (has same unit as flux) + g = fwhm / x0 + x0 = central wavelength + + The integrated power (Fnu dnu) of the profile is + + P = (pi * c / 2) * (b * g / x0) + + Which can be solved for the amplitude b. + + b = (P * 2 * x0) / (pi * c * g) = 2P / (pi nu0 g). + + According to the above equations, without additional + conversions, the resulting amplitude unit will be unit(P) * + Hz-1. This will result in very small values for for Inu(lambda), + or very large values for P. To avoid numerical problems with the + fitting algorithm, we apply conversions so that Inu(lambda) and + P are in the internal units. + + Parameters + ---------- + power : float + fwhm : float + central intensity (x_0) : float + + """ + # The equation and unit conversion for the amplitude: + # b = (2 * P * x_0 / (pi * c * g)).to(output_unit).value + + # Use predetermined factor that deals with units and constants. + # factor = (2 * unit(power) * unit(wavelength) / (pi * c)).to(unit(intensity)) + + g = fwhm / x_0 + b = power * x_0 / g * PowerDrude1D.intensity_amplitude_factor + return b * g**2 / ((x / x_0 - x_0 / x) ** 2 + g**2) + + +class PowerGaussian1D(Fittable1DModel): + """Gaussian profile with amplitude derived from power. + + Implementation analogous to PowerDrude1D. + + The amplitude of a gaussian profile given its power P, is P / + (stddev sqrt(2 pi)). Since stddev is given in wavelength units, this + equation gives the peak density per wavelength interval, Alambda. + The profile Flambda is then + + Flambda(lambda) = Alambda * G(lambda; mean, stddev) + + where G is a gaussian profile with amplitude 1.Converting this to + Fnu units yields + + Fnu(lambda) = lambda**2 / c * Flambda(lambda) + + Approximating the lambda**2 factor as a constant (the central + wavelength = 'mean'), then yields + + Fnu(lambda) = mean**2 / c * Alambda * G(lambda; mean, stddev) + + In other words, for narrow lines, the per-frequency profile is + approximately a gaussian with amplitude + + Anu = P * mean**2 / (c * stddev sqrt(2 pi)). + + So the constant factor we can set is + (unit(power) * unit(wavelength)**2 / (c * unit(wavelength) * sqrt(2 pi))).to(intensity) + + """ + + power = Parameter(min=0.0) + mean = Parameter() + stddev = Parameter(default=1, min=0.0) + + intensity_amplitude_factor = ( + ( + units.intensity_power + * (units.wavelength) ** 2 + / (constants.c * units.wavelength * np.sqrt(2 * np.pi)) + ) + .to(units.intensity) + .value + ) + + @staticmethod + def evaluate(x, power, mean, stddev): + """Evaluate F_nu(lambda) given the power. + + See class description for equations and unit notes.""" + + # amplitude in intensity units + Anu = power * mean**2 / stddev * PowerGaussian1D.intensity_amplitude_factor + return Anu * np.exp(-0.5 * np.square((x - mean) / stddev)) From 76d1d2033af14ad067c91e14ea055dddbfb093c4 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 25 Jun 2024 19:40:48 -0400 Subject: [PATCH 60/68] Switch from amplitude to power units and adjust guess accordingly --- pahfit/fitters/ap_fitter.py | 23 +++++----- pahfit/model.py | 89 ++++++++++++++++++++++--------------- 2 files changed, 64 insertions(+), 48 deletions(-) diff --git a/pahfit/fitters/ap_fitter.py b/pahfit/fitters/ap_fitter.py index 92760bf0..af189341 100644 --- a/pahfit/fitters/ap_fitter.py +++ b/pahfit/fitters/ap_fitter.py @@ -5,9 +5,9 @@ ModifiedBlackBody1D, S07_attenuation, att_Drude1D, - Drude1D, + PowerDrude1D, + PowerGaussian1D ) -from astropy.modeling.functional_models import Gaussian1D from astropy.modeling.fitting import LevMarLSQFitter import numpy as np @@ -161,10 +161,10 @@ def add_feature_line(self, name, power, wavelength, fwhm): kwargs = self._astropy_model_kwargs( name, - ["amplitude", "mean", "stddev"], + ["power", "mean", "stddev"], [power, wavelength, fwhm / 2.355], ) - self._add_component(Gaussian1D, **kwargs) + self._add_component(PowerGaussian1D, **kwargs) def add_feature_dust_feature(self, name, power, wavelength, fwhm): """Register a PowerDrude1D. @@ -175,9 +175,9 @@ def add_feature_dust_feature(self, name, power, wavelength, fwhm): """ self.feature_types[name] = "dust_feature" kwargs = self._astropy_model_kwargs( - name, ["amplitude", "x_0", "fwhm"], [power, wavelength, fwhm] + name, ["power", "x_0", "fwhm"], [power, wavelength, fwhm] ) - self._add_component(Drude1D, **kwargs) + self._add_component(PowerDrude1D, **kwargs) def add_feature_attenuation(self, name, tau, model="S07", geometry="screen"): """Register the S07 attenuation component. @@ -286,10 +286,6 @@ def get_result(self, component_name): (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 @@ -300,7 +296,8 @@ def get_result(self, component_name): ------- dict with Parameters according to the PAHFIT definitions. - e.g. {'power': converted from amplitude, 'fwhm': converted from + e.g., for a feature with amplitude, stddev, and mean parameters: + {'power': converted from amplitude, 'fwhm': converted from stddev, 'mean': wavelength} """ @@ -322,13 +319,13 @@ def get_result(self, component_name): } elif c_type == "line": return { - "power": component.amplitude.value, + "power": component.power.value, "wavelength": component.mean.value, "fwhm": component.stddev.value * 2.355, } elif c_type == "dust_feature": return { - "power": component.amplitude.value, + "power": component.power.value, "wavelength": component.x_0.value, "fwhm": component.fwhm.value, } diff --git a/pahfit/model.py b/pahfit/model.py index 83380b07..347b5538 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -1,5 +1,6 @@ from specutils import Spectrum1D from astropy import units as u +from astropy import constants import copy import matplotlib as mpl from matplotlib import pyplot as plt @@ -187,8 +188,13 @@ def guess( # parse spectral data self.features.meta["user_unit"]["flux"] = spec.flux.unit _, _, _, lam, flux, _ = self._convert_spec_data(spec, z) - wmin = min(lam) - wmax = max(lam) + lam_min = min(lam) + lam_max = max(lam) + + # Some useful quantities for guessing + median_flux = np.median(flux) + Flambda = flux * units.intensity * (lam * units.wavelength) ** -2 * constants.c + total_power = np.trapz(Flambda, lam * units.wavelength) # simple linear interpolation function for spectrum sp = interpolate.interp1d(lam, flux) @@ -205,7 +211,7 @@ def loop_over_non_fixed(kind, parameter, estimate_function, force=False): # guess starting point of bb def starlight_guess(row): bb = BlackBody1D(1, row["temperature"][0]) - w = wmin + 0.1 # the wavelength used to compare + w = lam_min + 0.1 # the wavelength used to compare if w < 5: # wavelength is short enough to not have numerical # issues. Evaluate both at w. @@ -227,12 +233,12 @@ def dust_continuum_guess(row): temp = row["temperature"][0] fmax_lam = 2898.0 / temp bb = BlackBody1D(1, temp) - if fmax_lam >= wmin and fmax_lam <= wmax: + if fmax_lam >= lam_min and fmax_lam <= lam_max: lam_ref = fmax_lam - elif fmax_lam > wmax: - lam_ref = wmax + elif fmax_lam > lam_max: + lam_ref = lam_max else: - lam_ref = wmin + lam_ref = lam_min flux_ref = np.median(flux[(lam > lam_ref - 0.2) & (lam < lam_ref + 0.2)]) amp_guess = flux_ref / bb(lam_ref) @@ -241,47 +247,59 @@ def dust_continuum_guess(row): loop_over_non_fixed("dust_continuum", "tau", dust_continuum_guess) def line_fwhm_guess(row): - w = row["wavelength"][0] - if not instrument.within_segment(w, inst): + lam_line = row["wavelength"][0] + if not instrument.within_segment(lam_line, inst): return 0 - fwhm = instrument.fwhm(inst, w, as_bounded=True)[0][0] + fwhm = instrument.fwhm(inst, lam_line, as_bounded=True)[0][0] return fwhm - def amp_guess(row, fwhm): - w = row["wavelength"][0] - if not instrument.within_segment(w, inst): + def power_guess(row, fwhm): + # local integration for the lines + lam_line = row["wavelength"][0] + if not instrument.within_segment(lam_line, inst): return 0 factor = 1.5 - wmin = w - factor * fwhm - wmax = w + factor * fwhm - lam_window = (lam > wmin) & (lam < wmax) + lam_min = lam_line - factor * fwhm + lam_max = lam_line + factor * fwhm + lam_window = (lam > lam_min) & (lam < lam_max) xpoints = lam[lam_window] ypoints = flux[lam_window] if np.count_nonzero(lam_window) >= 2: # difference between flux in window and flux around it - power_guess = integrate.trapezoid(flux[lam_window], lam[lam_window]) + Fnu_dlambda = integrate.trapezoid(flux[lam_window], lam[lam_window]) # subtract continuum estimate, but make sure we don't go negative continuum = (ypoints[0] + ypoints[-1]) / 2 * (xpoints[-1] - xpoints[0]) - if continuum < power_guess: - power_guess -= continuum + if continuum < Fnu_dlambda: + Fnu_dlambda -= continuum else: - power_guess = 0 - return power_guess / fwhm + Fnu_dlambda = 0 + + # this is an unphysical power (Fnu * dlambda), but we + # convert to Fnu dnu = Fnu dnu/dlambda dlambda = Fnu c / + # lambda **2 dlambda + Fnu_dlambda *= units.intensity * units.wavelength + Fnu_dnu = Fnu_dlambda * constants.c / (lam_line * units.wavelength) ** 2 + return Fnu_dnu.to(units.intensity_power).value - # Same logic as in the old function: just use same amp for all - # dust features. - some_flux = 0.5 * np.median(flux) - loop_over_non_fixed("dust_feature", "power", lambda row: some_flux) + def drude_power_guess(row): + # multiply total power by some fraction to guess Drude power + fwhm = row["fwhm"][0] * units.wavelength + delta_w = spec.spectral_axis[-1] - spec.spectral_axis[0] + return (total_power * fwhm / delta_w).to(units.intensity_power).value + + loop_over_non_fixed("dust_feature", "power", drude_power_guess) if integrate_line_flux: - # calc line amplitude using instrumental fwhm and integral over data + # calc line power using instrumental fwhm and integral over data loop_over_non_fixed( - "line", "power", lambda row: amp_guess(row, line_fwhm_guess(row)) + "line", "power", lambda row: power_guess(row, line_fwhm_guess(row)) ) else: - loop_over_non_fixed("line", "power", lambda row: some_flux) + loop_over_non_fixed( + "line", "power", lambda row: median_flux * line_fwhm_guess(row) + ) # Set the fwhms in the features table. Slightly different logic, # as the fwhm for lines are masked by default. TODO: leave FWHM @@ -745,14 +763,16 @@ def tabulate( if wavelengths is None: ranges = instrument.wave_range(instrumentname) if isinstance(ranges[0], float): - wmin, wmax = ranges + lam_min, lam_max = ranges else: # In case of multiple ranges (multiple segments), choose # the min and max instead - wmin = min(r[0] for r in ranges) - wmax = max(r[1] for r in ranges) - wfwhm = instrument.fwhm(instrumentname, wmin, as_bounded=True)[0, 0] - lam = np.arange(wmin, wmax, wfwhm / 2) * u.micron + lam_min = min(r[0] for r in ranges) + lam_max = max(r[1] for r in ranges) + + wfwhm = instrument.fwhm(instrumentname, lam_min, as_bounded=True)[0, 0] + lam = np.arange(lam_min, lam_max, wfwhm / 2) * u.micron + elif isinstance(wavelengths, Spectrum1D): lam = wavelengths.spectral_axis.to(u.micron) / (1 + z) else: @@ -902,13 +922,12 @@ def cleaned(features_tuple3): elif kind == "line": # be careful with lines that have masked FWHM values here - if use_instrument_fwhm or row['fwhm'] is np.ma.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 # 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 From ddec47b92dc9a0fab069b533371a08a434334635 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Thu, 27 Jun 2024 15:47:17 -0400 Subject: [PATCH 61/68] Remove equation with nu0 from PowerDrude docstring --- pahfit/fitters/ap_components.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pahfit/fitters/ap_components.py b/pahfit/fitters/ap_components.py index 7604a61a..68c05494 100644 --- a/pahfit/fitters/ap_components.py +++ b/pahfit/fitters/ap_components.py @@ -197,14 +197,14 @@ def evaluate(x, power, x_0, fwhm): Which can be solved for the amplitude b. - b = (P * 2 * x0) / (pi * c * g) = 2P / (pi nu0 g). + b = (P * 2 * x0) / (pi * c * g) According to the above equations, without additional conversions, the resulting amplitude unit will be unit(P) * - Hz-1. This will result in very small values for for Inu(lambda), - or very large values for P. To avoid numerical problems with the - fitting algorithm, we apply conversions so that Inu(lambda) and - P are in the internal units. + Hz-1. The factor x0 / c = nu0 (Hz-1) will result in very small + values for for Inu(lambda), or very large values for P. To avoid + numerical problems, we apply a conversion that ensures + Inu(lambda) and P are in internal units. Parameters ---------- From aa155145a6b4365bca179b8eaf9700a552907e88 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Fri, 28 Jun 2024 18:31:29 -0400 Subject: [PATCH 62/68] Set power unit in Features table to intensity_power --- 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 0efed64b..b563bbef 100644 --- a/pahfit/features/features.py +++ b/pahfit/features/features.py @@ -41,7 +41,7 @@ PARAM_UNITS = {'temperature': pahfit.units.temperature, 'wavelength': pahfit.units.wavelength, 'fwhm': pahfit.units.wavelength, - 'power': pahfit.units.intensity} + 'power': pahfit.units.intensity_power} class UniqueKeyLoader(yaml.SafeLoader): From 79cac6f865a0a95f46707258f9e231b9996fef5d Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 9 Jul 2024 16:55:42 -0400 Subject: [PATCH 63/68] Change trapz to trapezoid to fix deprecation warning --- pahfit/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pahfit/model.py b/pahfit/model.py index 347b5538..aa75bc4c 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -194,7 +194,7 @@ def guess( # Some useful quantities for guessing median_flux = np.median(flux) Flambda = flux * units.intensity * (lam * units.wavelength) ** -2 * constants.c - total_power = np.trapz(Flambda, lam * units.wavelength) + total_power = np.trapezoid(Flambda, lam * units.wavelength) # simple linear interpolation function for spectrum sp = interpolate.interp1d(lam, flux) From 8310da10475a5b434e562123ec5303e6f0b6e38a Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 9 Jul 2024 16:59:06 -0400 Subject: [PATCH 64/68] Formatting --- pahfit/fitters/ap_components.py | 9 +++++---- pahfit/fitters/ap_fitter.py | 2 +- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/pahfit/fitters/ap_components.py b/pahfit/fitters/ap_components.py index 68c05494..b05352c5 100644 --- a/pahfit/fitters/ap_components.py +++ b/pahfit/fitters/ap_components.py @@ -22,12 +22,11 @@ class BlackBody1D(Fittable1DModel): @staticmethod def evaluate(x, amplitude, temperature): - """ - """ + """ """ return ( amplitude * 3.97289e13 - / x ** 3 + / x**3 / (np.exp(1.4387752e4 / x / temperature) - 1.0) ) @@ -84,7 +83,9 @@ def kvt(in_x): # Extend kvt profile to shorter wavelengths if min(in_x) < min(kvt_wav): kvt_wav_short = in_x[in_x < min(kvt_wav)] - kvt_int_short_tmp = min(kvt_int) * np.exp(2.03 * (kvt_wav_short - min(kvt_wav))) + kvt_int_short_tmp = min(kvt_int) * np.exp( + 2.03 * (kvt_wav_short - min(kvt_wav)) + ) # Since kvt_int_shoft_tmp does not reach min(kvt_int), # we scale it to stitch it. kvt_int_short = kvt_int_short_tmp * (kvt_int[0] / max(kvt_int_short_tmp)) diff --git a/pahfit/fitters/ap_fitter.py b/pahfit/fitters/ap_fitter.py index af189341..af8be31a 100644 --- a/pahfit/fitters/ap_fitter.py +++ b/pahfit/fitters/ap_fitter.py @@ -6,7 +6,7 @@ S07_attenuation, att_Drude1D, PowerDrude1D, - PowerGaussian1D + PowerGaussian1D, ) from astropy.modeling.fitting import LevMarLSQFitter import numpy as np From c500ae8a8c21d1fce9a49d5f7c39edd0e1585c62 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 13 Aug 2024 15:58:18 -0400 Subject: [PATCH 65/68] Fix typo in PowerDrude1D docstring --- pahfit/fitters/ap_components.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pahfit/fitters/ap_components.py b/pahfit/fitters/ap_components.py index b05352c5..2940a375 100644 --- a/pahfit/fitters/ap_components.py +++ b/pahfit/fitters/ap_components.py @@ -202,7 +202,7 @@ def evaluate(x, power, x_0, fwhm): According to the above equations, without additional conversions, the resulting amplitude unit will be unit(P) * - Hz-1. The factor x0 / c = nu0 (Hz-1) will result in very small + Hz-1. The factor x0 / c = 1 / nu0 (Hz-1) will result in very small values for for Inu(lambda), or very large values for P. To avoid numerical problems, we apply a conversion that ensures Inu(lambda) and P are in internal units. From 4c222073ce9cb8aec15004bb889d062b916f6a63 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 13 Aug 2024 15:59:38 -0400 Subject: [PATCH 66/68] Use scipy.integrate.trapezoid instead of np.trapezoid --- pahfit/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pahfit/model.py b/pahfit/model.py index aa75bc4c..70b55ae1 100644 --- a/pahfit/model.py +++ b/pahfit/model.py @@ -194,7 +194,7 @@ def guess( # Some useful quantities for guessing median_flux = np.median(flux) Flambda = flux * units.intensity * (lam * units.wavelength) ** -2 * constants.c - total_power = np.trapezoid(Flambda, lam * units.wavelength) + total_power = integrate.trapezoid(Flambda, lam * units.wavelength) # simple linear interpolation function for spectrum sp = interpolate.interp1d(lam, flux) From 936a2fa5748dfe13596976de0c4aab040299e7da Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 13 Aug 2024 16:02:22 -0400 Subject: [PATCH 67/68] Access intensity_amplitude_factor via self in power features --- pahfit/fitters/ap_components.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/pahfit/fitters/ap_components.py b/pahfit/fitters/ap_components.py index 2940a375..ed306c73 100644 --- a/pahfit/fitters/ap_components.py +++ b/pahfit/fitters/ap_components.py @@ -177,8 +177,7 @@ class PowerDrude1D(Fittable1DModel): .value ) - @staticmethod - def evaluate(x, power, x_0, fwhm): + def evaluate(self, x, power, x_0, fwhm): """Smith, et al. (2007) dust features model. Calculation is for a Drude profile (equation in section 4.1.4). @@ -221,7 +220,7 @@ def evaluate(x, power, x_0, fwhm): # factor = (2 * unit(power) * unit(wavelength) / (pi * c)).to(unit(intensity)) g = fwhm / x_0 - b = power * x_0 / g * PowerDrude1D.intensity_amplitude_factor + b = power * x_0 / g * self.intensity_amplitude_factor return b * g**2 / ((x / x_0 - x_0 / x) ** 2 + g**2) @@ -271,12 +270,11 @@ class PowerGaussian1D(Fittable1DModel): .value ) - @staticmethod - def evaluate(x, power, mean, stddev): + def evaluate(self, x, power, mean, stddev): """Evaluate F_nu(lambda) given the power. See class description for equations and unit notes.""" # amplitude in intensity units - Anu = power * mean**2 / stddev * PowerGaussian1D.intensity_amplitude_factor + Anu = power * mean**2 / stddev * self.intensity_amplitude_factor return Anu * np.exp(-0.5 * np.square((x - mean) / stddev)) From a6434175987f3e219229058379d2bcbe7b5163b3 Mon Sep 17 00:00:00 2001 From: Dries Van De Putte Date: Tue, 13 Aug 2024 16:52:41 -0400 Subject: [PATCH 68/68] Clean up some more docstrings in ap_components --- pahfit/fitters/ap_components.py | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/pahfit/fitters/ap_components.py b/pahfit/fitters/ap_components.py index ed306c73..921e23b3 100644 --- a/pahfit/fitters/ap_components.py +++ b/pahfit/fitters/ap_components.py @@ -140,7 +140,8 @@ def evaluate(x, tau, x_0, fwhm): class PowerDrude1D(Fittable1DModel): - """Drude profile with amplitude determined by power. + """ + Drude profile with amplitude determined by power. Special attention is needed for the units. An implementation for this is 'unitful' because 'power' is defined as integral over @@ -178,11 +179,11 @@ class PowerDrude1D(Fittable1DModel): ) def evaluate(self, x, power, x_0, fwhm): - """Smith, et al. (2007) dust features model. Calculation is for - a Drude profile (equation in section 4.1.4). + """ + Smith, et al. (2007) dust features model. Calculation is for a + Drude profile (equation in section 4.1.4). - The intensity profile as a function of - wavelength is + The intensity profile as a function of wavelength is Inu(lambda) = (b * g**2) / ((lambda / x0 - x0 / lambda)**2 + g**2) @@ -201,9 +202,9 @@ def evaluate(self, x, power, x_0, fwhm): According to the above equations, without additional conversions, the resulting amplitude unit will be unit(P) * - Hz-1. The factor x0 / c = 1 / nu0 (Hz-1) will result in very small - values for for Inu(lambda), or very large values for P. To avoid - numerical problems, we apply a conversion that ensures + Hz-1. The factor x0 / c = 1 / nu0 (Hz-1) will result in very + small values for for Inu(lambda), or very large values for P. To + avoid numerical problems, we apply a conversion that ensures Inu(lambda) and P are in internal units. Parameters @@ -225,7 +226,8 @@ def evaluate(self, x, power, x_0, fwhm): class PowerGaussian1D(Fittable1DModel): - """Gaussian profile with amplitude derived from power. + """ + Gaussian profile with amplitude derived from power. Implementation analogous to PowerDrude1D. @@ -236,7 +238,7 @@ class PowerGaussian1D(Fittable1DModel): Flambda(lambda) = Alambda * G(lambda; mean, stddev) - where G is a gaussian profile with amplitude 1.Converting this to + where G is a gaussian profile with amplitude 1. Converting this to Fnu units yields Fnu(lambda) = lambda**2 / c * Flambda(lambda) @@ -271,7 +273,8 @@ class PowerGaussian1D(Fittable1DModel): ) def evaluate(self, x, power, mean, stddev): - """Evaluate F_nu(lambda) given the power. + """ + Evaluate F_nu(lambda) given the power. See class description for equations and unit notes."""