Repository: StingraySoftware/stingray

Issue:
`_isattribute` in `io.py` not general enough?
I am moderately concerned that our `_isattribute(data)` function in `io.py is not general enough.

This function checks whether the data is a single number or an array, by running `isinstance` for `numpy.ndarray` and `list`. 

My feeling is that `isinstance` is perhaps not the best way to go about this? What, for example, if `data` is a tuple or a set? 

I'm wondering whether one can e.g. transform `data` into a `numpy.ndarray` and then run e.g. `np.size` on it instead? 

Discussion:

Patch:
diff --git a/setup.cfg b/setup.cfg
index 2200bd09..06301011 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -34,7 +34,7 @@ all =
     tbb
     numba
     jplephem
-    emcee>=3.0.0
+    emcee
     statsmodels
     corner
     lightcurve
@@ -64,6 +64,10 @@ doctest_plus = enabled
 text_file_format = rst
 addopts = --doctest-rst
 filterwarnings =
+    ignore:SIMON says. Errorbars on cross spectra are not thoroughly tested.:UserWarning
+    ignore:SIMON says. Stingray only uses poisson err_dist at the moment.:UserWarning
+    ignore:SIMON says. Looks like your lightcurve:UserWarning
+    ignore:Matplotlib is currently using agg, which is a:UserWarning
     ignore:Using or importing the ABCs from 'collections':DeprecationWarning
     ignore:unclosed file:ResourceWarning
     ignore:numpy.ufunc size changed:RuntimeWarning
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 9ca4cb14..42bad3fe 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -25,7 +25,7 @@ except ImportError:
 try:
     from pyfftw.interfaces.scipy_fftpack import fft, fftfreq
 except ImportError:
-    warnings.warn("pyfftw not installed. Using standard scipy fft")
+    warnings.warn("Using standard scipy fft")
     from scipy.fftpack import fft, fftfreq
 
 __all__ = [
@@ -614,7 +614,6 @@ class Crossspectrum(object):
         binfreq, bincs, binerr, step_size = \
             rebin_data(self.freq, self.power, df, self.power_err,
                        method=method, dx=self.df)
-
         # make an empty cross spectrum object
         # note: syntax deliberate to work with subclass Powerspectrum
         bin_cs = copy.copy(self)
@@ -995,6 +994,11 @@ class AveragedCrossspectrum(Crossspectrum):
     large_data : bool, default False
         Use only for data larger than 10**7 data points!! Uses zarr and dask for computation.
 
+    save_all : bool, default False
+        Save all intermediate PDSs used for the final average. Use with care.
+        This is likely to fill up your RAM on medium-sized datasets, and to
+        slow down the computation when rebinning.
+
     Attributes
     ----------
     freq: numpy.ndarray
@@ -1033,7 +1037,7 @@ class AveragedCrossspectrum(Crossspectrum):
 
     def __init__(self, data1=None, data2=None, segment_size=None, norm='none',
                  gti=None, power_type="real", silent=False, lc1=None, lc2=None,
-                 dt=None, large_data=False):
+                 dt=None, large_data=False, save_all=False):
 
         if lc1 is not None or lc2 is not None:
             warnings.warn("The lcN keywords are now deprecated. Use dataN "
@@ -1089,6 +1093,7 @@ class AveragedCrossspectrum(Crossspectrum):
 
         self.show_progress = not silent
         self.dt = dt
+        self.save_all = save_all
 
         if isinstance(data1, EventList):
             lengths = data1.gti[:, 1] - data1.gti[:, 0]
@@ -1128,14 +1133,14 @@ class AveragedCrossspectrum(Crossspectrum):
                                               segment_size=self.segment_size,
                                               norm='none', gti=self.gti,
                                               power_type=self.power_type,
-                                              dt=self.dt)
+                                              dt=self.dt, save_all=self.save_all)
             self.pds2 = AveragedCrossspectrum(lc2, lc2,
                                               segment_size=self.segment_size,
                                               norm='none', gti=self.gti,
                                               power_type=self.power_type,
-                                              dt=self.dt)
+                                              dt=self.dt, save_all=self.save_all)
 
-    def _make_segment_spectrum(self, lc1, lc2, segment_size):
+    def _make_segment_spectrum(self, lc1, lc2, segment_size, silent=False):
         """
         Split the light curves into segments of size ``segment_size``, and calculate a cross spectrum for
         each.
@@ -1148,6 +1153,11 @@ class AveragedCrossspectrum(Crossspectrum):
         segment_size : ``numpy.float``
             Size of each light curve segment to use for averaging.
 
+        Other parameters
+        ----------------
+        silent : bool, default False
+            Suppress progress bars
+
         Returns
         -------
         cs_all : list of :class:`Crossspectrum`` objects
