Repository: StingraySoftware/stingray
Pull request #279 (2017-11-14 18:43:56)

Fixed rebinning issues for complex numbers
This should fix the issue whereby the `rebin_data_log` in `utils.py` casts complex numbers to real ones, which means we loose the imaginary parts of the cross spectrum when rebinning. 

It now computes the real and imaginary parts of the rebinned cross spectrum separately and combines them back in the end into a set of complex values. 


Patch:
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 8e56b129..5366648f 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -469,7 +469,6 @@ class Crossspectrum(object):
 
     def _phase_lag(self):
         """Return the fourier phase lag of the cross spectrum."""
-
         return np.angle(self.unnorm_power)
 
     def time_lag(self):
@@ -758,9 +757,18 @@ class AveragedCrossspectrum(Crossspectrum):
         """Calculate time lag and uncertainty.
 
         Formula from Bendat & Piersol 1986
+
+        Returns
+        -------
+        lag : np.ndarray
+            The time lag
+
+        lag_err : np.ndarray
+            The uncertainty in the time lag
+
         """
         lag = super(AveragedCrossspectrum, self).time_lag()
         coh, uncert = self.coherence()
         dum = (1. - coh) / (2. * coh)
         lag_err = np.sqrt(dum / self.m) / (2 * np.pi * self.freq)
-        return (lag, lag_err)
+        return lag, lag_err
diff --git a/stingray/utils.py b/stingray/utils.py
index c56dd7e5..c303b8aa 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -7,6 +7,7 @@ from six import string_types
 
 import warnings
 import numpy as np
+import scipy
 
 # If numba is installed, import jit. Otherwise, define an empty decorator with
 # the same name.
@@ -222,7 +223,6 @@ def rebin_data_log(x, y, f, y_err=None, dx=None):
     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)))
@@ -238,16 +238,34 @@ 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]
 
+    real = y.real
+    real_err = y_err.real
     # 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),
+        x.astype(np.double), real.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),
+        x.astype(np.double), real_err.astype(np.double),
         statistic=_root_squared_mean, bins=binx)
 
+    if isinstance(y[0], np.complex):
+        imag = y.imag
+        biny_imag, bin_edges, binno = scipy.stats.binned_statistic(
+            x.astype(np.double), imag.astype(np.double),
+            statistic="mean", bins=binx)
+
+        biny = biny + 1j*biny_imag
+
+    if isinstance(y_err[0], np.complex):
+        imag_err = y_err.imag
+
+        biny_err_imag, bin_edges, binno = scipy.stats.binned_statistic(
+            x.astype(np.double), imag_err.astype(np.double),
+            statistic=_root_squared_mean, bins=binx)
+        biny_err = biny_err + 1j * biny_err_imag
+
     # 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))])


Test patch:
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index 64e326fd..a9df0bda 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -376,6 +376,16 @@ class TestAveragedCrossspectrum(object):
         assert type(new_cs) == type(self.cs)
         new_cs.time_lag()
 
+    def test_rebin_log_returns_complex_values(self):
+        # For now, just verify that it doesn't crash
+        new_cs = self.cs.rebin_log(f=0.1)
+        assert isinstance(new_cs.power[0], np.complex)
+
+    def test_rebin_log_returns_complex_errors(self):
+        # For now, just verify that it doesn't crash
+        new_cs = self.cs.rebin_log(f=0.1)
+        assert isinstance(new_cs.power_err[0], np.complex)
+
     def test_timelag(self):
         from ..simulator.simulator import Simulator
         dt = 0.1