Repository: astropy/astropy

Issue:
LombScargle.false_alarm_probability fails for Quantities
Trying to calculate the Lomb Scragle false alarm probability while using Quantities with units throws an error while executing `fap_baluev`, since the exponent argument has units.
(It works ok when using `method='naive'`)

Using Ubuntu 18.04.3, python 3.7.4, astropy 3.1.2, numpy 1.16.2.

Minimal working example:
```python
import numpy as np
import astropy.units as u
from astropy.timeseries import LombScargle

t = 100 * np.random.rand(10) * u.day
dy = 0.1 * u.mag
y = (np.sin(2 * np.pi * t.value) + dy.value * np.random.randn(10)) * u.mag

ls = LombScargle(t, y, dy)
frequency, power = ls.autopower()
ls.false_alarm_probability(power.max())
```

EDIT: Syntax highlighting

Discussion:

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index b606e3c044..c065ceb98f 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -210,10 +210,6 @@ astropy.io.fits
 
 - Masked column handling has changed, see ``astropy.table`` entry below. [#8789]
 
-- ``io.fits.Header`` has been made safe for subclasses for copying and slicing.
-  As a result of this change, the private subclass ``CompImageHeader`` now always
-  should be passed an explicit ``image_header``. [#9229]
-
 astropy.io.registry
 ^^^^^^^^^^^^^^^^^^^
 
@@ -533,6 +529,9 @@ astropy.time
 astropy.timeseries
 ^^^^^^^^^^^^^^^^^^
 
+- Fixed handling of ``Quantity`` input data for all methods of
+  ``LombScarge.false_alarm_probabilty``. [#9246]
+
 astropy.uncertainty
 ^^^^^^^^^^^^^^^^^^^
 
diff --git a/astropy/io/fits/hdu/compressed.py b/astropy/io/fits/hdu/compressed.py
index d1a922e72e..f288a8ae52 100644
--- a/astropy/io/fits/hdu/compressed.py
+++ b/astropy/io/fits/hdu/compressed.py
@@ -85,9 +85,6 @@ class CompImageHeader(Header):
     This essentially wraps the image header, so that all values are read from
     and written to the image header.  However, updates to the image header will
     also update the table header where appropriate.
-
-    Note that if no image header is passed in, the code will instantiate a
-    regular `~astropy.io.fits.Header`.
     """
 
     # TODO: The difficulty of implementing this screams a need to rewrite this
@@ -108,19 +105,9 @@ class CompImageHeader(Header):
     # the schema system, but it's not quite ready for that yet.  Also it still
     # makes more sense to change CompImageHDU to subclass ImageHDU :/
 
-    def __new__(cls, table_header, image_header=None):
-        # 2019-09-14 (MHvK): No point wrapping anything if no image_header is
-        # given.  This happens if __getitem__ and copy are called - our super
-        # class will aim to initialize a new, possibly partially filled
-        # header, but we cannot usefully deal with that.
-        # TODO: the above suggests strongly we should *not* subclass from
-        # Header.  See also comment above about the need for reorganization.
+    def __init__(self, table_header, image_header=None):
         if image_header is None:
-            return Header(table_header)
-        else:
-            return super().__new__(cls)
-
-    def __init__(self, table_header, image_header):
+            image_header = Header()
         self._cards = image_header._cards
         self._keyword_indices = image_header._keyword_indices
         self._rvkc_indices = image_header._rvkc_indices
diff --git a/astropy/io/fits/header.py b/astropy/io/fits/header.py
index f60ef93049..abd253e4ad 100644
--- a/astropy/io/fits/header.py
+++ b/astropy/io/fits/header.py
@@ -132,10 +132,10 @@ class Header:
 
     def __getitem__(self, key):
         if isinstance(key, slice):
-            return self.__class__([copy.copy(c) for c in self._cards[key]])
+            return Header([copy.copy(c) for c in self._cards[key]])
         elif self._haswildcard(key):
-            return self.__class__([copy.copy(self._cards[idx])
-                                   for idx in self._wildcardmatch(key)])
+            return Header([copy.copy(self._cards[idx])
+                           for idx in self._wildcardmatch(key)])
         elif (isinstance(key, str) and
               key.upper() in Card._commentary_keywords):
             key = key.upper()
@@ -824,7 +824,7 @@ class Header:
             A new :class:`Header` instance.
         """
 
-        tmp = self.__class__((copy.copy(card) for card in self._cards))
+        tmp = Header(copy.copy(card) for card in self._cards)
         if strip:
             tmp._strip()
         return tmp
@@ -1294,7 +1294,7 @@ class Header:
             new cards to the header.
         """
 
-        temp = self.__class__(cards)
+        temp = Header(cards)
         if strip:
             temp._strip()
 
diff --git a/astropy/timeseries/periodograms/lombscargle/_statistics.py b/astropy/timeseries/periodograms/lombscargle/_statistics.py
index 78c87ec9ba..b5aa8fbb43 100644
--- a/astropy/timeseries/periodograms/lombscargle/_statistics.py
+++ b/astropy/timeseries/periodograms/lombscargle/_statistics.py
@@ -10,6 +10,8 @@ from functools import wraps
 
 import numpy as np
 
+from astropy import units as u
+
 
 def _weighted_sum(val, dy):
     if dy is not None:
@@ -22,7 +24,7 @@ def _weighted_mean(val, dy):
     if dy is None:
         return val.mean()
     else:
-        return _weighted_sum(val, dy) / _weighted_sum(np.ones_like(val), dy)
+        return _weighted_sum(val, dy) / _weighted_sum(np.ones(val.shape), dy)
 
 
 def _weighted_var(val, dy):
@@ -342,25 +344,30 @@ def inv_fap_baluev(p, fmax, t, y, dy, normalization='standard'):
     return res.x
 
 
-def _bootstrap_max(t, y, dy, fmax, normalization, random_seed):
+def _bootstrap_max(t, y, dy, fmax, normalization, random_seed, n_bootstrap=1000):
     """Generate a sequence of bootstrap estimates of the max"""
     from .core import LombScargle
     rng = np.random.RandomState(random_seed)
-    while True:
+    power_max = []
+    for _ in range(n_bootstrap):
         s = rng.randint(0, len(y), len(y))  # sample with replacement
         ls_boot = LombScargle(t, y[s], dy if dy is None else dy[s],
                               normalization=normalization)
         freq, power = ls_boot.autopower(maximum_frequency=fmax)
-        yield power.max()
+        power_max.append(power.max())
+
+    power_max = u.Quantity(power_max)
+    power_max.sort()
+
+    return power_max
 
 
 def fap_bootstrap(Z, fmax, t, y, dy, normalization='standard',
                   n_bootstraps=1000, random_seed=None):
     """Bootstrap estimate of the false alarm probability"""
-    pmax = np.fromiter(_bootstrap_max(t, y, dy, fmax,
-                                      normalization, random_seed),
-                       float, n_bootstraps)
-    pmax.sort()
+    pmax = _bootstrap_max(t, y, dy, fmax, normalization, random_seed,
+                          n_bootstraps)
+
     return 1 - np.searchsorted(pmax, Z) / len(pmax)
 
 
@@ -368,10 +375,9 @@ def inv_fap_bootstrap(fap, fmax, t, y, dy, normalization='standard',
                       n_bootstraps=1000, random_seed=None):
     """Bootstrap estimate of the inverse false alarm probability"""
     fap = np.asarray(fap)
-    pmax = np.fromiter(_bootstrap_max(t, y, dy, fmax,
-                                      normalization, random_seed),
-                       float, n_bootstraps)
-    pmax.sort()
+    pmax = _bootstrap_max(t, y, dy, fmax, normalization, random_seed,
+                          n_bootstraps)
+
     return pmax[np.clip(np.floor((1 - fap) * len(pmax)).astype(int),
                         0, len(pmax) - 1)]
 


Test patch:
diff --git a/astropy/io/fits/tests/test_header.py b/astropy/io/fits/tests/test_header.py
index e1c4df57cc..54013cd199 100644
--- a/astropy/io/fits/tests/test_header.py
+++ b/astropy/io/fits/tests/test_header.py
@@ -2838,36 +2838,3 @@ class TestRecordValuedKeywordCards(FitsTestCase):
         with open(self.temp('mode.fits'), mode=mode) as ff:
             hdu = fits.ImageHDU(data=np.ones(5))
             hdu.writeto(ff)
-
-
-def test_subclass():
-    """Check that subclasses don't get ignored on slicing and copying."""
-    class MyHeader(fits.Header):
-        def append(self, card, *args, **kwargs):
-            if isinstance(card, tuple) and len(card) == 2:
-                # Just for our checks we add a comment if there is none.
-                card += ('no comment',)
-
-            return super().append(card, *args, **kwargs)
-
-    my_header = MyHeader((('a', 1., 'first'),
-                          ('b', 2., 'second'),
-                          ('c', 3.,)))
-
-    assert my_header.comments['a'] == 'first'
-    assert my_header.comments['b'] == 'second'
-    assert my_header.comments['c'] == 'no comment'
-
-    slice_ = my_header[1:]
-    assert type(slice_) is MyHeader
-    assert slice_.comments['b'] == 'second'
-    assert slice_.comments['c'] == 'no comment'
-    selection = my_header['c*']
-    assert type(selection) is MyHeader
-    assert selection.comments['c'] == 'no comment'
-    copy_ = my_header.copy()
-    assert type(copy_) is MyHeader
-    assert copy_.comments['b'] == 'second'
-    assert copy_.comments['c'] == 'no comment'
-    my_header.extend((('d', 4.),))
-    assert my_header.comments['d'] == 'no comment'
diff --git a/astropy/timeseries/periodograms/lombscargle/tests/test_statistics.py b/astropy/timeseries/periodograms/lombscargle/tests/test_statistics.py
index 8adec8cd6c..90a3849a63 100644
--- a/astropy/timeseries/periodograms/lombscargle/tests/test_statistics.py
+++ b/astropy/timeseries/periodograms/lombscargle/tests/test_statistics.py
@@ -9,6 +9,7 @@ except ImportError:
 else:
     HAS_SCIPY = True
 
+import astropy.units as u
 from astropy.timeseries.periodograms.lombscargle import LombScargle
 from astropy.timeseries.periodograms.lombscargle._statistics import (cdf_single, pdf_single, fap_single, inv_fap_single,
                            METHODS)
@@ -18,7 +19,7 @@ METHOD_KWDS = dict(bootstrap={'n_bootstraps': 20, 'random_seed': 42})
 NORMALIZATIONS = ['standard', 'psd', 'log', 'model']
 
 
-def make_data(N=100, period=1, theta=[10, 2, 3], dy=1, rseed=0):
+def make_data(N=100, period=1, theta=[10, 2, 3], dy=1, rseed=0, units=False):
     """Generate some data for testing"""
     rng = np.random.RandomState(rseed)
     t = 5 * period * rng.rand(N)
@@ -27,27 +28,37 @@ def make_data(N=100, period=1, theta=[10, 2, 3], dy=1, rseed=0):
     dy = dy * (0.5 + rng.rand(N))
     y += dy * rng.randn(N)
 
-    return t, y, dy
+    fmax = 5
+
+    if units:
+        return t * u.day, y * u.mag, dy * u.mag, fmax / u.day
+    else:
+        return t, y, dy, fmax
 
 
-@pytest.fixture
-def null_data(N=1000, dy=1, rseed=0):
+def null_data(N=1000, dy=1, rseed=0, units=False):
     """Generate null hypothesis data"""
     rng = np.random.RandomState(rseed)
     t = 100 * rng.rand(N)
     dy = 0.5 * dy * (1 + rng.rand(N))
     y = dy * rng.randn(N)
-    return t, y, dy
+    fmax = 40
+
+    if units:
+        return t * u.day, y * u.mag, dy * u.mag, fmax / u.day
+    else:
+        return t, y, dy, fmax
 
 
 @pytest.mark.parametrize('normalization', NORMALIZATIONS)
 @pytest.mark.parametrize('with_errors', [True, False])
-def test_distribution(null_data, normalization, with_errors, fmax=40):
-    t, y, dy = null_data
+@pytest.mark.parametrize('units', [False, True])
+def test_distribution(normalization, with_errors, units):
+    t, y, dy, fmax = null_data(units=units)
+
     if not with_errors:
         dy = None
 
-    N = len(t)
     ls = LombScargle(t, y, dy, normalization=normalization)
     freq, power = ls.autopower(maximum_frequency=fmax)
     z = np.linspace(0, power.max(), 1000)
@@ -57,6 +68,8 @@ def test_distribution(null_data, normalization, with_errors, fmax=40):
     z_mid = z[:-1] + 0.5 * dz
     pdf = ls.distribution(z_mid)
     cdf = ls.distribution(z, cumulative=True)
+    if isinstance(dz, u.Quantity):
+        dz = dz.value
     assert_allclose(pdf, np.diff(cdf) / dz, rtol=1E-5, atol=1E-8)
 
     # psd normalization without specified errors produces bad results
@@ -80,8 +93,9 @@ def test_inverse_single(N, normalization):
 
 @pytest.mark.parametrize('normalization', NORMALIZATIONS)
 @pytest.mark.parametrize('use_errs', [True, False])
-def test_inverse_bootstrap(null_data, normalization, use_errs, fmax=5):
-    t, y, dy = null_data
+@pytest.mark.parametrize('units', [False, True])
+def test_inverse_bootstrap(normalization, use_errs, units):
+    t, y, dy, fmax = null_data(units=units)
     if not use_errs:
         dy = None
 
@@ -105,11 +119,12 @@ def test_inverse_bootstrap(null_data, normalization, use_errs, fmax=5):
 @pytest.mark.parametrize('normalization', NORMALIZATIONS)
 @pytest.mark.parametrize('use_errs', [True, False])
 @pytest.mark.parametrize('N', [10, 100, 1000])
-def test_inverses(method, normalization, use_errs, N, T=5, fmax=5):
+@pytest.mark.parametrize('units', [False, True])
+def test_inverses(method, normalization, use_errs, N, units, T=5):
     if not HAS_SCIPY and method in ['baluev', 'davies']:
         pytest.skip("SciPy required")
 
-    t, y, dy = make_data(N, rseed=543)
+    t, y, dy, fmax = make_data(N, rseed=543, units=units)
     if not use_errs:
         dy = None
     method_kwds = METHOD_KWDS.get(method, None)
@@ -128,13 +143,13 @@ def test_inverses(method, normalization, use_errs, N, T=5, fmax=5):
 
 @pytest.mark.parametrize('method', sorted(METHODS))
 @pytest.mark.parametrize('normalization', NORMALIZATIONS)
-def test_false_alarm_smoketest(method, normalization):
+@pytest.mark.parametrize('units', [False, True])
+def test_false_alarm_smoketest(method, normalization, units):
     if not HAS_SCIPY and method in ['baluev', 'davies']:
         pytest.skip("SciPy required")
 
     kwds = METHOD_KWDS.get(method, None)
-    t, y, dy = make_data()
-    fmax = 5
+    t, y, dy, fmax = make_data(units=units)
 
     ls = LombScargle(t, y, dy, normalization=normalization)
     freq, power = ls.autopower(maximum_frequency=fmax)
@@ -152,7 +167,8 @@ def test_false_alarm_smoketest(method, normalization):
 @pytest.mark.parametrize('method', sorted(METHODS))
 @pytest.mark.parametrize('use_errs', [True, False])
 @pytest.mark.parametrize('normalization', sorted(set(NORMALIZATIONS) - {'psd'}))
-def test_false_alarm_equivalence(method, normalization, use_errs):
+@pytest.mark.parametrize('units', [False, True])
+def test_false_alarm_equivalence(method, normalization, use_errs, units):
     # Note: the PSD normalization is not equivalent to the others, in that it
     # depends on the absolute errors rather than relative errors. Because the
     # scaling contributes to the distribution, it cannot be converted directly
@@ -161,10 +177,9 @@ def test_false_alarm_equivalence(method, normalization, use_errs):
         pytest.skip("SciPy required")
 
     kwds = METHOD_KWDS.get(method, None)
-    t, y, dy = make_data()
+    t, y, dy, fmax = make_data(units=units)
     if not use_errs:
         dy = None
-    fmax = 5
 
     ls = LombScargle(t, y, dy, normalization=normalization)
     freq, power = ls.autopower(maximum_frequency=fmax)