@@ -1158,7 +1168,6 @@ class AveragedCrossspectrum(Crossspectrum):
 
         """
 
-        # TODO: need to update this for making cross spectra.
         assert isinstance(lc1, Lightcurve)
         assert isinstance(lc2, Lightcurve)
 
@@ -1199,7 +1208,7 @@ class AveragedCrossspectrum(Crossspectrum):
               "Please report any inconsistencies.")
 
         local_show_progress = show_progress
-        if not self.show_progress:
+        if not self.show_progress or silent:
             local_show_progress = lambda a: a
 
         for start_ind, end_ind in \
@@ -1259,47 +1268,49 @@ class AveragedCrossspectrum(Crossspectrum):
                 isinstance(lc2, Lightcurve):
 
             if self.type == "crossspectrum":
-                self.cs_all, nphots1_all, nphots2_all = \
+                cs_all, nphots1_all, nphots2_all = \
                     self._make_segment_spectrum(lc1, lc2, self.segment_size)
 
             elif self.type == "powerspectrum":
-                self.cs_all, nphots1_all = \
+                cs_all, nphots1_all = \
                     self._make_segment_spectrum(lc1, self.segment_size)
 
             else:
                 raise ValueError("Type of spectrum not recognized!")
 
         else:
-            self.cs_all, nphots1_all, nphots2_all = [], [], []
+            cs_all, nphots1_all, nphots2_all = [], [], []
 
             for lc1_seg, lc2_seg in local_show_progress(zip(lc1, lc2)):
                 if self.type == "crossspectrum":
                     cs_sep, nphots1_sep, nphots2_sep = \
                         self._make_segment_spectrum(lc1_seg, lc2_seg,
-                                                    self.segment_size)
+                                                    self.segment_size,
+                                                    silent=True)
                     nphots2_all.append(nphots2_sep)
                 elif self.type == "powerspectrum":
                     cs_sep, nphots1_sep = \
-                        self._make_segment_spectrum(lc1_seg, self.segment_size)
+                        self._make_segment_spectrum(lc1_seg, self.segment_size,
+                            silent=True)
 
                 else:
                     raise ValueError("Type of spectrum not recognized!")
-                self.cs_all.append(cs_sep)
+                cs_all.append(cs_sep)
                 nphots1_all.append(nphots1_sep)
 
-            self.cs_all = np.hstack(self.cs_all)
+            cs_all = np.hstack(cs_all)
             nphots1_all = np.hstack(nphots1_all)
 
             if self.type == "crossspectrum":
                 nphots2_all = np.hstack(nphots2_all)
 
-        m = len(self.cs_all)
+        m = len(cs_all)
         nphots1 = np.mean(nphots1_all)
 
-        power_avg = np.zeros_like(self.cs_all[0].power)
-        power_err_avg = np.zeros_like(self.cs_all[0].power_err)
-        unnorm_power_avg = np.zeros_like(self.cs_all[0].unnorm_power)
-        for cs in self.cs_all:
+        power_avg = np.zeros_like(cs_all[0].power)
+        power_err_avg = np.zeros_like(cs_all[0].power_err)
+        unnorm_power_avg = np.zeros_like(cs_all[0].unnorm_power)
+        for cs in cs_all:
             power_avg += cs.power
             unnorm_power_avg += cs.unnorm_power
             power_err_avg += (cs.power_err) ** 2
@@ -1308,14 +1319,16 @@ class AveragedCrossspectrum(Crossspectrum):
         power_err_avg = np.sqrt(power_err_avg) / m
         unnorm_power_avg /= np.float(m)
 
-        self.freq = self.cs_all[0].freq
+        self.freq = cs_all[0].freq
         self.power = power_avg
         self.unnorm_power = unnorm_power_avg
         self.m = m
         self.power_err = power_err_avg
-        self.df = self.cs_all[0].df
-        self.n = self.cs_all[0].n
+        self.df = cs_all[0].df
+        self.n = cs_all[0].n
         self.nphots1 = nphots1
+        if self.save_all:
+            self.cs_all = cs_all
 
         if self.type == "crossspectrum":
             self.nphots1 = nphots1
diff --git a/stingray/io.py b/stingray/io.py
index 40566161..25a787a9 100644
--- a/stingray/io.py
+++ b/stingray/io.py
@@ -990,7 +990,7 @@ def _isattribute(data):
         True if the data is a single number, False if it is an iterable.
     """
 
