Repository: StingraySoftware/stingray (python)

Issue:
Shouldn't rebin and rebin_log both return `*Spectrum`s?
As things are now, `Crossspectrum.rebin` returns a `Crossspectrum` object while `Crossspectrum.rebin_log` returns a tuple (frequency, power, number). Is there a specific reason for this? Shouldn't they return the same kind of value (preferably a `Crossspectrum`?)


Pull request description:
rebin_log returns a Powerspectrum/Crossspectrum instance, allows rebin factor
Closes #150, reboots #233

~Two changes here~:

+ `rebin_log` returns a *spectrum instance, consistently with `rebin`

+ `rebin` now accept a rebin _factor_ instead of just the new resolution, both in `Lightcurve` and in {cross,power} spectra

EDIT: Actually things turned out to be far more complicated than this. Tests did not contain the case of calculating lags from rebinned cross spectra. Well, this was were things went terribly wrong: `coherence` wanted some power spectra, that needed to be binned exactly like the main spectrum. So, in the end, I did this:

+ In `{Averaged,}Crossspectrum`, calculate the cross spectrum _and_ the power spectra of the input light curves. This operation (`_make_auxil_spectra`) was put into an `if lc1 is not lc2` condition so that this doesn't become an infinite recursion...

+ When rebinning, rebin the cross spectrum plus each spectrum in the `cs_all` attribute _and_ the input lc's power spectra.

Patch:
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 2a4aee53..0f41724a 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -7,17 +7,14 @@ import scipy.fftpack
 import scipy.optimize
 
 from stingray.lightcurve import Lightcurve
-from stingray.utils import rebin_data, simon
+from stingray.utils import rebin_data, simon, rebin_data_log
 from stingray.exceptions import StingrayError
 from stingray.gti import cross_two_gtis, bin_intervals_from_gtis, check_gtis
+import copy
 
 __all__ = ["Crossspectrum", "AveragedCrossspectrum", "coherence"]
 
 
-def _root_squared_mean(array):
-    return np.sqrt(np.sum(array**2)) / len(array)
-
-
 def coherence(lc1, lc2):
     """
     Estimate coherence function of two light curves.
@@ -135,6 +132,13 @@ class Crossspectrum(object):
         self.lc2 = lc2
 
         self._make_crossspectrum(lc1, lc2)
+        # These are needed to calculate coherence
+        self._make_auxil_pds(lc1, lc2)
+
+    def _make_auxil_pds(self, lc1, lc2):
+        if lc1 is not lc2 and isinstance(lc1, Lightcurve):
+            self.pds1 = Crossspectrum(lc1, lc1, norm='none')
+            self.pds2 = Crossspectrum(lc2, lc2, norm='none')
 
     def _make_crossspectrum(self, lc1, lc2):
 
@@ -254,7 +258,7 @@ class Crossspectrum(object):
 
         return freqs[freqs > 0], cross
 
-    def rebin(self, df, method="mean"):
+    def rebin(self, df=None, f=None, method="mean"):
         """
         Rebin the cross spectrum to a new frequency resolution df.
 
@@ -263,20 +267,31 @@ class Crossspectrum(object):
         df: float
             The new frequency resolution
 
