Merge pull request #257 from pllim/tynt-poc

Initial port of tynt to filter_parameterization subpackage

CHANGES.rst

@@ -1,6 +1,9 @@
 1.2.0 (unreleased)
 ==================
 
+- New ``filter_parameterization`` subpackage to handle filter parameterization,
+  adapted from ``tynt`` package written by Brett Morris. [#257]
+
 - Dropped support for Python 3.6 and 3.7. Minimum supported Python
   version is now 3.8. [#330]
 

docs/index.rst

@@ -287,6 +287,7 @@ Using **synphot**
    synphot/observation
    synphot/formulae
    synphot/units
+   synphot/filter_par
    synphot/tutorials
 
 .. _synphot_history:
@@ -325,6 +326,8 @@ Also imports this C-extension to local namespace:
 
 .. automodapi:: synphot.exceptions
 
+.. automodapi:: synphot.filter_parameterization
+
 .. automodapi:: synphot.models
 
 .. automodapi:: synphot.observation

docs/synphot/filter_par.rst

@@ -0,0 +1,124 @@
+.. _synphot_par_filters:
+
+Parameterized Filters
+=====================
+
+.. note::
+
+    The algorithm for parameterized filters here was originally developed by
+    Brett Morris for the `tynt <https://github.com/bmorris3/tynt/>`_ package.
+
+Filter responses can be approximated using Fast Fourier Transform (FFT).
+If a filter is approximated this way, one only needs to store its FFT
+parameters instead of all the sampled data points. This reduces the
+storage size and increases performance, at the cost of reduced accuracy.
+If you decide to use the parameterization functions provided here,
+it is up to you to decide whether the results are good enough for your
+use cases or not.
+
+.. _filter_fft_generation:
+
+Generating FFT
+--------------
+
+.. testsetup::
+
+    >>> import os
+    >>> from astropy.utils.data import get_pkg_data_filename
+    >>> filename = get_pkg_data_filename(
+    ...     os.path.join('data', 'hst_acs_hrc_f555w.fits'),
+    ...     package='synphot.tests')
+
+You could parameterize a given filter using
+:func:`~synphot.filter_parameterization.filter_to_fft` as follows.
+By default, 10 FFT parameters are returned as complex numbers::
+
+    >>> from synphot import SpectralElement
+    >>> from synphot.filter_parameterization import filter_to_fft
+    >>> filename = 'hst_acs_hrc_f555w.fits'  # doctest: +SKIP
+    >>> bp = SpectralElement.from_file(filename)
+    >>> n_lambda, lambda_0, delta_lambda, tr_max, fft_pars = filter_to_fft(bp)
+    >>> n_lambda  # Number of elements in wavelengths
+    10000
+    >>> lambda_0  # Starting value of wavelengths  # doctest: +FLOAT_CMP
+    <Quantity 3479.999 Angstrom>
+    >>> delta_lambda  # Median wavelength separation  # doctest: +FLOAT_CMP
+    <Quantity 0.66748047 Angstrom>
+    >>> tr_max  # Peak value of throughput  # doctest: +FLOAT_CMP
+    <Quantity 0.241445>
+    >>> fft_pars  # FFT parameters  # doctest: +SKIP
+    [(407.5180314841658+7.494005416219807e-16j),
+     (-78.52240189503877-376.53990235136575j),
+     (-294.86589196496584+127.25464850352665j),
+     (130.20273803287864+190.84263652863257j),
+     (96.62299079012317-91.70087676328245j),
+     (-32.572468348727654-34.227696019221035j),
+     (-8.051741476066471-21.354793540998294j),
+     (-51.708676896903725+6.883836090870033j),
+     (13.08719675518801+54.48177212720124j),
+     (38.635087381362396-13.02803811279449j)]
+
+.. TODO: Only skipping the fft_pars comparison above because output is very
+   different for NUMPY_LT_1_17. Unskip it and replace with +FLOAT_CMP when
+   Numpy minversion is 1.17.
+
+It is up to you to decide how to store this data, though storing it in a
+table format is recommended. In fact, if you have many filters to parameterize,
+:func:`~synphot.filter_parameterization.filters_to_fft_table`
+will store the results in a table for you::
+
+    >>> from synphot.filter_parameterization import filters_to_fft_table
+    >>> mapping = {'HST/ACS/HRC/F555W': (bp, None)}
+    >>> filter_pars_table = filters_to_fft_table(mapping)
+    >>> filter_pars_table  # doctest: +SKIP
+    <Table length=1>
+          filter      n_lambda ...                  fft_9
+                               ...
+          str17        int...  ...                complex128
+    ----------------- -------- ... ---------------------------------------
+    HST/ACS/HRC/F555W    10000 ... (38.635087381362396-13.02803811279449j)
+    >>> filter_pars_table.write('my_filter_pars.fits')  # doctest: +SKIP
+
+.. TODO: Only skipping the filter_pars_table comparison above because output
+   is slightly different for NUMPY_LT_1_17. Unskip it and replace with
+   +FLOAT_CMP +ELLIPSIS when Numpy minversion is 1.17.
+
+.. _filter_fft_construction:
+
+Reconstructing Filter from FFT
+------------------------------
+
+Once you have a parameterized filter (see :ref:`filter_fft_generation`),
+you can reconstruct it for use using
+:func:`~synphot.filter_parameterization.filter_from_fft`.
+Following from the example above::
+
+    >>> from synphot.filter_parameterization import filter_from_fft
+    >>> reconstructed_bp = filter_from_fft(
+    ...     n_lambda, lambda_0, delta_lambda, tr_max, fft_pars)
+
+For this particular example using HST ACS/HRC F555W filter, perhaps 10
+parameters are not quite sufficient. Therefore, caution needs to be exercised
+if you opt to parameterize your filters using this method.
+
+.. plot::
+
+    import os
+    import matplotlib.pyplot as plt
+    from astropy.utils.data import get_pkg_data_filename
+    from synphot import SpectralElement
+    from synphot.filter_parameterization import filter_to_fft, filter_from_fft
+    filename = get_pkg_data_filename(
+        os.path.join('data', 'hst_acs_hrc_f555w.fits'),
+        package='synphot.tests')
+    bp = SpectralElement.from_file(filename)
+    fit_result = filter_to_fft(bp)
+    reconstructed_bp = filter_from_fft(*fit_result)
+    w = bp.waveset
+    plt.plot(w, bp(w), 'b-', label='Original')
+    plt.plot(w, reconstructed_bp(w), 'r--', label='Reconstructed')
+    plt.xlim(3500, 8000)
+    plt.xlabel('Wavelength (Angstrom)')
+    plt.ylabel('Throughput')
+    plt.title('HST ACS/HRC F555W')
+    plt.legend(loc='upper right', numpoints=1)

licenses/TYNT_LICENSE.txt

@@ -0,0 +1,29 @@
+BSD 3-Clause License
+
+Copyright (c) 2019, Brett M. Morris
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+   contributors may be used to endorse or promote products derived from
+   this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

setup.cfg

@@ -61,6 +61,7 @@ docs =
 
 [options.package_data]
 synphot.tests = data/*
+synphot.filter_parameterization.tests = data/*
 
 [flake8]
 # Ignoring these for now:

synphot/compat.py

@@ -2,6 +2,7 @@
 """Module to handle backward-compatibility."""
 
 import astropy
+import numpy
 from astropy.utils.introspection import minversion
 
 try:
@@ -18,11 +19,12 @@
 else:
     HAS_DUST_EXTINCTION = True
 
-
 __all__ = ['ASTROPY_LT_5_0', 'ASTROPY_LT_4_3', 'ASTROPY_LT_4_1',
-           'ASTROPY_LT_4_0', 'HAS_SPECUTILS', 'HAS_DUST_EXTINCTION']
+           'ASTROPY_LT_4_0', 'NUMPY_LT_1_17', 'HAS_SPECUTILS',
+           'HAS_DUST_EXTINCTION']
 
 ASTROPY_LT_5_0 = not minversion(astropy, '4.99')  # astropy<5 but includes 5.0.dev  # noqa
 ASTROPY_LT_4_3 = not minversion(astropy, '4.3')
 ASTROPY_LT_4_1 = not minversion(astropy, '4.1')
 ASTROPY_LT_4_0 = not minversion(astropy, '4.0')
+NUMPY_LT_1_17 = not minversion(numpy, '1.17')

synphot/filter_parameterization/__init__.py

@@ -0,0 +1,9 @@
+"""This subpackage handles filter parameterization.
+
+The algorithms in this subpackage were originally developed by
+Brett Morris as part of the `tynt <https://github.com/bmorris3/tynt>`_
+package.
+
+"""
+
+from .filter_fft import *  # noqa