-    if isinstance(data, np.ndarray) or isinstance(data, list):
+    if isinstance(data, Iterable) and not isinstance(data, (str, bytes)):
         return False
     else:
         return True
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index cf4893ff..f651f806 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -1163,7 +1163,6 @@ class Lightcurve(object):
         >>> split_lc = lc.split(1.5)
 
         """
-
         # calculate the difference between time bins
         tdiff = np.diff(self.time)
         # find all distances between time bins that are larger than `min_gap`
@@ -1180,8 +1179,9 @@ class Lightcurve(object):
         gti = np.vstack([gti_start, gti_stop]).T
         if hasattr(self, 'gti') and self.gti is not None:
             gti = cross_two_gtis(self.gti, gti)
+        self.gti = gti
 
-        lc_split = self.split_by_gti(gti, min_points=min_points)
+        lc_split = self.split_by_gti(min_points=min_points)
         return lc_split
 
     def sort(self, reverse=False):
@@ -1302,17 +1302,7 @@ class Lightcurve(object):
         >>> count[2:4] = 1
         >>> lc = Lightcurve(time, count, dt=1)
         >>> lc.estimate_chunk_length(min_total_counts=3, min_time_bins=1)
-        3.0
-        >>> # A slightly more complex example
-        >>> dt=0.2
-        >>> time = np.arange(0, 1000, dt)
-        >>> counts = np.random.poisson(100, size=len(time))
-        >>> lc = Lightcurve(time, counts, dt=dt)
-        >>> lc.estimate_chunk_length(100, 2)
-        0.4
-        >>> min_total_bins = 40
-        >>> lc.estimate_chunk_length(100, 40)
-        8.0
+        4.0
         """
 
         rough_estimate = np.ceil(min_total_counts / self.meancounts) * self.dt
@@ -1327,7 +1317,7 @@ class Lightcurve(object):
             if mincounts >= min_total_counts:
                 keep_searching = False
             else:
-                chunk_length += self.dt 
+                chunk_length *= np.ceil(min_total_counts / mincounts) * self.dt
 
         return chunk_length
 
@@ -1382,7 +1372,7 @@ class Lightcurve(object):
                                               self.time,
                                               fraction_step=fraction_step,
                                               dt=self.dt)
-        start_times = self.time[start] - 0.5 * self.dt
+        start_times = self.time[start] - self.dt * 0.5
 
         # Remember that stop is one element above the last element, because
         # it's defined to be used in intervals start:stop
@@ -1658,7 +1648,7 @@ class Lightcurve(object):
                           err_dist=data['err_dist'],
                           mjdref=data['mjdref'])
 
-    def split_by_gti(self, gti=None, min_points=2):
+    def split_by_gti(self, min_points=2):
         """
         Split the current :class:`Lightcurve` object into a list of :class:`Lightcurve` objects, one
         for each continuous GTI segment as defined in the ``gti`` attribute.
@@ -1674,13 +1664,9 @@ class Lightcurve(object):
         list_of_lcs : list
             A list of :class:`Lightcurve` objects, one for each GTI segment
         """
-
-        if gti is None:
-            gti = self.gti
-
         list_of_lcs = []
 
-        start_bins, stop_bins = gti_border_bins(gti, self.time, self.dt)
+        start_bins, stop_bins = gti_border_bins(self.gti, self.time, self.dt)
         for i in range(len(start_bins)):
             start = start_bins[i]
             stop = stop_bins[i]
@@ -1694,7 +1680,7 @@ class Lightcurve(object):
             # Note: GTIs are consistent with default in this case!
             new_lc = Lightcurve(self.time[start:stop], self.counts[start:stop],
                                 err=self.counts_err[start:stop],
-                                mjdref=self.mjdref, gti=[gti[i]],
+                                mjdref=self.mjdref, gti=[self.gti[i]],
                                 dt=self.dt, err_dist=self.err_dist,
                                 skip_checks=True)
             list_of_lcs.append(new_lc)