+        Other Parameters
+        ----------------
+        f: float
+            the rebin factor. If specified, it substitutes df with f*self.df
+
         Returns
         -------
         bin_cs = Crossspectrum object
             The newly binned cross spectrum
         """
 
+        if f is None and df is None:
+            raise ValueError('You need to specify at least one between f and '
+                             'df')
+        elif f is not None:
+            df = f * self.df
+
         # rebin cross spectrum to new resolution
         binfreq, bincs, binerr, step_size = \
             rebin_data(self.freq, self.power, df, self.power_err,
-                       method=method)
+                       method=method, dx=self.df)
 
         # make an empty cross spectrum object
         # note: syntax deliberate to work with subclass Powerspectrum
-        bin_cs = self.__class__()
+        bin_cs = copy.copy(self)
 
         # store the binned periodogram in the new object
         bin_cs.freq = binfreq
@@ -287,6 +302,23 @@ class Crossspectrum(object):
         bin_cs.nphots1 = self.nphots1
         bin_cs.power_err = binerr
 
+        if hasattr(self, 'unnorm_power'):
+            _, binpower_unnorm, _, _ = \
+                rebin_data(self.freq, self.unnorm_power, df,
+                           method=method, dx=self.df)
+
+            bin_cs.unnorm_power = binpower_unnorm
+
+        if hasattr(self, 'cs_all'):
+            cs_all = []
+            for c in self.cs_all:
+                cs_all.append(c.rebin(df=df, f=f, method=method))
+            bin_cs.cs_all = cs_all
+        if hasattr(self, 'pds1'):
+            bin_cs.pds1 = self.pds1.rebin(df=df, f=f, method=method)
+        if hasattr(self, 'pds2'):
+            bin_cs.pds2 = self.pds2.rebin(df=df, f=f, method=method)
+
         try:
             bin_cs.nphots2 = self.nphots2
         except AttributeError:
@@ -381,39 +413,42 @@ class Crossspectrum(object):
             frequency bin
         """
 
-        minfreq = self.freq[1] * 0.5  # frequency to start from
-        maxfreq = self.freq[-1]  # maximum frequency to end
-        binfreq = [minfreq, minfreq + self.df]  # first
-        df = self.freq[1]  # the frequency resolution of the first bin
-
-        # until we reach the maximum frequency, increase the width of each
-        # frequency bin by f
-        while binfreq[-1] <= maxfreq:
-            binfreq.append(binfreq[-1] + df*(1.0+f))
-            df = binfreq[-1] - binfreq[-2]
-
-        # compute the mean of the powers that fall into each new frequency bin.
-        # we cast to np.double due to scipy's bad handling of longdoubles
-        binpower, bin_edges, binno = scipy.stats.binned_statistic(
-            self.freq.astype(np.double), self.power.astype(np.double),
-            statistic="mean", bins=binfreq)
-
-        binpower_err, bin_edges, binno = scipy.stats.binned_statistic(
-            self.freq.astype(np.double), self.power_err.astype(np.double),
-            statistic=_root_squared_mean, bins=binfreq)
-
-        # compute the number of powers in each frequency bin
-        nsamples = np.array([len(binno[np.where(binno == i)[0]])
-                             for i in range(np.max(binno))])
+        binfreq, binpower, binpower_err, nsamples = \
+            rebin_data_log(self.freq, self.power, f,
+                           y_err=self.power_err, dx=self.df)
 
         # the frequency resolution
         df = np.diff(binfreq)
 
         # shift the lower bin edges to the middle of the bin and drop the
         # last right bin edge
-        binfreq = binfreq[:-1] + df/2
+        binfreq = binfreq[:-1] + df / 2
+
+        new_spec = copy.copy(self)
+        new_spec.freq = binfreq
+        new_spec.power = binpower
+        new_spec.power_err = binpower_err
+        new_spec.m = nsamples * self.m
+
+
+        if hasattr(self, 'unnorm_power'):
+            _, binpower_unnorm, _, _ = \
+                rebin_data_log(self.freq, self.unnorm_power, f, dx=self.df)
+
+            new_spec.unnorm_power = binpower_unnorm
 
-        return binfreq, binpower, binpower_err, nsamples
+        if hasattr(self, 'pds1'):
+            new_spec.pds1 = self.pds1.rebin_log(f)
+        if hasattr(self, 'pds2'):
+            new_spec.pds2 = self.pds2.rebin_log(f)
+
+        if hasattr(self, 'cs_all'):
+            cs_all = []
+            for c in self.cs_all:
+                cs_all.append(c.rebin_log(f))
+            new_spec.cs_all = cs_all
+
+        return new_spec
 
     def coherence(self):
         """
@@ -433,10 +468,8 @@ class Crossspectrum(object):
         """
         # this computes the averaged power spectrum, but using the
         # cross spectrum code to avoid circular imports
