Repository: StingraySoftware/stingray
Pull request #306 (2018-02-01 17:31:38)

Can now run classical_significances on log-binned periodograms
I realized it wasn't possible to run the method `classical_significances` on log-binned periodograms, because it expected `self.m` to be a single integer, rather than a `numpy.ndarray`. 

In the process of enabling that functionality, I discovered a bug in the way `self.m` was calculated for the log-binned periodogram and fixed that, too.


Patch:
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index e49eaf74..df370741 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -314,10 +314,15 @@ class Powerspectrum(Crossspectrum):
             raise ValueError("This method only works on "
                              "Leahy-normalized power spectra!")
 
-        # calculate p-values for all powers
-        # leave out zeroth power since it just encodes the number of photons!
-        pv = np.array([classical_pvalue(power, self.m)
-                       for power in self.power])
+        if np.size(self.m) == 1:
+            # calculate p-values for all powers
+            # leave out zeroth power since it just encodes the number of photons!
+            pv = np.array([classical_pvalue(power, self.m)
+                           for power in self.power])
+        else:
+            pv = np.array([classical_pvalue(power, m)
+                           for power, m in zip(self.power, self.m)])
+
 
         # if trial correction is used, then correct the threshold for
         # the number of powers in the power spectrum
diff --git a/stingray/utils.py b/stingray/utils.py
index 48e7ae61..996da68f 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -222,14 +222,21 @@ def rebin_data_log(x, y, f, y_err=None, dx=None):
     step_size: float
         The size of the binning step
     """
+
     dx_init = assign_value_if_none(dx, np.median(np.diff(x)))
+    x = np.asarray(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
+    if x.shape[0] != y.shape[0]:
+        raise ValueError("x and y must be of the same length!")
+    if y.shape[0] != y_err.shape[0]:
+        raise ValueError("y and y_err must be of the same length!")
+
+    minx = x[0] * 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
+    dx = dx_init  # the frequency resolution of the first bin
 
     # until we reach the maximum frequency, increase the width of each
     # frequency bin by f
@@ -237,6 +244,8 @@ def rebin_data_log(x, y, f, y_err=None, dx=None):
         binx.append(binx[-1] + dx * (1.0 + f))
         dx = binx[-1] - binx[-2]
 
+    binx = np.asarray(binx)
+
     real = y.real
     real_err = y_err.real
     # compute the mean of the ys that fall into each new frequency bin.
@@ -267,7 +276,7 @@ def rebin_data_log(x, y, f, y_err=None, dx=None):
 
     # 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))])
+                         for i in range(1, np.max(binno)+1, 1)])
 
     return binx, biny, biny_err, nsamples
 


Test patch:
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index 5c4803f6..0c99b952 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -109,8 +109,6 @@ class TestPowerspectrum(object):
         ps = Powerspectrum(lc=self.lc, norm="frac")
         ps_int = np.sum(ps.power[:-1] * ps.df) + ps.power[-1] * ps.df / 2
         std_lc = np.var(self.lc.counts) / np.mean(self.lc.counts) ** 2
-        print(ps_int)
-        print(std_lc)
         assert np.isclose(ps_int, std_lc, atol=0.01, rtol=0.01)
 
     def test_fractional_rms_in_frac_norm_is_consistent(self):
@@ -287,6 +285,14 @@ class TestPowerspectrum(object):
                                           trial_correction=True)
         assert np.size(pval) == 0
 
+
+    def test_classical_significances_with_logbinned_psd(self):
+        ps = Powerspectrum(lc=self.lc, norm="leahy")
+        ps_log = ps.rebin_log()
+        pval = ps_log.classical_significances(threshold=1.1, trial_correction=False)
+
+        assert len(pval[0]) == len(ps_log.power)
+
     def test_pvals_is_numpy_array(self):
         ps = Powerspectrum(lc=self.lc, norm="leahy")
         # change the powers so that just one exceeds the threshold
diff --git a/stingray/tests/test_utils.py b/stingray/tests/test_utils.py
index 2c666973..863873e1 100644
--- a/stingray/tests/test_utils.py
+++ b/stingray/tests/test_utils.py
@@ -65,6 +65,90 @@ class TestRebinData(object):
                              method='not_allowed_method')
 
 
+class TestRebinDataLog(object):
+
+    @classmethod
+    def setup_class(cls):
+        cls.dx = 1
+        cls.xmax = 21
+        cls.xmin = 1
+        cls.x = np.arange(cls.xmin, cls.xmax, cls.dx)
+        cls.y = np.arange(cls.xmin, cls.xmax, cls.dx)
+        cls.y_err = np.ones_like(cls.y)
+
+        cls.true_bins = np.array([ 1., 1.1, 1.21, 1.331, 1.4641, 1.61051, 1.771561,
+                                     1.9487171, 2.14358881, 2.35794769, 2.59374246, 2.85311671])
+
+        cls.true_bin_edges = np.array([0.5, 1.5, 2.6000000000000001, 3.81, 5.141, 6.6051, 8.21561,
+                                         9.987171, 11.9358881, 14.07947691, 16.437424601, 19.0311670611,
+                                         21.88428376721])
+
+        cls.true_values = np.array([1. , 2., 3., 4.5, 6., 7.5, 9., 10.5, 13., 15.5, 18., 20.])
+        cls.true_nsamples = np.array([1, 1, 1, 2, 1, 2, 1, 2, 3, 2, 3, 1])
+
+        cls.f = 0.1
+
+    def test_rebin_data_log_runs(self):
+        _, _, _, _ = utils.rebin_data_log(self.x, self.y, self.f, y_err=self.y_err, dx=self.dx)
+
+    def test_method_fails_if_x_and_y_of_unequal_length(self):
+        with pytest.raises(ValueError):
+            _, _, _, _ = utils.rebin_data_log(self.x[1:], self.y, self.f, y_err=self.y_err, dx=self.dx)
+
+    def test_method_fails_if_y_and_yerr_of_unequal_length(self):
+        with pytest.raises(ValueError):
+            _, _, _, _ = utils.rebin_data_log(self.x, self.y, self.f, y_err=self.y_err[1:], dx=self.dx)
+
+
+    def test_all_outputs_have_the_same_dimension_except_binx(self):
+        binx, biny, binyerr, nsamples = utils.rebin_data_log(self.x, self.y, self.f, y_err=self.y_err, dx=self.dx)
+
+        # binx describes the bin _edges_ rather than midpoints, so has one more entry
+        # than biny and the rest
+        assert binx.shape[0] == biny.shape[0]+1
+        assert biny.shape[0] == binyerr.shape[0]
+        assert binyerr.shape[0] == nsamples.shape[0]
+
+    def test_binning_works_correctly(self):
+        binx, _, _, _ = utils.rebin_data_log(self.x, self.y, self.f, y_err=self.y_err, dx=self.dx)
+
+        assert np.allclose(np.diff(binx), self.true_bins)
+
+    def test_bin_edges_are_correct(self):
+        binx, _, _, _ = utils.rebin_data_log(self.x, self.y, self.f, y_err=self.y_err, dx=self.dx)
+
+        assert np.allclose(binx, self.true_bin_edges)
+
+    def test_bin_values_are_correct(self):
+        _, biny, _, _ = utils.rebin_data_log(self.x, self.y, self.f, y_err=self.y_err, dx=self.dx)
+
+        assert np.allclose(biny, self.true_values)
+
+    def test_nsamples_are_correctly_calculated(self):
+        _, _, _, nsamples = utils.rebin_data_log(self.x, self.y, self.f, y_err=self.y_err, dx=self.dx)
+
+        assert np.allclose(nsamples, self.true_nsamples)
+
+    def test_method_works_on_complex_numbers(self):
+        re = np.arange(self.xmin, self.xmax, self.dx)
+        im = np.arange(self.xmin, self.xmax, self.dx)
+
+        y = np.zeros(re.shape[0], dtype=np.complex)
+        yerr = np.zeros(re.shape[0], dtype=np.complex)
+
+        for k, (r, i) in enumerate(zip(re, im)):
+            y[k] = r + i * 1j
+            yerr[k] = r + i * 1j
+
+        real_binned = np.zeros(self.true_values.shape[0], dtype=np.complex)
+
+        for i in range(self.true_values.shape[0]):
+            real_binned[i] = self.true_values[i] + self.true_values[i] * 1j
+
+        _, biny, _, _ = utils.rebin_data_log(self.x, y, self.f, y_err=yerr, dx=self.dx)
+
+        assert np.allclose(biny, real_binned)
+
 class TestUtils(object):
 
     def test_optimal_bin_time(self):