Repository: StingraySoftware/stingray
Pull request #270 (2017-10-18 21:49:28)

New version of covariance spectrum class
I've written a new implementation of the `Covariancespectrum` class in response to Issue #262 to add the capability to include light curves instead of event files. I tried with the original implementation, but the changes touched so many parts that it was easier to just write it from scratch. On the plus side, the code has fewer duplications.

The tests pass, but I get a lot of warnings about the excess variance being <0, so I'd appreciate someone with more knowledge of covariance spectra to take a look and hunt for bugs in the implementation (especially the parts that actually calculate the covariance etc).


Patch:
diff --git a/stingray/covariancespectrum.py b/stingray/covariancespectrum.py
index 6e551df0..26e31c9a 100644
--- a/stingray/covariancespectrum.py
+++ b/stingray/covariancespectrum.py
@@ -1,9 +1,11 @@
 # -*- coding: utf-8 -*-
 from __future__ import division
 import collections
+
 import numpy as np
 
 from stingray import Lightcurve
+from stingray.events import EventList
 import stingray.utils as utils
 
 __all__ = ['Covariancespectrum', 'AveragedCovariancespectrum']
@@ -11,27 +13,59 @@ __all__ = ['Covariancespectrum', 'AveragedCovariancespectrum']
 
 class Covariancespectrum(object):
 