diff --git a/stingray/modeling/parameterestimation.py b/stingray/modeling/parameterestimation.py
index f2f4f442..e1109faf 100644
--- a/stingray/modeling/parameterestimation.py
+++ b/stingray/modeling/parameterestimation.py
@@ -63,16 +63,12 @@ class OptimizationResults(object):
     res: instance of ``scipy.OptimizeResult``
         The object containing the results from a optimization run
 
+    Attributes
+    ----------
     neg : bool, optional, default ``True``
         A flag that sets whether the log-likelihood or negative log-likelihood
         is being used
 
-    log : a logging.getLogger() object, default None
-        You can pass a pre-defined object for logging, else a new 
-        logger will be instantiated
-
-    Attributes
-    ----------
     result : float
         The result of the optimization, i.e. the function value at the
         minimum that the optimizer found
@@ -129,26 +125,17 @@ class OptimizationResults(object):
     .. [bic] https://projecteuclid.org/euclid.aos/1176344136
 
     """
-    def __init__(self, lpost, res, neg=True, log=None):
+    def __init__(self, lpost, res, neg=True):
         self.neg = neg
         self.result = res.fun
         self.p_opt = np.atleast_1d(res.x)
         self.model = lpost.model
 
-        if log is None:
-            self.log = logging.getLogger('Fitting summary')
-            self.log.setLevel(logging.DEBUG)
-            if not self.log.handlers:
-                ch = logging.StreamHandler()
-                ch.setLevel(logging.DEBUG)
-                self.log.addHandler(ch)
-
         self._compute_covariance(lpost, res)
         self._compute_model(lpost)
         self._compute_criteria(lpost)
         self._compute_statistics(lpost)
 
-
     def _compute_covariance(self, lpost, res):
         """
         Compute the covariance of the parameters using inverse of the Hessian, i.e.
@@ -179,7 +166,7 @@ class OptimizationResults(object):
         else:
             if comp_hessian:
                 # calculate Hessian approximating with finite differences
-                self.log.info("Approximating Hessian with finite differences ...")
+                logging.info("Approximating Hessian with finite differences ...")
 
                 phess = approx_hess(np.atleast_1d(self.p_opt), lpost)
 
@@ -260,7 +247,7 @@ class OptimizationResults(object):
         self.ssd = np.sqrt(2.0*self.sexp)
         self.sobs = np.sum(lpost.y-self.mfit)
 
-    def print_summary(self, lpost):
+    def print_summary(self, lpost, log=None):
         """
         Print a useful summary of the fitting procedure to screen or
         a log file.
@@ -270,9 +257,19 @@ class OptimizationResults(object):
         lpost : instance of :class:`Posterior` or one of its subclasses
             The object containing the function that is being optimized
             in the regression
+
+        log : logging handler, optional, default None
+            A handler used for logging the output properly
         """
+        if log is None:
+            log = logging.getLogger('Fitting summary')
+            log.setLevel(logging.DEBUG)
 
-        self.log.info("The best-fit model parameters plus errors are:")
+            ch = logging.StreamHandler()
+            ch.setLevel(logging.DEBUG)
+            log.addHandler(ch)
+
+        log.info("The best-fit model parameters plus errors are:")
 
         fixed = [lpost.model.fixed[n] for n in lpost.model.param_names]
         tied = [lpost.model.tied[n] for n in lpost.model.param_names]
@@ -284,7 +281,7 @@ class OptimizationResults(object):
 
         all_parnames = [n for n in lpost.model.param_names]
         for i, par in enumerate(all_parnames):
-            self.log.info("{:3}) Parameter {:<20}: ".format(i, par))
+            log.info("{:3}) Parameter {:<20}: ".format(i, par))
 
             if par in parnames:
                 idx = parnames.index(par)
@@ -292,29 +289,29 @@ class OptimizationResults(object):
                 err_info = " (no error estimate)"
                 if self.err is not None:
                     err_info = " +/- {:<20.5f}".format(self.err[idx])