-        ps1 = Crossspectrum(self.lc1, self.lc1)
-        ps2 = Crossspectrum(self.lc2, self.lc2)
 
-        return self.unnorm_power.real/(ps1.unnorm_power.real * ps2.unnorm_power.real)
+        return self.unnorm_power.real/(self.pds1.power.real * self.pds2.power.real)
 
     def _phase_lag(self):
         """Return the fourier phase lag of the cross spectrum."""
@@ -542,6 +575,16 @@ class AveragedCrossspectrum(Crossspectrum):
 
         return
 
+    def _make_auxil_pds(self, lc1, lc2):
+        # A way to say that this is actually not a power spectrum
+        if lc1 is not lc2 and isinstance(lc1, Lightcurve):
+            self.pds1 = AveragedCrossspectrum(lc1, lc1,
+                                              segment_size=self.segment_size,
+                                              norm='none', gti=lc1.gti)
+            self.pds2 = AveragedCrossspectrum(lc2, lc2,
+                                              segment_size=self.segment_size,
+                                              norm='none', gti=lc2.gti)
+
     def _make_segment_spectrum(self, lc1, lc2, segment_size):
 
         # TODO: need to update this for making cross spectra.
@@ -685,7 +728,7 @@ class AveragedCrossspectrum(Crossspectrum):
         .. [1] http://iopscience.iop.org/article/10.1086/310430/pdf
 
         """
-        if self.m < 50:
+        if np.any(self.m < 50):
             simon("Number of segments used in averaging is "
                   "significantly low. The result might not follow the "
                   "expected statistical distributions.")
@@ -698,22 +741,9 @@ class AveragedCrossspectrum(Crossspectrum):
         unnorm_power_avg /= self.m
         num = np.absolute(unnorm_power_avg)**2
 
-        # this computes the averaged power spectrum, but using the
-        # cross spectrum code to avoid circular imports
-        aps1 = AveragedCrossspectrum(self.lc1, self.lc1,
-                                     segment_size=self.segment_size)
-        aps2 = AveragedCrossspectrum(self.lc2, self.lc2,
-                                     segment_size=self.segment_size)
-
-        unnorm_powers_avg_1 = np.zeros_like(aps1.cs_all[0].unnorm_power.real)
-        for ps in aps1.cs_all:
-            unnorm_powers_avg_1 += ps.unnorm_power.real
-        unnorm_powers_avg_1 /= aps1.m
-
-        unnorm_powers_avg_2 = np.zeros_like(aps2.cs_all[0].unnorm_power.real)
-        for ps in aps2.cs_all:
-            unnorm_powers_avg_2 += ps.unnorm_power.real
-        unnorm_powers_avg_2 /= aps2.m
+        # The normalization was 'none'!
+        unnorm_powers_avg_1 = self.pds1.power.real
+        unnorm_powers_avg_2 = self.pds2.power.real
 
         coh = num / (unnorm_powers_avg_1 * unnorm_powers_avg_2)
 
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index 6ea063ed..6e98da9f 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -504,7 +504,7 @@ class Lightcurve(object):
 
         return Lightcurve(time, counts, gti=gti, mjdref=mjdref, dt=dt)
 
-    def rebin(self, dt_new, method='sum'):
+    def rebin(self, dt_new=None, f=None, method='sum'):
         """
         Rebin the light curve to a new time resolution. While the new
         resolution need not be an integer multiple of the previous time
@@ -521,6 +521,11 @@ class Lightcurve(object):
             This keyword argument sets whether the counts in the new bins
             should be summed or averaged.
 
+        Other Parameters
+        ----------------
+        f: float
+            the rebin factor. If specified, it substitutes dt_new with
+            f*self.dt
 
         Returns
         -------