-    def __init__(self, event_list, dt, band_interest=None,
+    def __init__(self, data, dt=None, band_interest=None,
                  ref_band_interest=None, std=None):
+
         """
+        Compute a covariance spectrum for the data. The input data can be
+        either in event data or pre-made light curves. Event data can either
+        be in the form of a numpy.ndarray with (time stamp, energy) pairs or
+        a `stingray.events.EventList` object. If light curves are formed ahead
+        of time, then a list of `Lightcurve` objects should be passed to the
+        object, ideally one light curve for each band of interest.
+
+        For the case where the data is input as a list of `Lightcurve` objects,
+        the reference band(s) should either be (1) a single `Lightcurve` object,
+        (2) a list of `Lightcurve` objects with the reference band for each band
+        of interest pre-made, or (3) `None`, in which case reference bands will
+        formed by combining all light curves *except* for the band of interest.
+
+        In the case of event data, `band_interest` and `ref_band_interest` can
+        be (multiple) pairs of energies, and the light curves for the bands of
+        interest and reference bands will be produced dynamically.
+
+
         Parameters
         ----------
-        event_list : numpy 2D array
-            A numpy 2D array with first column as time of arrival and second
-            column as photon energies associated.
+        data : {numpy.ndarray | EventList object | list of Lightcurve objects}
+            `data` contains the time series data, either in the form of a
+            2-D array of (time stamp, energy) pairs for event data, or as a
+            list of light curves.
             Note : The event list must be in sorted order with respect to the
             times of arrivals.
 
         dt : float
             The time resolution of the Lightcurve formed from the energy bin.
-
-        band_interest : iterable of tuples, default All
-            An iterable of tuples with minimum and maximum values of the range
-            in the band of interest. e.g list of tuples, tuple of tuples.
-
-        ref_band_interest : tuple of reference band range, default All
-            A tuple with minimum and maximum values of the range in the band
-            of interest in reference channel.
+            Only used if `data` is an event list.
+
+        band_interest : {None, iterable of tuples}
+            If None, all possible energy values will be assumed to be of
+            interest, and a covariance spectrum in the highest resolution
+            will be produced.
+            Note;
+            If the input is a list of Lightcurve objects, then the user may
+            supply their energy values here, for construction of a
+            reference band.
+
+        ref_band_interest : {None, tuple, Lightcurve, list of Lightcurves}
+            Defines the reference band to be used for comparison with the
+            bands of interest. If None, all bands *except* the band of
+            interest will be used for each band of interest, respectively.
+            Alternatively, a tuple can be given for event list data, which will
+            extract the reference band (always excluding the band of interest),
+            or one may put in a single Lightcurve object to be used (the same
+            for each band of interest) or a list of Lightcurve objects, one for
+            each band of interest.
 
         std : float or np.array or list of numbers
             The term std is used to calculate the excess variance of a band.
@@ -41,18 +75,8 @@ class Covariancespectrum(object):
             is also used as the std. And if the std is an iterable of
             numbers, their mean is used for the same purpose.
 
-
         Attributes
         ----------
-        energy_events : dictionary
-            A dictionary with energy bins as keys and time of arrivals of
-            photons with the same energy as value.
-
-        energy_covar : dictionary
-            A dictionary with mid point of band_interest and their covariance
-            computed with their individual reference band. The covariance
-            values are normalized.
-
         unnorm_covar : np.ndarray
             An array of arrays with mid point band_interest and their
             covariance. It is the array-form of the dictionary `energy_covar`.
@@ -64,18 +88,6 @@ class Covariancespectrum(object):
         covar_error : np.ndarray
             Errors of the normalized covariance spectrum.
 
-        min_time : int
-            Time of arrival of the earliest photon.
-
-        max_time : int
-            Time of arrival of the last photon.
-
-        min_energy : float
-            Energy of the photon with the minimum energy.
-
-        max_energy : float
-            Energy of the photon with the maximum energy.
-
         References
         ----------
         [1] Wilkinson, T. and Uttley, P. (2009), Accretion disc variability
@@ -87,202 +99,255 @@ class Covariancespectrum(object):
         --------
         See https://github.com/StingraySoftware/notebooks repository for
         detailed notebooks on the code.
-        """
-
-        # This parameter is used to identify whether the current object is
-        # an instance of Covariancespectrum or AveragedCovariancespectrum.
-        self.avg_covar = False
 
-        self._init_vars(event_list, dt, band_interest,
-                        ref_band_interest, std)
-
-        # A dictionary with energy bin as key and events as value of the key
-        self.energy_events = {}
+        """
 
-        self._construct_energy_events(self.energy_events)
+        self.dt = dt
+        self.std = std
 
-        self._update_energy_events(self.energy_events)
+        # check whether data is an EventList object:
+        if isinstance(data, EventList):
+            data = np.vstack([data.time, data.energy]).T
 
-        # The dictionary with covariance spectrum for each energy bin
-        self.energy_covar = {}
+        # check whether the data contains a list of Lightcurve objects
+        if isinstance(data[0], Lightcurve):
+            self.use_lc = True
+            self.lcs = data
+        else:
+            self.use_lc = False
 
-        self._construct_energy_covar(self.energy_events, self.energy_covar)
+        # if band_interest is None, extract the energy bins and make an array
+        # with the lower and upper bounds of the energy bins
+        if not band_interest:
+            if not self.use_lc:
+                self._create_band_interest(data)
+            else:
+                self.band_interest = np.vstack([np.arange(len(data)),
+                                                np.arange(1, len(data)+1, 1)]).T
+        else:
+            if np.size(band_interest) < 2:
+                raise ValueError('band_interest must contain at least 2 values '
+                                 '(minimum and maximum values for each band) '
+                                 'and be a 2D array!')
+
+            self.band_interest = np.atleast_2d(band_interest)
+
+        if self.use_lc is False and not dt:
+            raise ValueError("If the input data is event data, the dt keyword "
+                             "must be set and supply a time resolution for "
+                             "creating light curves!")
+
+        # if we don't have light curves already, make them:
+        if not self.use_lc:
+            if not np.all(np.diff(data, axis=0).T[0] >= 0):
+                utils.simon("The event list must be sorted with respect to "
+                            "times of arrivals.")
+                data = data[data[:, 0].argsort()]
+
+            self.lcs = self._make_lightcurves(data)
+
+        # check whether band of interest contains a Lightcurve object:
+        if np.size(ref_band_interest) == 1  or isinstance(ref_band_interest,
+                                                          Lightcurve):
+            if isinstance(ref_band_interest, Lightcurve):
+                self.ref_band_lcs = ref_band_interest
+            # ref_band_interest must either be a Lightcurve, or must have
+            # multiple entries
+
+            elif ref_band_interest is None:
+                if self.use_lc:
+                    self.ref_band_lcs = \
+                        self._make_reference_bands_from_lightcurves(ref_band_interest)
+                else:
+                    self.ref_band_lcs = \
+                        self._make_reference_bands_from_event_data(data)
+            else:
+                print("This is where I should be!")
+                raise ValueError("ref_band_interest must contain either "
+                                 "a Lightcurve object, a list of Lightcurve "
+                                 "objects or a tuple of length 2.")
+        else:
+            print("I have multiple ref_band_interest!")
+            # check whether ref_band_interest is a list of light curves
+            if isinstance(ref_band_interest[0], Lightcurve):
+                self.ref_band_lcs = ref_band_interest
+                assert len(ref_band_interest) == len(self.lcs), "The list of " \
+                                                                "reference light " \
+                                                                "curves must have " \
+                                                                "the same length as " \
+                                                                "the list of light curves" \
+                                                                "of interest."
+            # if not, it must be a tuple, so we're going to make a list of light
+            # curves
+            else:
+                if self.use_lc:
+                    print("I am in use_lc")
+                    print("bounds in intro: " + str(ref_band_interest))
+                    self.ref_band_lcs = \
+                        self._make_reference_bands_from_lightcurves(bounds=
+                                                                    ref_band_interest)
+                else:
+                    print("I am in the wrong place!")
+                    self.ref_band_lcs = \
+                        self._make_reference_bands_from_event_data(data)
+
+        self._construct_covar()
+
+    def _make_reference_bands_from_event_data(self, data, bounds=None):
+
+        if not bounds:
+            bounds = [np.min(data[:, 1]), np.max(data[:, 1])]
+
+        print("Band interest: " + str(self.band_interest))
+        print("Bounds: " + str(bounds))
+        if bounds[1] <= np.min(self.band_interest[:, 0]) or \
+           bounds[0] >= np.max(self.band_interest[:, 1]):
+            elow = bounds[0]
+            ehigh = bounds[1]
+
+            toa = data[np.logical_and(
+                data[:, 1] >= elow,
+                data[:, 1] <= ehigh)]
+
+            lc_all = Lightcurve.make_lightcurve(toa, self.dt,
+                                                tstart=self.tstart,
+                                                tseg=self.tseg)
 
-    def _init_vars(self, event_list, dt, band_interest,
-                   ref_band_interest, std):
-        """
-        Check for consistency with input variables and declare public ones.
-        """
-        if not np.all(np.diff(event_list, axis=0).T[0] >= 0):
-            utils.simon("The event list must be sorted with respect to "
-                        "times of arrivals.")
-            event_list = event_list[event_list[:, 0].argsort()]
+        else:
 
-        self.event_list = event_list
+            lc_all = []
+            for i, b in enumerate(self.band_interest):
+                elow = b[0]
+                ehigh = b[1]
 
-        self.event_list_T = event_list.T
+                emask1 = data[np.logical_and(
+                    data[:, 1] <= elow,
+                    data[:, 1] >= bounds[0])]
 
-        self._init_special_vars()
+                emask2 = data[np.logical_and(
+                    data[:, 1] <= bounds[1],
+                    data[:, 1] >= ehigh)]
 
-        if ref_band_interest is None:
-            ref_band_interest = (self.min_energy, self.max_energy)
+                toa = np.vstack([emask1, emask2])
+                lc = Lightcurve.make_lightcurve(toa, self.dt,
+                                                tstart=self.tstart,
+                                                tseg=self.tseg)
+                lc_all.append(lc)
 
-        assert type(ref_band_interest) in (list, tuple), "Ref Band interest " \
-                                                         "should be either " \
-                                                         "tuple or list."
+        return lc_all
 
-        assert len(ref_band_interest) == 2, "Band interest should be a tuple" \
-                                            " with min and max energy value " \
-                                            "for the reference band."
-        self.ref_band_interest = ref_band_interest
+    def _make_reference_bands_from_lightcurves(self, bounds=None):
 
-        if band_interest is not None:
-            for element in list(band_interest):
-                assert type(element) in (list, tuple), \
-                    "band_interest should be iterable of either tuple or list."
-                assert len(element) == 2, "Band interest should be a tuple " \
-                                          "with min and max energy values."
 
-        self.band_interest = band_interest
-        self.dt = dt
+        print("bounds: " + str(bounds))
+        if not bounds:
+            bounds_idx = [0, len(self.band_interest)]
 
-        self.std = std
+        else:
+            low_bound = self.band_interest.searchsorted(bounds[0])
+            high_bound = self.band_interest.searchsorted(bounds[1])
 
-    def _init_special_vars(self, T_start=None, T_end=None):
-        """
-        Method to set mininum and maximum time and energy parameters. It has
-        been separated from the main init method due to multiple calls from
-        AveragedCovariancespectrum.
-        """
-        self.min_energy = np.min(self.event_list_T[1][T_start:T_end])
-        self.max_energy = np.max(self.event_list_T[1][T_start:T_end])
-        self.min_time = np.min(self.event_list_T[0][T_start:T_end])
-        self.max_time = np.max(self.event_list_T[0][T_start:T_end])
-
-    def _construct_energy_events(self, energy_events, T_start=None, T_end=None):
-        # The T_start and T_end parameters are for the purpose of
-        # AveragedCovariancespectrum where the range of consideration
-        # is defined.
-        event_list_T = np.array([self.event_list_T[0][T_start: T_end],
-                                 self.event_list_T[1][T_start: T_end]])
-        least_count = np.diff(np.unique(event_list_T[1])).min()
-
-        # An array of unique energy values
-        unique_energy = np.unique(event_list_T[1])
-
-        for i in range(len(unique_energy) - 1):
-            energy_events[unique_energy[i] + least_count*0.5] = []
-
-        # Add time of arrivals to corresponding energy bins
-        # For each bin except the last one, the lower bound is included and
-        # the upper bound is excluded.
-        for energy in energy_events.keys():
-            # The last energy bin
-            if energy == self.max_energy - least_count*0.5:
-                toa = event_list_T[0][np.logical_and(
-                    event_list_T[1] >= energy - least_count*0.5,
-                    event_list_T[1] <= energy + least_count*0.5)]
-                energy_events[energy] = sorted(toa)
-            else:
-                toa = event_list_T[0][np.logical_and(
-                    event_list_T[1] >= energy - least_count*0.5,
-                    event_list_T[1] < energy + least_count*0.5)]
-                energy_events[energy] = sorted(toa)
+            bounds_idx = [low_bound, high_bound]
 
-    def _update_energy_events(self, energy_events):
-        """
-        In case of a specific band interest, merge the required energy bins
-        into one with the new key as the mid-point of the band interest.
-        """
-        if self.band_interest is not None:
-            energy_events_ = {}
-            for band in list(self.band_interest):
-                mid_bin = (band[0] + band[1]) / 2
-                energy_events_[mid_bin] = []
+        lc_all = []
+        for i, b in enumerate(self.band_interest):
 
-                # Modify self.energy_events to form a band with one key
-                for key in list(energy_events.keys()):
-                    if key >= band[0] and key <= band[1]:
-                        energy_events_[mid_bin] += energy_events[key]
-                        del energy_events[key]
+            # initialize empty counts array
+            counts = np.zeros_like(self.lcs[0].counts)
+            for j in range(bounds_idx[0], bounds_idx[1], 1):
+                if i == j:
+                    continue
+                else:
+                    counts += self.lcs[j].counts
 
-            energy_events.update(energy_events_)
+            # make a combined light curve
+            lc = Lightcurve(self.lcs[0].time, counts)
 
-    def _init_energy_covar(self, energy_events, energy_covar):
-        """
-        Initialize the energy_covar dictionary for further computations.
-        """
-        # Initialize it with empty mapping
-        if self.band_interest is None:
-            for key in energy_events.keys():
-                energy_covar[key] = []
-        else:
-            for band in list(self.band_interest):
-                mid_bin = (band[0] + band[1]) / 2
-                energy_covar[mid_bin] = []
+            # add to list of reference light curves
+            lc_all.append(lc)
 
-        if not self.avg_covar:
-            # Error in covariance
-            self.covar_error = {}
+        return lc_all
 
-    def _construct_energy_covar(self, energy_events, energy_covar,
-                                xs_var=None):
-        """Form the actual output covariance dictionary and array."""
-        self._init_energy_covar(energy_events, energy_covar)
+    def _construct_covar(self):
 
-        if not self.avg_covar:
-            xs_var = dict()
+        self.avg_covar = False
+        covar = np.zeros(len(self.lcs))
+        covar_err = np.zeros(len(self.lcs))
+        xs_var = np.zeros(len(self.lcs))
 
-        for energy in energy_covar.keys():
-            lc, lc_ref = self._create_lc_and_lc_ref(energy, energy_events)
+        for i in range(len(self.lcs)):
+            lc = self.lcs[i]
 
-            covar = self._compute_covariance(lc, lc_ref)
+            if np.size(self.ref_band_lcs) == 1 or isinstance(self.ref_band_lcs,
+                                                             Lightcurve):
+                lc_ref = self.ref_band_lcs
+            else:
+                lc_ref = self.ref_band_lcs[i]
 
-            energy_covar[energy] = covar
-            if not self.avg_covar:
-                self.covar_error[energy] = self._calculate_covariance_error(
-                                                lc, lc_ref)
+            cv = self._compute_covariance(lc, lc_ref)
+            cv_err = self._calculate_covariance_error(lc, lc_ref)
 
-            # Excess variance in ref band
-            xs_var[energy] = self._calculate_excess_variance(lc_ref)
+            covar[i] = cv
+            covar_err[i] = cv_err
 
-        for key, value in energy_covar.items():
-            if not xs_var[key] > 0:
+            xs = self._calculate_excess_variance(lc_ref)
+            if not xs > 0:
                 utils.simon("The excess variance in the reference band is "
                             "negative. This implies that the reference "
                             "band was badly chosen. Beware that the "
                             "covariance spectra will have NaNs!")
 
-        if not self.avg_covar:
-            self.unnorm_covar = np.vstack(energy_covar.items())
-            energy_covar[key] = value / (xs_var[key])**0.5
+            xs_var[i] = xs
+
+        self.unnorm_covar = covar
+        energy_covar = covar / xs_var**0.5
+
+        self.covar = energy_covar
+
+        self.covar_error = covar_err
+
+        return
+
+    def _make_lightcurves(self, data):
+
+        self.tstart = np.min(data[:, 0])
+        self.tend = np.max(data[:, 0])
+
+        self.tseg = self.tend - self.tstart
+
+        lc_all = []
+
+        for i, b in enumerate(self.band_interest):
+            elow = b[0]
+            ehigh = b[1]
+
+            toa = data[np.logical_and(
+                data[:, 1] >= elow,
+                data[:, 1] <= ehigh)]
+
+            lc = Lightcurve.make_lightcurve(toa, self.dt, tstart=self.tstart,
+                                            tseg=self.tseg)
+            lc_all.append(lc)
 
-            self.covar = np.vstack(energy_covar.items())
+        return lc_all
 
-            self.covar_error = np.vstack(self.covar_error.items())
+    def _create_band_interest(self, data):
 
-    def _create_lc_and_lc_ref(self, energy, energy_events):
-        lc = Lightcurve.make_lightcurve(
-                energy_events[energy], self.dt, tstart=self.min_time,
-                tseg=self.max_time - self.min_time)
+        unique_energy = np.unique(data[:, 1])
+        energ_diff = np.diff(unique_energy)
 
-        # Calculating timestamps for lc_ref
-        toa_ref = []
-        for key, value in energy_events.items():
-            if key >= self.ref_band_interest[0] and \
-                    key <= self.ref_band_interest[1]:
-                if key != energy:
-                    toa_ref.extend(value)
+        energy_low = np.zeros_like(unique_energy)
+        energy_high = np.zeros_like(unique_energy)
 
-        toa_ref = np.array(sorted(toa_ref))
+        energy_low[:-1] = unique_energy[:-1] - 0.5 * energ_diff
+        energy_high[:-1] = unique_energy[:-1] + 0.5 * energ_diff
 
-        lc_ref = Lightcurve.make_lightcurve(
-                toa_ref, self.dt, tstart=self.min_time,
-                tseg=self.max_time - self.min_time)
+        energy_low[-1] = unique_energy[-1] - 0.5 * energ_diff[-1]
+        energy_high[-1] = unique_energy[-1] + 0.5 * energ_diff[-1]
 
-        assert len(lc.time) == len(lc_ref.time)
+        energy_list = np.vstack([energy_low, energy_high]).T
 
-        return lc, lc_ref
+        self.band_interest = energy_list
 
     def _calculate_excess_variance(self, lc):
         """Calculate excess variance in a band with the standard deviation."""
@@ -315,52 +380,63 @@ class Covariancespectrum(object):
         # Standard deviation of light curve
         err_x = self._calculate_std(lc_x)
         # Number of time bins in lightcurve
-        N = lc_x.n
+        nn = lc_x.n
         # Number of segments averaged
         if not self.avg_covar:
-            M = 1
+            mm = 1
         else:
-            M = self.nbins
+            mm = self.nbins
 
         num = xs_x*err_y + xs_y*err_x + err_x*err_y
-        denom = N * M * xs_y
+        denom = nn * mm * xs_y
 
         return (num / denom)**0.5
 
 
 class AveragedCovariancespectrum(Covariancespectrum):
-    def __init__(self, event_list, dt, segment_size, band_interest=None,
+
+    def __init__(self, data, segment_size, dt=None, band_interest=None,
                  ref_band_interest=None, std=None):
+
         """
-        Make an averaged covariance spectrum by segmenting the light curve
-        formed, calculating covariance for each segment and then averaging
-        the resulting covariance spectra.
+        Compute a covariance spectrum for the data.
 
         Parameters
         ----------
-        event_list : numpy 2D array
-            A numpy 2D array with first column as time of arrival and second
-            column as photon energies associated.
+        data : {numpy.ndarray | list of Lightcurve objects}
+            `data` contains the time series data, either in the form of a
+            2-D array of (time stamp, energy) pairs for event data, or as a
+            list of light curves.
             Note : The event list must be in sorted order with respect to the
             times of arrivals.
 
-        dt : float
-            The time resolution of the Lightcurve formed from the energy bin.
-
         segment_size : float
-            The size of each segment to average. Note that if the total
-            duration of each Lightcurve object formed is not an integer
-            multiple of the segment_size, then any fraction left-over at the
-            end of the time series will be lost.
-
+            The length of each segment in the averaged covariance spectrum.
+            The number of segments will be calculated automatically using the
+            total length of the data set and the segment_size defined here.
 
-        band_interest : iterable of tuples, default All
-            An iterable of tuples with minimum and maximum values of the range
-            in the band of interest. e.g list of tuples, tuple of tuples.
-
-        ref_band_interest : tuple of reference band range, default All
-            A tuple with minimum and maximum values of the range in the band
-            of interest in reference channel.
+        dt : float
+            The time resolution of the Lightcurve formed from the energy bin.
+            Only used if `data` is an event list.
+
+        band_interest : {None, iterable of tuples}
+            If None, all possible energy values will be assumed to be of
+            interest, and a covariance spectrum in the highest resolution
+            will be produced.
+            Note;
+            If the input is a list of Lightcurve objects, then the user may
+            supply their energy values here, for construction of a
+            reference band.
+
+        ref_band_interest : {None, tuple, Lightcurve, list of Lightcurve}
+            Defines the reference band to be used for comparison with the
+            bands of interest. If None, all bands *except* the band of
+            interest will be used for each band of interest, respectively.
+            Alternatively, a tuple can be given for event list data, which will
+            extract the reference band (always excluding the band of interest),
+            or one may put in a single Lightcurve object to be used (the same
+            for each band of interest) or a list of Lightcurve objects, one for
+            each band of interest.
 
         std : float or np.array or list of numbers
             The term std is used to calculate the excess variance of a band.
@@ -370,18 +446,8 @@ class AveragedCovariancespectrum(Covariancespectrum):
             is also used as the std. And if the std is an iterable of
             numbers, their mean is used for the same purpose.
 
-
         Attributes
         ----------
-        energy_events : dictionary
-            A dictionary with energy bins as keys and time of arrivals of
-            photons with the same energy as value.
-
-        energy_covar : dictionary
-            A dictionary with mid point of band_interest and their covariance
-            computed with their individual reference band. The covariance
-            values are normalized.
-
         unnorm_covar : np.ndarray
             An array of arrays with mid point band_interest and their
             covariance. It is the array-form of the dictionary `energy_covar`.
@@ -393,108 +459,73 @@ class AveragedCovariancespectrum(Covariancespectrum):
         covar_error : np.ndarray
             Errors of the normalized covariance spectrum.
 
-        min_time : int
-            Time of arrival of the earliest photon.
-
-        max_time : int
-            Time of arrival of the last photon.
-
-        min_energy : float
-            Energy of the photon with the minimum energy.
-
-        max_energy : float
-            Energy of the photon with the maximum energy.
+        References
+        ----------
+        http://arxiv.org/pdf/1405.6575v2.pdf Equation 15
 
         """
-        # Set parameter to distinguish between parent class and derived class.
-        self.avg_covar = True
 
-        self._init_vars(event_list, dt, band_interest, ref_band_interest, std)
         self.segment_size = segment_size
 
-        self._make_averaged_covar_spectrum()
-
-        self._init_covar_error()
+        Covariancespectrum.__init__(self, data, dt=dt,
+                                    band_interest=band_interest,
+                                    ref_band_interest=ref_band_interest,
+                                    std=std)
 
-        self._calculate_covariance_error()
-
-    def _make_averaged_covar_spectrum(self):
-        """
-        Calls methods from base class for every segment and calculates averaged
-        covariance and error.
-        """
-        self.nbins = int((self.max_time - self.min_time + 1) / self.segment_size)
+    def _construct_covar(self):
 
-        for n in range(self.nbins):
-            tstart = self.min_time + n*self.segment_size
-            tend = self.min_time + self.segment_size*(n+1) - 1
-            indices = np.intersect1d(np.where(self.event_list_T[0] >= tstart),
-                                     np.where(self.event_list_T[0] <= tend))
-
-            # Set minimum and maximum values for the specified indices value
-            self._init_special_vars(T_start=indices[0], T_end=indices[-1]+1)
-
-            energy_events = {}
-
-            self._construct_energy_events(energy_events, T_start=indices[0],
-                                          T_end=indices[-1]+1)
+        self.avg_covar = True
 
-            self._update_energy_events(energy_events)
+        start_time = self.lcs[0].time[0]
 
-            energy_covar = {}
-            xs_var = {}
-            self._construct_energy_covar(energy_events, energy_covar,
-                                         xs_var)
+        covar = np.zeros(len(self.lcs))
+        covar_err = np.zeros(len(self.lcs))
+        xs_var = np.zeros(len(self.lcs))
 
-            if n == 0:  # Declare
-                self.energy_covar = energy_covar
-                self.xs_var = xs_var
-            else:  # Sum up
-                for key in energy_covar.keys():
-                    self.energy_covar[key] = self.energy_covar.get(key, 0) + \
-                                             energy_covar[key]
-                    self.xs_var[key] = self.xs_var.get(key, 0) + xs_var[key]
+        for i in range(len(self.lcs)):
+            lc = self.lcs[i]
 
-        # Now divide with total number of bins for averaging
-        for key in self.energy_covar.keys():
-            self.energy_covar[key] /= self.nbins
-            self.xs_var[key] /= self.nbins
+            if np.size(self.ref_band_lcs) == 1:
+                lc_ref = self.ref_band_lcs
+            else:
+                lc_ref = self.ref_band_lcs[i]
 
-        self.unnorm_covar = np.vstack(self.energy_covar.items())
+            tstart = start_time
+            tend = start_time + self.segment_size
+            cv = 0.0
+            cv_err = 0.0
+            xs = 0.0
 
-        for key, value in self.energy_covar.items():
-            self.energy_covar[key] = value / (self.xs_var[key])**0.5
+            self.nbins = int((tend - tstart)/self.segment_size)
+            for k in range(self.nbins):
+                start_ind = lc.time.searchsorted(tstart)
+                end_ind = lc.time.searchsorted(tend)
 
-        self.covar = np.vstack(self.energy_covar.items())
+                lc_seg = lc.truncate(start=start_ind, stop=end_ind)
+                lc_ref_seg = lc_ref.truncate(start=start_ind, stop=end_ind)
 
-    def _init_covar_error(self):
-        """Initialize dictionaries separately for the calculation of error."""
-        self.energy_events = {}
-        self._construct_energy_events(self.energy_events)
-        self._update_energy_events(self.energy_events)
-        self.covar_error = {}
-        self._init_energy_covar(self.energy_events, self.covar_error)
+                cv += self._compute_covariance(lc_seg, lc_ref_seg)
+                cv_err += self._calculate_covariance_error(lc_seg, lc_ref_seg)
+                xs += self._calculate_excess_variance(lc_ref_seg)
+                if not xs > 0:
+                    utils.simon("The excess variance in the reference band is "
+                                "negative. This implies that the reference "
+                                "band was badly chosen. Beware that the "
+                                "covariance spectra will have NaNs!")
 
-    def _calculate_covariance_error(self):
-        """
-        Calculate Covariance error on the averaged quantities.
+                tstart += self.segment_size
+                tend += self.segment_size
 
-        References
-        ----------
-        http://arxiv.org/pdf/1405.6575v2.pdf Equation 15
-
-        """
-        for energy in self.covar_error.keys():
-            lc, lc_ref = self._create_lc_and_lc_ref(energy, self.energy_events)
 
-            xs_y = self._calculate_excess_variance(lc_ref)
+            covar[i] = cv/self.nbins
+            covar_err[i] = cv_err/self.nbins
+            xs_var[i] = xs/self.nbins
 
-            err_x = self._calculate_std(lc)
-            err_y = self._calculate_std(lc_ref)
+        self.unnorm_covar = covar
+        energy_covar = covar / xs_var**0.5
 
-            covar = self.energy_covar[energy]
+        self.covar = energy_covar
 
-            num = (covar**2)*err_y + xs_y*err_x + err_x*err_y
-            denom = 2*self.nbins*xs_y
+        self.covar_error = covar_err
 
-            self.covar_error[energy] = (num / denom)**0.5
+        return


Test patch:
diff --git a/stingray/tests/test_covariancespectrum.py b/stingray/tests/test_covariancespectrum.py
index 5dd04760..5091f9a2 100644
--- a/stingray/tests/test_covariancespectrum.py
+++ b/stingray/tests/test_covariancespectrum.py
@@ -5,9 +5,9 @@ import pytest
 import warnings
 
 from stingray import AveragedCovariancespectrum, Covariancespectrum, Lightcurve
+from stingray.events import EventList
 
-
-class TestCovariancespectrum(object):
+class TestCovariancespectrumwithEvents(object):
 
     def setup_class(self):
         event_list = np.array([[1, 2], [1, 4], [1, 8], [3, 6], [3, 1], [4, 2],
@@ -17,35 +17,139 @@ class TestCovariancespectrum(object):
                                [9, 2], [9, 5], [9, 9]])
         self.event_list = event_list
 
+
+    def test_class_fails_if_events_is_set_but_dt_is_not(self):
+        with pytest.raises(ValueError):
+            c = Covariancespectrum(self.event_list)
+
+    def test_set_lc_flag_correctly(self):
+        c = Covariancespectrum(self.event_list, dt=1)
+        assert c.use_lc is False
+
     def test_covar_without_any_band_interest(self):
-        c = Covariancespectrum(self.event_list, 1)
-        assert len(c.covar) == 8
+        c = Covariancespectrum(self.event_list, dt=1)
+        assert len(c.covar) == 9
+
+    def test_creates_band_interest_correctly(self):
+        c = Covariancespectrum(self.event_list, dt=1)
+        energies = np.arange(1, 10, 1)
+        band_low = np.zeros_like(energies)
+        band_high = np.zeros_like(energies)
+        ediff = np.diff(energies)
+
+        band_low[:-1] = energies[:-1] - 0.5*ediff
+        band_high[:-1] = energies[:-1] + 0.5*ediff
+        band_low[-1] = energies[-1] - 0.5*ediff[-1]
+        band_high[-1] = energies[-1] + 0.5*ediff[-1]
+
+        band_interest = np.vstack([band_low, band_high]).T
+
+        assert np.allclose(band_interest, c.band_interest)
 
     def test_init_with_invalid_band_interest(self):
-        with pytest.raises(AssertionError):
-            c = Covariancespectrum(self.event_list, 1, band_interest=(1, 4))
+        with pytest.raises(ValueError):
+            c = Covariancespectrum(self.event_list, dt=1, band_interest=(1))
 
     def test_init_with_invalid_ref_band_interest(self):
-        with pytest.raises(AssertionError):
-            c = Covariancespectrum(self.event_list, 1, ref_band_interest=(0))
+        with pytest.raises(ValueError):
+            c = Covariancespectrum(self.event_list, dt=1, ref_band_interest=(0))
 
     def test_covar_with_both_bands(self):
-        c = Covariancespectrum(self.event_list, 1,
+        c = Covariancespectrum(self.event_list, dt=1,
                                band_interest=[(2, 6), (8, 9)],
                                ref_band_interest=(1, 10))
         assert len(c.covar) == 2
 
+    def test_with_separate_reference_band(self):
+        c = Covariancespectrum(self.event_list, dt=1,
+                               band_interest=[(2, 6)],
+                               ref_band_interest=((7, 10)))
+
     def test_with_unsorted_event_list(self):
         event_list = np.array([[2, 1], [2, 3], [1, 2], [5, 2], [4, 1]])
         with warnings.catch_warnings(record=True) as w:
-            c = Covariancespectrum(event_list, 1)
+            c = Covariancespectrum(event_list, dt=1)
             assert "must be sorted" in str(w[0].message)
 
     def test_with_std_as_iterable(self):
-        c = Covariancespectrum(self.event_list, 1, std=[1, 2])
+        c = Covariancespectrum(self.event_list, dt=1, std=[1, 2])
 
     def test_with_std_as_a_single_float(self):
-        c = Covariancespectrum(self.event_list, 1, std=2.55)
+        c = Covariancespectrum(self.event_list, dt=1, std=2.55)
+
+    def test_with_eventlist_object(self):
+        e = EventList(time = self.event_list[:,0],
+                      energy = self.event_list[:,1])
+
+        c = Covariancespectrum(e, dt=1)
+
+class TestCovariancewithLightcurves(object):
+    def setup_class(self):
+
+        time = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])
+        counts1 = np.array([0, 0, 1, 0, 1, 1, 0, 0, 0])
+        counts2 = np.array([1, 0, 0, 1, 0, 0, 0, 0, 1])
+        counts3 = np.array([0, 0, 0, 0, 1, 1, 1, 0, 0])
+        counts4 = np.array([1, 0, 0, 1, 0, 0, 1, 0, 0])
+        counts5 = np.array([0, 0, 0, 1, 0, 0, 1, 0, 1])
+        counts6 = np.array([0, 0, 1, 1, 1, 0, 0, 0, 0])
+        counts7 = np.array([0, 0, 0, 1, 1, 0, 1, 0, 0])
+        counts8 = np.array([1, 0, 0, 0, 1, 0, 0, 1, 0])
+        counts9 = np.array([0, 0, 0, 1, 1, 0, 0, 0, 1])
+
+        lc1 = Lightcurve(time, counts1)
+        lc2 = Lightcurve(time, counts2)
+        lc3 = Lightcurve(time, counts3)
+        lc4 = Lightcurve(time, counts4)
+        lc5 = Lightcurve(time, counts5)
+        lc6 = Lightcurve(time, counts6)
+        lc7 = Lightcurve(time, counts7)
+        lc8 = Lightcurve(time, counts8)
+        lc9 = Lightcurve(time, counts9)
+
+        self.lcs = [lc1, lc2, lc3, lc4, lc5, lc6, lc7, lc8, lc9]
+
+        self.ref_band_interest = Lightcurve(time,
+                                            np.array([2, 3, 1,
+                                                      4, 4, 0,
+                                                      1, 3, 2]))
+
+    def test_class_initializes_with_lightcurves(self):
+        c = Covariancespectrum(self.lcs)
+
+    def test_set_lc_flag_correctly(self):
+        c = Covariancespectrum(self.lcs)
+        assert c.use_lc is True
+        assert c.lcs == self.lcs
+
+    def test_single_lc_as_ref_band_interest(self):
+        c = Covariancespectrum(self.lcs,
+                               ref_band_interest=self.ref_band_interest)
+
+    def test_lc_with_no_band_interest(self):
+        c = Covariancespectrum(self.lcs,
+                               ref_band_interest=None)
+
+        band_interest = np.vstack([np.arange(len(self.lcs)),
+                        np.arange(1, len(self.lcs) + 1, 1)]).T
+
+        assert np.all(c.band_interest == band_interest)
+
+    def test_failure_without_valid_band_interest(self):
+        with pytest.raises(ValueError):
+            c = Covariancespectrum(self.lcs, band_interest=1.0)
+
+    def test_ref_band_interest_created_correctly(self):
+        c = Covariancespectrum(self.lcs, ref_band_interest=None)
+        assert len(c.ref_band_lcs) == len(self.lcs)
+
+    def test_class_fails_with_wrong_number_of_ref_band_interest_lcs(self):
+        c = Covariancespectrum(self.lcs, ref_band_interest=self.lcs[::-1])
+        assert len(c.lcs) == len(c.ref_band_lcs)
+
+    def test_ref_band_interest_get_saved_correctly(self):
+        with pytest.raises(AssertionError):
+            c = Covariancespectrum(self.lcs, ref_band_interest=self.lcs[:-1])
 
 
 class TestAveragedCovariancespectrum(object):
@@ -60,8 +164,9 @@ class TestAveragedCovariancespectrum(object):
 
     def test_with_full_segment(self):
         c = Covariancespectrum(self.event_list, 1)
-        avg_c = AveragedCovariancespectrum(self.event_list, 1, segment_size=10)
+        avg_c = AveragedCovariancespectrum(self.event_list, segment_size=10,
+                                           dt=1)
         assert np.all(avg_c.unnorm_covar == c.unnorm_covar)
 
     def test_with_two_segments(self):
-        avg_c = AveragedCovariancespectrum(self.event_list, 1, segment_size=5)
+        avg_c = AveragedCovariancespectrum(self.event_list, segment_size=5, dt=1)