-                self.log.info("{:<20.5f}{} ".format(self.p_opt[idx], err_info))
-                self.log.info("[{:>10} {:>10}]".format(str(bounds[i][0]),
+                log.info("{:<20.5f}{} ".format(self.p_opt[idx], err_info))
+                log.info("[{:>10} {:>10}]".format(str(bounds[i][0]),
                                                str(bounds[i][1])))
             elif fixed[i]:
-                self.log.info("{:<20.5f} (Fixed) ".format(lpost.model.parameters[i]))
+                log.info("{:<20.5f} (Fixed) ".format(lpost.model.parameters[i]))
             elif tied[i]:
-                self.log.info("{:<20.5f} (Tied) ".format(lpost.model.parameters[i]))
+                log.info("{:<20.5f} (Tied) ".format(lpost.model.parameters[i]))
 
-        self.log.info("\n")
+        log.info("\n")
 
-        self.log.info("Fitting statistics: ")
-        self.log.info(" -- number of data points: %i"%(len(lpost.x)))
+        log.info("Fitting statistics: ")
+        log.info(" -- number of data points: %i"%(len(lpost.x)))
 
         try:
             self.deviance
         except AttributeError:
             self._compute_criteria(lpost)
 
-        self.log.info(" -- Deviance [-2 log L] D = %f.3"%self.deviance)
-        self.log.info(" -- The Akaike Information Criterion of the model is: " +
+        log.info(" -- Deviance [-2 log L] D = %f.3"%self.deviance)
+        log.info(" -- The Akaike Information Criterion of the model is: " +
               str(self.aic) + ".")
 
-        self.log.info(" -- The Bayesian Information Criterion of the model is: " +
+        log.info(" -- The Bayesian Information Criterion of the model is: " +
               str(self.bic) + ".")
 
         try:
@@ -322,13 +319,13 @@ class OptimizationResults(object):
         except AttributeError:
             self._compute_statistics(lpost)
 
-        self.log.info(" -- The figure-of-merit function for this model " +
+        log.info(" -- The figure-of-merit function for this model " +
               " is: %f.5f"%self.merit +
               " and the fit for %i dof is %f.3f"%(self.dof,
                                                   self.merit/self.dof))
 
-        self.log.info(" -- Summed Residuals S = %f.5f"%self.sobs)
-        self.log.info(" -- Expected S ~ %f.5 +/- %f.5"%(self.sexp, self.ssd))
+        log.info(" -- Summed Residuals S = %f.5f"%self.sobs)
+        log.info(" -- Expected S ~ %f.5 +/- %f.5"%(self.sexp, self.ssd))
 
         return
 
@@ -904,9 +901,7 @@ class SamplingResults(object):
     ci_max: float out of [0,100]
         The upper bound percentile for printing credible intervals
         on the parameters
-    log : a logging.getLogger() object, default None
-        You can pass a pre-defined object for logging, else a new 
-        logger will be instantiated
+
 
     Attributes
     ----------
@@ -956,16 +951,7 @@ class SamplingResults(object):
     .. [gelman-rubin] https://projecteuclid.org/euclid.ss/1177011136
     """
 
-    def __init__(self, sampler, ci_min=5, ci_max=95, log=None):
-
-        if log is None:
-            self.log = logging.getLogger('MCMC summary')
-            self.log.setLevel(logging.DEBUG)
-
-            if not self.log.handlers:
-                ch = logging.StreamHandler()
-                ch.setLevel(logging.DEBUG)
-                self.log.addHandler(ch)
+    def __init__(self, sampler, ci_min=5, ci_max=95):
 
         # store all the samples
         self.samples = sampler.get_chain(flat=True)
@@ -1008,7 +994,7 @@ class SamplingResults(object):
         try:
             self.acor = sampler.get_autocorr_time()
         except emcee.autocorr.AutocorrError:
-            self.log.info("Chains too short to compute autocorrelation lengths.")
+            logging.info("Chains too short to compute autocorrelation lengths.")
 
         self.rhat = self._compute_rhat(sampler)
 
@@ -1068,26 +1054,37 @@ class SamplingResults(object):
         self.std = np.std(self.samples, axis=0)
         self.ci = np.percentile(self.samples, [ci_min, ci_max], axis=0)
 
-    def print_results(self):
+    def print_results(self, log=None):
         """
         Print results of the MCMC run on screen or to a log-file.
 
+        Parameters
+        ----------
+        log : a ``logging.getLogger()`` object
+            Object to handle logging output
 
         """
+        if log is None:
+            log = logging.getLogger('MCMC summary')
+            log.setLevel(logging.DEBUG)
+
+            ch = logging.StreamHandler()
+            ch.setLevel(logging.DEBUG)
+            log.addHandler(ch)
 
-        self.log.info("-- The acceptance fraction is: %f.5"%self.acceptance)
+        log.info("-- The acceptance fraction is: %f.5"%self.acceptance)
         try:
-            self.log.info("-- The autocorrelation time is: {}".format(self.acor))
+            log.info("-- The autocorrelation time is: {}".format(self.acor))
         except AttributeError:
             pass
 
-        self.log.info("R_hat for the parameters is: " + str(self.rhat))
+        log.info("R_hat for the parameters is: " + str(self.rhat))
 
-        self.log.info("-- Posterior Summary of Parameters: \n")
-        self.log.info("parameter \t mean \t\t sd \t\t 5% \t\t 95% \n")
-        self.log.info("---------------------------------------------\n")
+        log.info("-- Posterior Summary of Parameters: \n")
+        log.info("parameter \t mean \t\t sd \t\t 5% \t\t 95% \n")
+        log.info("---------------------------------------------\n")
         for i in range(self.ndim):
-            self.log.info("theta[" + str(i) + "] \t " +
+            log.info("theta[" + str(i) + "] \t " +
                   str(self.mean[i]) + "\t" + str(self.std[i]) + "\t" +
                   str(self.ci[0, i]) + "\t" + str(self.ci[1, i]) + "\n")
 
@@ -1154,7 +1151,7 @@ class SamplingResults(object):
 
                     if i == j:
                         ntemp, binstemp, patchestemp = \
-                            ax.hist(samples[:, i], 30, density=True,
+                            ax.hist(samples[:, i], 30, normed=True,
                                     histtype='stepfilled')
                         ax.axis([ymin, ymax, 0, np.max(ntemp)*1.2])
 
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index f30102dc..ab804cdc 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -323,6 +323,11 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
     large_data : bool, default False
         Use only for data larger than 10**7 data points!! Uses zarr and dask for computation.
 
+    save_all : bool, default False
+        Save all intermediate PDSs used for the final average. Use with care.
+        This is likely to fill up your RAM on medium-sized datasets, and to
+        slow down the computation when rebinning.
+
     Attributes
     ----------
     norm: {``leahy`` | ``frac`` | ``abs`` | ``none`` }
@@ -356,7 +361,8 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
 
     """
     def __init__(self, data=None, segment_size=None, norm="frac", gti=None,
-                 silent=False, dt=None, lc=None, large_data=False):
+                 silent=False, dt=None, lc=None, large_data=False,
+                 save_all=False):
 
         if lc is not None:
             warnings.warn("The lc keyword is now deprecated. Use data "
@@ -402,6 +408,7 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
 
         self.type = "powerspectrum"
         self.dt = dt
+        self.save_all = save_all
 
         if isinstance(data, EventList):
             lengths = data.gti[:, 1] - data.gti[:, 0]
@@ -414,7 +421,7 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
 
         return
 
-    def _make_segment_spectrum(self, lc, segment_size):
+    def _make_segment_spectrum(self, lc, segment_size, silent=False):
         """
         Split the light curves into segments of size ``segment_size``, and
         calculate a power spectrum for each.
@@ -427,6 +434,11 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
         segment_size : ``numpy.float``
             Size of each light curve segment to use for averaging.
 
+        Other parameters
+        ----------------
+        silent : bool, default False
+            Suppress progress bars
+
         Returns
         -------
         power_all : list of :class:`Powerspectrum` objects
@@ -455,7 +467,7 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
         nphots_all = []
 
         local_show_progress = show_progress
-        if not self.show_progress:
+        if not self.show_progress or silent:
             local_show_progress = lambda a: a
 
         for start_ind, end_ind in \
diff --git a/stingray/utils.py b/stingray/utils.py
index cfd18a9d..7ba63fe1 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -175,7 +175,7 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
     >>> y = np.ones(x.size)
     >>> yerr = np.ones(x.size)
     >>> xbin, ybin, ybinerr, step_size = rebin_data(
-    ...     x, y, 4, yerr=yerr, method='sum')
+    ...     x, y, 4, yerr=yerr, method='sum', dx=0.01)
     >>> np.allclose(ybin, 400)
     True
     >>> np.allclose(ybinerr, 20)
@@ -189,7 +189,10 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
     """
 
     y = np.asarray(y)
-    yerr = np.asarray(apply_function_if_none(yerr, y, np.zeros_like))
+    if yerr is None:
+        yerr = np.zeros_like(y)
+    else:
+        yerr = np.asarray(yerr)
 
     if not dx:
         dx_old = np.diff(x)
@@ -214,16 +217,17 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
     outputerr = np.zeros(xbin.shape[0] - 1, dtype=type(y[0]))
     step_size = np.zeros(xbin.shape[0] - 1)
 
-    for i in range(len(xbin)-1):
-
-        xmin = xbin[i]
-        xmax = xbin[i+1]
-        min_ind = xedges.searchsorted(xmin)
-        max_ind = xedges.searchsorted(xmax)
-
-        output[i] = np.sum(y[min_ind:max_ind-1])
-        outputerr[i] = np.sum(yerr[min_ind:max_ind-1])
-        step_size[i] = len(y[min_ind:max_ind-1])
+    all_x = np.searchsorted(xedges, xbin)
+    min_inds = all_x[:-1]
+    max_inds = all_x[1:]
+    xmins = xbin[:-1]
+    xmaxs = xbin[1:]
+    for i, (xmin, xmax, min_ind, max_ind) in enumerate(zip(xmins, xmaxs, min_inds, max_inds)):
+        filtered_y = y[min_ind:max_ind-1]
+        filtered_yerr = yerr[min_ind:max_ind-1]
+        output[i] = np.sum(filtered_y)
+        outputerr[i] = np.sum(filtered_yerr)
+        step_size[i] = max_ind -1 - min_ind
 
         prev_dx = xedges[min_ind] - xedges[min_ind-1]
         prev_frac = (xedges[min_ind] - xmin)/prev_dx
@@ -231,7 +235,7 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
         outputerr[i] += yerr[min_ind-1]*prev_frac
         step_size[i] += prev_frac
 
-        if not max_ind == len(xedges):
+        if not max_ind == xedges.size:
             dx_post = xedges[max_ind] - xedges[max_ind-1]
             post_frac = (xmax-xedges[max_ind-1])/dx_post
             output[i] += y[max_ind-1]*post_frac


Test patch:
diff --git a/stingray/modeling/tests/test_parameterestimation.py b/stingray/modeling/tests/test_parameterestimation.py
index 70fe1f3c..752c474e 100644
--- a/stingray/modeling/tests/test_parameterestimation.py
+++ b/stingray/modeling/tests/test_parameterestimation.py
@@ -43,15 +43,7 @@ class LogLikelihoodDummy(LogLikelihood):
 
 class OptimizationResultsSubclassDummy(OptimizationResults):
 
-    def __init__(self, lpost, res, neg, log=None):
-        if log is None:
-            self.log = logging.getLogger('Fitting summary')
-            self.log.setLevel(logging.DEBUG)
-            if not self.log.handlers:
-                ch = logging.StreamHandler()
-                ch.setLevel(logging.DEBUG)
-                self.log.addHandler(ch)
-
+    def __init__(self, lpost, res, neg):
         self.neg = neg
         if res is not None:
             self.result = res.fun
@@ -184,18 +176,6 @@ class TestParameterEstimation(object):
         assert sample_res.acceptance > 0.25
         assert isinstance(sample_res, SamplingResults)
 
-# TODO: Fix pooling with the current setup of logprior
-#    @pytest.mark.skipif("not can_sample", "not can_plot")
-#    def test_sampler_pooling(self):
-#        pe = ParameterEstimation()
-#        if os.path.exists("test_corner.pdf"):
-#            os.unlink("test_corner.pdf")
-#        with catch_warnings(RuntimeWarning):
-#            sample_res = pe.sample(self.lpost, [2.0], nwalkers=50, niter=10,
-#                                   burnin=50, print_results=True, plot=True, 
-#                                   pool=True)
-
-
     @pytest.mark.skipif("can_sample")
     def test_sample_raises_error_without_emcee(self):
         pe = ParameterEstimation()
@@ -401,16 +381,7 @@ class TestOptimizationResults(object):
 
 if can_sample:
     class SamplingResultsDummy(SamplingResults):
-        def __init__(self, sampler, ci_min=0.05, ci_max=0.95, log=None):
-
-            if log is None:
-                self.log = logging.getLogger('Fitting summary')
-                self.log.setLevel(logging.DEBUG)
-                if not self.log.handlers:
-                    ch = logging.StreamHandler()
-                    ch.setLevel(logging.DEBUG)
-                    self.log.addHandler(ch)
-
+        def __init__(self, sampler, ci_min=0.05, ci_max=0.95):
             # store all the samples
             self.samples = sampler.get_chain(flat=True)
 
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index 7186f784..258da4ca 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -166,14 +166,10 @@ class TestAveragedCrossspectrumEvents(object):
                                        segment_size=np.inf, dt=self.dt)
 
     def test_coherence(self):
-        with warnings.catch_warnings(record=True) as w:
-            warnings.simplefilter("always")
-
+        with pytest.warns(UserWarning) as w:
             coh = self.acs.coherence()
-
-            assert len(coh[0]) == 4999
-            assert len(coh[1]) == 4999
-            assert issubclass(w[-1].category, UserWarning)
+        assert len(coh[0]) == 4999
+        assert len(coh[1]) == 4999
 
     def test_failure_when_normalization_not_recognized(self):
         with pytest.raises(ValueError):
@@ -623,6 +619,11 @@ class TestAveragedCrossspectrum(object):
         with pytest.warns(UserWarning) as record:
             self.cs = AveragedCrossspectrum(self.lc1, self.lc2, segment_size=1)
 
+    def test_save_all(self):
+        cs = AveragedCrossspectrum(self.lc1, self.lc2, segment_size=1,
+                                   save_all=True)
+        assert hasattr(cs, 'cs_all')
+
     def test_lc_keyword_deprecation(self):
         cs1 = AveragedCrossspectrum(data1=self.lc1, data2=self.lc2,
                                     segment_size=1)
@@ -728,8 +729,12 @@ class TestAveragedCrossspectrum(object):
         lc1 = Lightcurve(x, y1)
         lc2 = Lightcurve(x, y2)
 
-        acs = AveragedCrossspectrum(lc1, lc2, segment_size=5.0, norm="leahy")
+        with pytest.warns(UserWarning) as record:
+            acs = AveragedCrossspectrum(
+                lc1, lc2, segment_size=5.0, norm="leahy")
         assert acs.m == 1
+        assert np.any(["No counts in interval" in r.message.args[0]
+                       for r in record])
 
 
     def test_rebin_with_invalid_type_attribute(self):
@@ -865,14 +870,15 @@ class TestAveragedCrossspectrum(object):
         test_lc1 = simulator.simulate(2)
         test_lc1.counts -= np.min(test_lc1.counts)
 
-        test_lc1 = Lightcurve(test_lc1.time,
-                              test_lc1.counts,
-                              err_dist=test_lc1.err_dist,
-                              dt=dt)
-        test_lc2 = Lightcurve(test_lc1.time,
-                              np.array(np.roll(test_lc1.counts, 2)),
-                              err_dist=test_lc1.err_dist,
-                              dt=dt)
+        with pytest.warns(UserWarning):
+            test_lc1 = Lightcurve(test_lc1.time,
+                                  test_lc1.counts,
+                                  err_dist=test_lc1.err_dist,
+                                  dt=dt)
+            test_lc2 = Lightcurve(test_lc1.time,
+                                  np.array(np.roll(test_lc1.counts, 2)),
+                                  err_dist=test_lc1.err_dist,
+                                  dt=dt)
 
         with warnings.catch_warnings(record=True) as w:
             cs = AveragedCrossspectrum(test_lc1, test_lc2, segment_size=5,
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index 932a9d1f..f3722396 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -280,6 +280,7 @@ class TestLightcurve(object):
                     "the moment")
 
         with warnings.catch_warnings(record=True) as w:
+            warnings.filterwarnings("always")
             lc = Lightcurve(times, counts, err_dist='gauss')
             assert np.any([warn_str in str(wi.message) for wi in w])
 
@@ -558,7 +559,7 @@ class TestLightcurve(object):
         lc1 = Lightcurve(self.times, self.counts)
         with warnings.catch_warnings():
             warnings.simplefilter("ignore", category=UserWarning)
- 
+
             lc2 = Lightcurve(self.times, self.counts, err=self.counts / 2,
                              err_dist="gauss")
         with warnings.catch_warnings(record=True) as w:
@@ -651,7 +652,7 @@ class TestLightcurve(object):
             newlc = lc1.join(lc2)
             # The join operation *averages* the overlapping arrays
             assert np.allclose(newlc.counts, lc1.counts)
-            assert np.any(["MJDref is different in the two light curves" 
+            assert np.any(["MJDref is different in the two light curves"
                            in str(wi.message) for wi in w])
             assert np.any(["The two light curves have overlapping time ranges"
                            in str(wi.message) for wi in w])
@@ -807,7 +808,7 @@ class TestLightcurve(object):
             warnings.simplefilter("ignore", category=UserWarning)
             lc_test = Lightcurve(test_time, test_counts)
         slc = lc_test.split(1.5)
- 
+
         assert len(slc) == 3
 
     def test_threeway_split_has_correct_data_points(self):