@@ -528,6 +533,12 @@ class Lightcurve(object):
             The :class:`Lightcurve` object with the new, binned light curve.
         """
 
+        if f is None and dt_new is None:
+            raise ValueError('You need to specify at least one between f and '
+                             'dt_new')
+        elif f is not None:
+            dt_new = f * self.dt
+
         if dt_new < self.dt:
             raise ValueError("New time resolution must be larger than "
                              "old time resolution!")
@@ -925,7 +936,7 @@ class Lightcurve(object):
             new_lc = Lightcurve(self.time[start:stop], self.counts[start:stop],
                                 err=self.counts_err[start:stop],
                                 mjdref=self.mjdref, gti=[self.gti[i]],
-                                dt=self.dt)
+                                dt=self.dt, err_dist=self.err_dist)
             list_of_lcs.append(new_lc)
 
         return list_of_lcs
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index cedf7b15..0cd61a5d 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -176,8 +176,8 @@ class Powerspectrum(Crossspectrum):
         Crossspectrum.__init__(self, lc1=lc, lc2=lc, norm=norm, gti=gti)
         self.nphots = self.nphots1
 
-    def rebin(self, df, method="mean"):
-        bin_ps = Crossspectrum.rebin(self, df=df, method=method)
+    def rebin(self, df=None, f=None, method="mean"):
+        bin_ps = Crossspectrum.rebin(self, df=df, f=f, method=method)
         bin_ps.nphots = bin_ps.nphots1
 
         return bin_ps
diff --git a/stingray/utils.py b/stingray/utils.py
index 2dfdf50a..74868918 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -29,6 +29,10 @@ except ImportError:
             return wrapped_f
 
 
+def _root_squared_mean(array):
+    return np.sqrt(np.sum(array**2)) / len(array)
+
+
 def simon(message, **kwargs):
     """The Statistical Interpretation MONitor.
 
@@ -47,7 +51,7 @@ def simon(message, **kwargs):
     warnings.warn("SIMON says: {0}".format(message), **kwargs)
 
 
