Standardize units

pahfit/features/features.py

@@ -21,16 +21,29 @@
 import numpy as np
 from astropy.table import vstack, Table, TableAttribute
 from astropy.io.misc.yaml import yaml
-import astropy.units as u
 from pkg_resources import resource_filename
 from pahfit.errors import PAHFITFeatureError
 from pahfit.features.features_format import BoundedMaskedColumn, BoundedParTableFormatter
+import pahfit.units
+
+# Feature kinds and associated parameters
+KIND_PARAMS = {'starlight': {'temperature', 'tau'},
+               'dust_continuum': {'temperature', 'tau'},
+               'line': {'wavelength', 'power'},  # 'fwhm', Instrument Pack detail!
+               'dust_feature': {'wavelength', 'fwhm', 'power'},
+               'attenuation': {'model', 'tau', 'geometry'},
+               'absorption': {'wavelength', 'fwhm', 'tau', 'geometry'}}
+
+# Parameter default units: flux density/intensity/power (other units determined on fit)
+PARAM_UNITS = {'temperature': pahfit.units.temperature,
+               'wavelength': pahfit.units.wavelength,
+               'fwhm': pahfit.units.wavelength}
 
 
 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.")
@@ -41,44 +54,46 @@ def construct_mapping(self, node, deep=False):
 def value_bounds(val, bounds):
     """Compute bounds for a bounded value.
 
-    Arguments:
+    Parameters
     ----------
-      val: The value to bound.
+      val : float
+          The value to bound.
 
-      bounds: Either None for no relevant bounds (i.e. fixed), or a
-        two element iterable specifying (min, max) bounds.  Each of
-        min and max can be a numerical value, None (for infinite
-        bounds, either negative or positive, as appropriate), or a
-        string ending in:
+      bounds : float or str iterable, or None
+          Either None for no relevant bounds (i.e. fixed), or a two
+          element iterable specifying (min, max) bounds.  Each of min
+          and max can be a numerical value, None (for infinite bounds,
+          either negative or positive, as appropriate), or a string
+          ending in:
 
-          #: an absolute offset from the value
-          %: a percentage offset from the value
+            #: an absolute offset from the value
+            %: a percentage offset from the value
 
-        Offsets are necessarily negative for min bound, positive for
-        max bounds.
+          Offsets are necessarily negative for min bound, positive for
+          max bounds.
 
-    Examples:
-    ---------
+    Returns:
+    -------
+
+      The value and bounds as a 3 element tuple (value, min, max).
+          Any missing bound is replaced with the numpy `masked' value.
+
+    Notes
+    --------
 
       A bound of ('-1.5%', '0%') would indicate a minimum bound
         1.5% below the value, and a max bound at the value itself.
 
       A bound of ('-0.1#', None) would indicate a minimum bound 0.1 below
         the value, and infinite maximum bound.
 
-    Returns:
-    -------
-
-      The value, if unbounded, or a 3 element tuple (value, min, max).
-        Any missing bound is replaced with the numpy `masked' value.
-
     Raises:
     -------
 
       ValueError: if bounds are specified and the value does not fall
-        between them.
-
+          between them.
     """
+
     if val is None:
         val = np.ma.masked
     if not bounds:
@@ -102,34 +117,46 @@ 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}
+    _param_attrs = set(('value', 'bounds'))  # params can have these attributes
     _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)
 
     @classmethod
     def read(cls, file, *args, **kwargs):
         """Read a table from file.
 
-        If reading a YAML file, read it in as a science pack and
-        return the new table. Otherwise, use astropy's normal Table
+        Parameters
+        ----------
+
+          file : str
+              The name of the file to read, either a full valid path,
+              or named file in the PAHFIT science_packs directory.
+
+        Returns
+        -------
+          table : Features
+              A filled `~pahfit.features.Features` table.
+
+        Notes
+        -----
+        If reading a YAML file, reads it in as a science pack and
+        return the new table. Otherwise, uses astropy's normal Table
         reader.
         """
         if file.endswith(".yaml") or file.endswith(".yml"):
@@ -143,16 +170,7 @@ def read(cls, file, *args, **kwargs):
     def _read_scipack(cls, file):
         """Read a science pack specification from YAML file.
 
-        Arguments:
-        ----------
-
-          file: the name of the file, either a full valid path, or
-            named file in the PAHFIT science_packs directory.!
-
-        Returns:
-        --------
-
-          table: A filled pahfit.features.Features table.
+        For parameters and return, see `read`.
         """
 
         feat_tables = dict()
@@ -179,7 +197,7 @@ def _read_scipack(cls, file):
                 raise PAHFITFeatureError(f"No kind found for {name}\n\t{file}")
 
             try:
-                valid_params = cls._kind_params[kind]
+                valid_params = KIND_PARAMS[kind]
             except KeyError:
                 raise PAHFITFeatureError(f"Unknown kind {kind} for {name}\n\t{file}")
             unknown_params = [x for x in keys
@@ -264,7 +282,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]
@@ -299,36 +317,46 @@ def _add_feature(cls, kind: str, t: dict, name: str, *,
 
     @classmethod
     def _construct_table(cls, inp: dict):
-        """Construct a masked table with units from input dictionary
-        INP.  INP is a dictionary with feature names as the key, and a
+        """Construct a masked table from input dictionary INP.
+        INP is a dictionary with feature names as the keys, and a
         dictionary of feature parameters as value.  Each value in the
         feature parameter dictionary is either a value or tuple of 3
-        values for bounds.
+        values, expressing bounds.
         """
         tables = []
         for (kind, features) in inp.items():
-            kind_params = cls._kind_params[kind]  # All params for this kind
+            if kind == "_ratios":
+                continue
+            kp = KIND_PARAMS[kind]  # All params for this kind
             rows = []
             for (name, params) in features.items():
-                for missing in kind_params - params.keys():
+                for missing in kp - params.keys():
                     if missing in cls._no_bounds:
                         params[missing] = 0.0
                     else:
                         params[missing] = value_bounds(0.0, bounds=(0.0, None))
                 rows.append(dict(name=name, **params))
             table_columns = rows[0].keys()
             t = cls(rows, names=table_columns)
-            for p in cls._kind_params[kind]:
+            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
 
+    @staticmethod
+    def _index_table(tbl):
+        for indx in ('name', 'group'):
+            tbl.add_index(indx)
+
     def mask_feature(self, name, mask_value=True):
         """Mask all the parameters of a feature.
 
@@ -345,7 +373,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

pahfit/units.py

@@ -0,0 +1,17 @@
+import astropy.units as u
+from astropy.units import CompositeUnit
+
+
+# Working/default unitParameter default units: flux density/intensity/power
+# These are PAHFITs default science packs parameter and output units
+temperature = u.K
+wavelength = u.um
+flux_density = u.mJy
+flux_power = CompositeUnit(1e-22, (u.W, u.m), (1, -2))
+solid_angle = u.sr
+intensity = u.MJy / u.sr
+intensity_power = CompositeUnit(1e-10, (u.W, u.m, u.sr), (1, -2, -1))
+
+# Note: integrated power units of 1e-22 W/m^2 (from flux) corresponds
+# to the unit 1e-10 W/m^2/sr (from intensity) if it occurs uniformly
+# over a solid angle 0.21" on a side (about a small JWST IFU pixel)