-def rebin_data(x, y, dx_new, yerr=None, method='sum'):
+def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
 
     """Rebin some data to an arbitrary new data resolution. Either sum
     the data points in the new bins or average them.
@@ -60,16 +64,21 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum'):
     y: iterable
         The independent variable to be binned
 
-    yerr: iterable, optional
-        The uncertainties of y, to be propagated during binning.
-
     dx_new: float
         The new resolution of the dependent variable x
 
+    Other parameters
+    ----------------
+    yerr: iterable, optional
+        The uncertainties of y, to be propagated during binning.
+
     method: {"sum" | "average" | "mean"}, optional, default "sum"
         The method to be used in binning. Either sum the samples y in
         each new bin of x, or take the arithmetic mean.
 
+    dx: float
+        The old resolution (otherwise, calculated from median diff)
+
 
     Returns
     -------
@@ -89,7 +98,7 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum'):
     y = np.asarray(y)
     yerr = np.asarray(assign_value_if_none(yerr, np.zeros_like(y)))
 
-    dx_old = x[1] - x[0]
+    dx_old = assign_value_if_none(dx, np.median(np.diff(x)))
 
     if dx_new < dx_old:
         raise ValueError("New frequency resolution must be larger than "
@@ -148,6 +157,83 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum'):
     return xbin, ybin, ybinerr, step_size
 
 
+def rebin_data_log(x, y, f, y_err=None, dx=None):
+    """Logarithmic rebin of the periodogram.
+
+    The new frequency depends on the previous frequency modified by a factor f:
+
+    dnu_j = dnu_{j-1}*(1+f)
+
+    Parameters
+    ----------
+    x: iterable
+        The dependent variable with some resolution dx_old = x[1]-x[0]
+
+    y: iterable
+        The independent variable to be binned
+
+    f: float
+        The factor of increase of each bin wrt the previous one.
+
+    Other Parameters
+    ----------------
+    yerr: iterable, optional
+        The uncertainties of y, to be propagated during binning.
+
+    method: {"sum" | "average" | "mean"}, optional, default "sum"
+        The method to be used in binning. Either sum the samples y in
+        each new bin of x, or take the arithmetic mean.
+
+    dx: float, optional
+        The binning step of the initial xs
+
+    Returns
+    -------
+    xbin: numpy.ndarray
+        The midpoints of the new bins in x
+
+    ybin: numpy.ndarray
+        The binned quantity y
+
+    ybin_err: numpy.ndarray
+        The uncertainties of the binned values of y.
+
+    step_size: float
+        The size of the binning step
+    """
+    import scipy
+    dx_init = assign_value_if_none(dx, np.median(np.diff(x)))
+    y = np.asarray(y)
+    y_err = np.asarray(assign_value_if_none(y_err, np.zeros_like(y)))
+
+    minx = x[1] * 0.5  # frequency to start from
+    maxx = x[-1]  # maximum frequency to end
+    binx = [minx, minx + dx_init]  # first
+    dx = x[1]  # the frequency resolution of the first bin
+
+    # until we reach the maximum frequency, increase the width of each
+    # frequency bin by f
+    while binx[-1] <= maxx:
+        binx.append(binx[-1] + dx*(1.0+f))
+        dx = binx[-1] - binx[-2]
+
+    # compute the mean of the ys that fall into each new frequency bin.
+    # we cast to np.double due to scipy's bad handling of longdoubles
+    biny, bin_edges, binno = scipy.stats.binned_statistic(
+        x.astype(np.double), y.astype(np.double),
+        statistic="mean", bins=binx)
+
+    biny_err, bin_edges, binno = scipy.stats.binned_statistic(
+        x.astype(np.double), y_err.astype(np.double),
+        statistic=_root_squared_mean, bins=binx)
+
+    # compute the number of powers in each frequency bin
+    nsamples = np.array([len(binno[np.where(binno == i)[0]])
+                         for i in range(np.max(binno))])
+
+    return binx, biny, biny_err, nsamples
+
+
 def assign_value_if_none(value, default):
     return default if value is None else value
 


Test patch:
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index b1cb7be1..897adf56 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -5,6 +5,7 @@ import warnings
 from stingray import Lightcurve, AveragedPowerspectrum
 from stingray import Crossspectrum, AveragedCrossspectrum, coherence
 from stingray import StingrayError
+import copy
 
 np.random.seed(20160528)
 
@@ -51,10 +52,13 @@ class TestCoherence(object):
         assert np.abs(np.mean(coh)) < 1
 
     def test_high_coherence(self):
+        import copy
         t = np.arange(1280)
         a = np.random.poisson(100, len(t))
         lc = Lightcurve(t, a)
-        c = AveragedCrossspectrum(lc, lc, 128)
+        lc2 = Lightcurve(t, copy.copy(a))
+        c = AveragedCrossspectrum(lc, lc2, 128)
+
         coh, _ = c.coherence()
         np.testing.assert_almost_equal(np.mean(coh).real, 1.0)
 
@@ -121,6 +125,26 @@ class TestCrossspectrum(object):
         with pytest.raises(StingrayError):
             cs = Crossspectrum(self.lc1, lc_)
 
+    def test_make_crossspectrum_diff_lc_stat(self):
+        lc_ = copy.copy(self.lc1)
+        lc_.err_dist = 'gauss'
+
+        with pytest.warns(UserWarning) as record:
+            cs = Crossspectrum(self.lc1, lc_)
+        assert np.any(["different statistics" in r.message.args[0]
+                       for r in record])
+
+    def test_make_crossspectrum_bad_lc_stat(self):
+        lc1 = copy.copy(self.lc1)
+        lc1.err_dist = 'gauss'
+        lc2 = copy.copy(self.lc1)
+        lc2.err_dist = 'gauss'
+
+        with pytest.warns(UserWarning) as record:
+            cs = Crossspectrum(lc1, lc2)
+        assert np.any(["is not poisson" in r.message.args[0]
+                       for r in record])
+
     def test_make_crossspectrum_diff_dt(self):
         counts = np.array([1]*10000)
         time = np.linspace(0.0, 2.0, 10000)
@@ -136,10 +160,18 @@ class TestCrossspectrum(object):
     def test_rebin(self):
         new_cs = self.cs.rebin(df=1.5)
         assert new_cs.df == 1.5
+        new_cs.time_lag()
+
+    def test_rebin_factor(self):
+        new_cs = self.cs.rebin(f=1.5)
+        assert new_cs.df == self.cs.df * 1.5
+        new_cs.time_lag()
 
     def test_rebin_log(self):
         # For now, just verify that it doesn't crash
-        _ = self.cs.rebin_log(f=0.01)
+        new_cs = self.cs.rebin_log(f=0.1)
+        assert type(new_cs) == type(self.cs)
+        new_cs.time_lag()
 
     def test_norm_leahy(self):
         cs = Crossspectrum(self.lc1, self.lc2, norm='leahy')
@@ -305,8 +337,6 @@ class TestAveragedCrossspectrum(object):
 
             assert len(coh[0]) == 4999
             assert len(coh[1]) == 4999
-
-            assert len(w) == 3
             assert issubclass(w[-1].category, UserWarning)
 
     def test_failure_when_normalization_not_recognized(self):
@@ -315,6 +345,22 @@ class TestAveragedCrossspectrum(object):
                                             segment_size=1,
                                             norm="wrong")
 
+    def test_rebin(self):
+        new_cs = self.cs.rebin(df=1.5)
+        assert new_cs.df == 1.5
+        new_cs.time_lag()
+
+    def test_rebin_factor(self):
+        new_cs = self.cs.rebin(f=1.5)
+        assert new_cs.df == self.cs.df * 1.5
+        new_cs.time_lag()
+
+    def test_rebin_log(self):
+        # For now, just verify that it doesn't crash
+        new_cs = self.cs.rebin_log(f=0.1)
+        assert type(new_cs) == type(self.cs)
+        new_cs.time_lag()
+
     def test_timelag(self):
         from ..simulator.simulator import Simulator
         dt = 0.1
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index ce38a830..47bbd00f 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -610,6 +610,15 @@ class TestLightcurveRebin(object):
             self.lc.counts[0]*dt_new/self.lc.dt
         assert np.allclose(lc_binned.counts, counts_test)
 
+    def test_rebin_even_factor(self):
+        f = 200
+        dt_new = f * self.lc.dt
+        lc_binned = self.lc.rebin(f=f)
+        assert np.isclose(dt_new, f * self.lc.dt)
+        counts_test = np.zeros_like(lc_binned.time) + \
+            self.lc.counts[0]*dt_new/self.lc.dt
+        assert np.allclose(lc_binned.counts, counts_test)
+
     def test_rebin_odd(self):
         dt_new = 1.5
         lc_binned = self.lc.rebin(dt_new)
@@ -619,6 +628,15 @@ class TestLightcurveRebin(object):
             self.lc.counts[0]*dt_new/self.lc.dt
         assert np.allclose(lc_binned.counts, counts_test)
 
+    def test_rebin_odd_factor(self):
+        f = 100.5
+        dt_new = f * self.lc.dt
+        lc_binned = self.lc.rebin(f=f)
+        assert np.isclose(dt_new, f * self.lc.dt)
+        counts_test = np.zeros_like(lc_binned.time) + \
+            self.lc.counts[0]*dt_new/self.lc.dt
+        assert np.allclose(lc_binned.counts, counts_test)
+
     def rebin_several(self, dt):
         lc_binned = self.lc.rebin(dt)
         assert len(lc_binned.time) == len(lc_binned.counts)
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index 41fdbdc5..87b488be 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -193,7 +193,7 @@ class TestPowerspectrum(object):
         had better be 'method'
         """
         ps = Powerspectrum(lc=self.lc, norm="Leahy")
-        assert ps.rebin.__defaults__[0] == "mean"
+        assert ps.rebin.__defaults__[2] == "mean"
 
     @pytest.mark.parametrize('df', [2, 3, 5, 1.5, 1, 85])
     def test_rebin(self, df):
@@ -344,32 +344,50 @@ class TestAveragedPowerspectrum(object):
         segment_size = 0.5
         assert AveragedPowerspectrum(lc_all, segment_size)
 
-        def rebin_several_averagedps(self, df):
-            """
-            TODO: Not sure how to write tests for the rebin method!
-            """
-
-            aps = AveragedPowerspectrum(lc=self.lc, segment_size=1,
-                                        norm="Leahy")
-            bin_aps = aps.rebin(df)
-            assert np.isclose(bin_aps.freq[1] - bin_aps.freq[0], bin_aps.df,
-                              atol=1e-4, rtol=1e-4)
-            assert np.isclose(bin_aps.freq[0],
-                              (aps.freq[0] - aps.df * 0.5 + bin_aps.df * 0.5),
-                              atol=1e-4, rtol=1e-4)
-
-        def test_rebin_averagedps(self):
-            df_all = [2, 3, 5, 1.5, 10, 45]
-            for df in df_all:
-                yield self.rebin_several, df
-
-        def test_rebin_with_invalid_type_attribute(self):
-            new_df = 2
-            aps = AveragedPowerspectrum(lc=self.lc, segment_size=1,
-                                        norm='leahy')
-            aps.type = 'invalid_type'
-            with pytest.raises(AttributeError):
-                assert aps.rebin(df=new_df)
+    @pytest.mark.parametrize('df', [2, 3, 5, 1.5, 1, 85])
+    def test_rebin(self, df):
+        """
+        TODO: Not sure how to write tests for the rebin method!
+        """
+
+        aps = AveragedPowerspectrum(lc=self.lc, segment_size=1,
+                                    norm="Leahy")
+        bin_aps = aps.rebin(df)
+        assert np.isclose(bin_aps.freq[1]-bin_aps.freq[0], bin_aps.df,
+                          atol=1e-4, rtol=1e-4)
+        assert np.isclose(bin_aps.freq[0],
+                          (aps.freq[0]-aps.df*0.5+bin_aps.df*0.5),
+                          atol=1e-4, rtol=1e-4)
+
+    @pytest.mark.parametrize('f', [2, 3, 5, 1.5, 1, 85])
+    def test_rebin_factor(self, f):
+        """
+        TODO: Not sure how to write tests for the rebin method!
+        """
+
+        aps = AveragedPowerspectrum(lc=self.lc, segment_size=1,
+                                    norm="Leahy")
+        bin_aps = aps.rebin(f=f)
+        assert np.isclose(bin_aps.freq[1]-bin_aps.freq[0], bin_aps.df,
+                          atol=1e-4, rtol=1e-4)
+        assert np.isclose(bin_aps.freq[0],
+                          (aps.freq[0]-aps.df*0.5+bin_aps.df*0.5),
+                          atol=1e-4, rtol=1e-4)
+
+    @pytest.mark.parametrize('df', [0.01, 0.1])
+    def test_rebin_log(self, df):
+        # For now, just verify that it doesn't crash
+        aps = AveragedPowerspectrum(lc=self.lc, segment_size=1,
+                                    norm="Leahy")
+        bin_aps = aps.rebin_log(df)
+
+    def test_rebin_with_invalid_type_attribute(self):
+        new_df = 2
+        aps = AveragedPowerspectrum(lc=self.lc, segment_size=1,
+                                    norm='leahy')
+        aps.type = 'invalid_type'
+        with pytest.raises(AttributeError):
+            assert aps.rebin(df=new_df)
 
     def test_list_with_nonsense_component(self):
         n_lcs = 10