Repository: StingraySoftware/stingray (python)

Issue:
Add cross-correlation function
Needed for [DAVE-105](https://timelabdev.com/jira/browse/DAVE-105)

Pull request:
Cross Correlation
Closes #172

Adds cross-correlation function for time series. Any normalization of data approach used in literature in welcomed to discuss, a simple method is to normalize data by subtracting mean and dividing by standard deviation before calculating cross-correlation. Some literature also divides correlation by length of input series.

Patch:
diff --git a/.travis.yml b/.travis.yml
index f8847b62..923b5e39 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -2,6 +2,7 @@ language: python
 
 python:
     - 2.7
+    - 3.5
     - 3.6
 
 # Setting sudo to false opts in to Travis-CI container-based builds.
@@ -19,10 +20,6 @@ addons:
             - texlive-latex-extra
             - dvipng
 
-branches:
-  only:
-    - "master"
-
 env:
     global:
         # The following versions are the 'default' for tests, unless
@@ -39,7 +36,7 @@ env:
         # Make sure that egg_info works without dependencies
         - SETUP_CMD='egg_info'
         # Try all python versions with the latest numpy
-        - SETUP_CMD='test --coverage'
+        - SETUP_CMD='test'
 
 matrix:
     fast_finish: true
@@ -56,15 +53,33 @@ matrix:
         - python: 3.5
           env: PIP_DEPENDENCIES='emcee statsmodels' SETUP_CMD='test --coverage'
 
+        # Do a coverage test in Python 2.
+        - python: 2.7
+          env: SETUP_CMD='test --coverage'
+
+        # Do a coverage test in Python 3.
+        - python: 3.5
+          env: SETUP_CMD='test --coverage'
+
+        # Try Astropy development version
+        - python: 2.7
+          env: ASTROPY_VERSION=development
+
         - python: 3.4
           env: NUMPY_VERSION=1.11
-
-        ## Here allowed failures start
+         # Check for sphinx doc build warnings - we do this first because it
+         # may run for a long time
+        - python: 2.7
+          env: SETUP_CMD='build_sphinx -w'
+        - python: 3.5
+          env: ASTROPY_VERSION=development
         # Try older numpy versions
         - python: 2.7
           env: NUMPY_VERSION=1.8
         - python: 2.7
           env: NUMPY_VERSION=1.7
+        - python: 2.7
+          env: NUMPY_VERSION=1.6
         # Try older scipy versions
         - python: 2.7
           env: SCIPY_VERSION=0.16.0
@@ -72,14 +87,8 @@ matrix:
           env: SCIPY_VERSION=0.14.0
         - python: 2.6
           env: SETUP_CMD='test'
-
-        - python: 2.7
-          env: SETUP_CMD='build_sphinx -w'
-        # Try Astropy development version
-        - python: 2.7
-          env: ASTROPY_VERSION=development
-        - python: 3.5
-          env: ASTROPY_VERSION=development
+        - python: 3.3
+          env: SETUP_CMD='test'
 
     allow_failures:
         # Try older numpy versions
@@ -87,6 +96,8 @@ matrix:
           env: NUMPY_VERSION=1.8
         - python: 2.7
           env: NUMPY_VERSION=1.7
+        - python: 2.7
+          env: NUMPY_VERSION=1.6
         # Try older scipy versions
         - python: 2.7
           env: SCIPY_VERSION=0.16.0
@@ -94,6 +105,8 @@ matrix:
           env: SCIPY_VERSION=0.14.0
         - python: 2.6
           env: SETUP_CMD='test'
+        - python: 3.3
+          env: SETUP_CMD='test'
 
         - python: 2.7
           env: SETUP_CMD='build_sphinx -w'
diff --git a/docs/simulator.rst b/docs/simulator.rst
index 2b1a8942..cb058f8e 100644
--- a/docs/simulator.rst
+++ b/docs/simulator.rst
@@ -1,7 +1,8 @@
 .. _stingray-simulator:
 
+***************************************
 Stingray Simulator (`stingray.simulator`)
-*****************************************
+***************************************
 
 Introduction
 ============
@@ -50,7 +51,7 @@ Simulate Method
 Stingray provides multiple ways to simulate a light curve. However, all these methods follow a common recipe::
 
   >>> sim = simulator.Simulator(N=1024, mean=0.5, dt=0.125)
-  >>> lc = sim.simulate(2)
+  >>> simulator.simulate(...)
 
 Using Power-Law Spectrum
 ------------------------
@@ -161,15 +162,13 @@ The `simulator` class provides the functionality to simulate light curves indepe
   1
   >>> lc = sim.get_channel('3.5 - 4.5')
   >>> sim.delete_channel('3.5 - 4.5')
+  0
 
 Alternatively, assume that we have light curves in the simulated energy channels `3.5 - 4.5`, `4.5 - 5.5` and `5.5 - 6.5`. These channels can be retreived or deleted in single commands. 
 
-  >>> sim.count_channels()
-  0
-  >>> sim.simulate_channel('3.5 - 4.5', 2)
-  >>> sim.simulate_channel('4.5 - 5.5', 2)
-  >>> sim.simulate_channel('5.5 - 6.5', 2)
-  >>> chans = sim.get_channels(['3.5 - 4.5','4.5 - 5.5','5.5 - 6.5'])
+  >>> sim.countchannels()
+  3
+  >>> sim.get_channels(['3.5 - 4.5','4.5 - 5.5','5.5 - 6.5'])
   >>> sim.delete_channels(['3.5 - 4.5','4.5 - 5.5','5.5 - 6.5'])
 
 Reference/API
diff --git a/stingray/__init__.py b/stingray/__init__.py
index f08d0c33..2f5c7d06 100644
--- a/stingray/__init__.py
+++ b/stingray/__init__.py
@@ -15,3 +15,4 @@ if not _ASTROPY_SETUP_:
     from stingray import *
     from stingray.exceptions import *
     from stingray.covariancespectrum import *
+    from stingray.crosscorrelation import *
diff --git a/stingray/crosscorrelation.py b/stingray/crosscorrelation.py
new file mode 100644
index 00000000..a4323407
--- /dev/null
+++ b/stingray/crosscorrelation.py
@@ -0,0 +1,209 @@
+from __future__ import division
+import numpy as np
+from scipy import signal
+
+from stingray import lightcurve
+from stingray.exceptions import StingrayError
+import stingray.utils as utils
+
+
+class CrossCorrelation(object):
+    def __init__(self, lc1=None, lc2=None, mode='same'):
+
+        """
+        Make a cross-correlation from a light curves.
+        You can also make an empty Crosscorrelation object to populate with your
+        own cross-correlation data.
+        Parameters
+        ----------
+        lc1: lightcurve.Lightcurve object, optional, default None
+            The first light curve data for correlation calculations.
+        lc2: lightcurve.Lightcurve object, optional, default None
+            The light curve data for the correlation calculations.
+        mode: {'full', 'valid', 'same'}, optional, default 'same'
+            A string indicating the size of the correlation output.
+            See https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.signal.correlate.html
+            for more details.
+            
+        Attributes
+        ----------
+        lc1: lightcurve.Lightcurve 
+            The first light curve data for correlation calculations.
+        lc2: lightcurve.Lightcurve
+            The light curve data for the correlation calculations.        
+        corr: numpy.ndarray
+             An array of correlation data calculated from two lighcurves
+        time_lags: numpy.ndarray
+             An array of all possible time lags against which each point in corr is calculated 
+        dt: float
+             The time resolution of each lightcurve (used in time_lag calculations)
+        time_shift: float
+             Time lag that gives maximum value of correlation between two lightcurves.
+             There will be maximum correlation between lightcurves if one of the lightcurve is shifted by time_shift.
+        n: int
+             Number of points in self.corr(Length of cross-correlation data) 
+        """
+
+        if isinstance(mode, str) is False:
+            raise TypeError("mode must be a string")
+
+        if mode.lower() not in ["full", "valid", "same"]:
+            raise ValueError("mode must be 'full', 'valid' or 'same'!")
+
+        self.mode = mode.lower()
+        self.lc1 = None
+        self.lc2 = None
+
+        # Populate all attributes by None if user passes no lightcurve data
+        if lc1 is None or lc2 is None:
+            if lc1 is not None or lc2 is not None:
+                raise TypeError("You can't do a cross correlation with just one "
+                                "light curve!")
+            else:
+                # both lc1 and lc2 are None.
+                self.corr = None
+                self.time_shift = None
+                self.time_lags = None
+                self.dt = None
+                self.n = None
+                return
+        else:
+            self._make_corr(lc1, lc2)
+
+    def _make_corr(self, lc1, lc2):
+
+        """
+        Creates Crosscorrelation Object.
+        Parameters
+        ----------
+        lc1: lightcurve.Lightcurve object
+            The first light curve data.
+        lc2: lightcurve.Lightcurve object
+            The second light curve data.
+        Returns
+        ----------
+        None
+        """
+
+        if not isinstance(lc1, lightcurve.Lightcurve):
+            raise TypeError("lc1 must be a lightcurve.Lightcurve object")
+        if not isinstance(lc2, lightcurve.Lightcurve):
+            raise TypeError("lc2 must be a lightcurve.Lightcurve object")
+
+        if not np.isclose(lc1.dt, lc2.dt):
+            raise StingrayError("Light curves do not have "
+                                "same time binning dt.")
+        else:
+            # ignore very small differences in dt neglected by np.isclose()
+            lc1.dt = lc2.dt
+            self.dt = lc1.dt
+
+        # self.lc1 and self.lc2 may get assigned values explicitly in which case there is no need to copy data
+        if self.lc1 is None:
+            self.lc1 = lc1
+        if self.lc2 is None:
+            self.lc2 = lc2
+
+        # Subtract means before passing scipy.signal.correlate into correlation
+        lc1_counts = self.lc1.counts - np.mean(self.lc1.counts)
+        lc2_counts = self.lc2.counts - np.mean(self.lc2.counts)
+
+        # Calculates cross-correlation of two lightcurves
+        self.corr = signal.correlate(lc1_counts, lc2_counts, self.mode)
+        self.n = len(self.corr)
+        self.time_shift, self.time_lags, self.n = self.cal_timeshift(dt=self.dt)
+
+    def cal_timeshift(self, dt=1.0):
+        """
+        Creates Crosscorrelation Object.
+        Parameters
+        ----------
+        dt: float , optional, default 1.0
+            Time resolution of lightcurve, should be passed when object is populated with correlation data
+            and no information about light curve can be extracted. Used to calculate time_lags.
+        Returns
+        ----------
+        self.time_shift: float
+             Value of time lag that gives maximum value of correlation between two lightcurves. 
+        self.time_lags: numpy.ndarray
+             An array of time_lags calculated from correlation data 
+        """
+        if self.dt is None:
+            self.dt = dt
+
+        if self.corr is None:
+            if self.lc1 is None or self.lc2 is None:
+                raise StingrayError('lc1 and lc2 should be provided to calculate correlation and time_shift')
+            else:
+                # This will cover very rare case of assigning self.lc1 and self.lc2 and self.corr = None.
+                # In this case, correlation is calculated using self.lc1 and self.lc2 and using that correlation data,
+                # time_shift is calculated.
+                self._make_corr(self.lc1, self.lc2)
+                return self.time_shift, self.time_lags, self.n
+
+        self.n = len(self.corr)
+        dur = int(self.n / 2)
+        # Correlation against all possible lags, positive as well as negative lags are stored
+        x_lags = np.linspace(-dur, dur, self.n)
+        self.time_lags = x_lags * self.dt
+        # time_shift is the time lag for max. correlation
+        self.time_shift = self.time_lags[np.argmax(self.corr)]
+
+        return self.time_shift, self.time_lags, self.n
+
+    def plot(self, labels=None, axis=None, title=None, marker='-', save=False, filename=None):
+        """
+        Plot the :class:`Crosscorrelation` as function using Matplotlib.
+        Plot the Crosscorrelation object on a graph ``self.time_lags`` on x-axis and
+        ``self.corr`` on y-axis
+        Parameters
+        ----------
+        labels : iterable, default None
+            A list of tuple with xlabel and ylabel as strings.
+        axis : list, tuple, string, default None
+            Parameter to set axis properties of Matplotlib figure. For example
+            it can be a list like ``[xmin, xmax, ymin, ymax]`` or any other
+            acceptable argument for `matplotlib.pyplot.axis()` function.
+        title : str, default None
+            The title of the plot.
+        marker : str, default '-'
+            Line style and color of the plot. Line styles and colors are
+            combined in a single format string, as in ``'bo'`` for blue
+            circles. See `matplotlib.pyplot.plot` for more options.
+        save : boolean, optional (default=False)
+            If True, save the figure with specified filename.
+        filename : str
+            File name of the image to save. Depends on the boolean ``save``.
+        """
+
+        try:
+            import matplotlib.pyplot as plt
+        except ImportError:
+            raise ImportError("Matplotlib required for plot()")
+
+        plt.plot(self.time_lags, self.corr, marker)
+        if labels is not None:
+            try:
+                plt.xlabel(labels[0])
+                plt.ylabel(labels[1])
+            except TypeError:
+                utils.simon("``labels`` must be either a list or tuple with "
+                            "x and y labels.")
+                raise
+            except IndexError:
+                utils.simon("``labels`` must have two labels for x and y "
+                            "axes.")
+                # Not raising here because in case of len(labels)==1, only
+                # x-axis will be labelled.
+
+        if axis is not None:
+            plt.axis(axis)
+
+        if title is not None:
+            plt.title(title)
+
+        if save:
+            if filename is None:
+                plt.savefig('corr.png')
+            else:
+                plt.savefig(filename)
diff --git a/stingray/events.py b/stingray/events.py
index 00b03d67..6c2f6d29 100644
--- a/stingray/events.py
+++ b/stingray/events.py
@@ -19,7 +19,7 @@ __all__ = ['EventList']
 
 
 class EventList(object):
-    def __init__(self, time=None, energy=None, ncounts=None, mjdref=0, dt=0, notes="",
+    def __init__(self, time=None, pha=None, ncounts=None, mjdref=0, dt=0, notes="", 
             gti=None, pi=None):
         """
         Make an event list object from an array of time stamps
@@ -35,7 +35,7 @@ class EventList(object):
             The time resolution of the events. Only relevant when using events
             to produce light curves with similar bin time.
 
-        energy: iterable
+        pha: iterable
             A list of array of photon energy values
 
         mjdref : float
@@ -56,7 +56,7 @@ class EventList(object):
             The array of event arrival times, in seconds from the reference
             MJD (self.mjdref)
 
-        energy: numpy.ndarray
+        pha: numpy.ndarray
             The array of photon energy values
 
         ncounts: int
@@ -77,7 +77,7 @@ class EventList(object):
 
         """
         
-        self.energy = None if energy is None else np.array(energy)
+        self.pha = None if pha is None else np.array(pha)
         self.notes = notes
         self.dt = dt
         self.mjdref = mjdref
@@ -91,9 +91,9 @@ class EventList(object):
         else:
             self.time = None
 
-        if (time is not None) and (energy is not None):
-            if len(time) != len(energy):
-                raise ValueError('Lengths of time and energy must be equal.')
+        if (time is not None) and (pha is not None):
+            if len(time) != len(pha):
+                raise ValueError('Lengths of time and pha must be equal.')
 
     def to_lc(self, dt, tstart=None, tseg=None):
         """
@@ -253,10 +253,10 @@ class EventList(object):
 
         if isinstance(spectrum, list) or isinstance(spectrum, np.ndarray):
             
-            energy = np.array(spectrum)[0]
+            pha = np.array(spectrum)[0]
             fluxes = np.array(spectrum)[1]
 
-            if not isinstance(energy, np.ndarray):
+            if not isinstance(pha, np.ndarray):
                 raise IndexError("Spectrum must be a 2-d array or list")
         
         else:
@@ -272,7 +272,7 @@ class EventList(object):
         R = ra.uniform(0, 1, self.ncounts)
 
         # Assign energies to events corresponding to the random numbers drawn
-        self.energy = np.array([energy[np.argwhere(cum_prob ==
+        self.pha = np.array([pha[np.argwhere(cum_prob == 
             min(cum_prob[(cum_prob - r) > 0]))] for r in R])
 
     def join(self, other):
@@ -331,16 +331,16 @@ class EventList(object):
             ev_new.pi = np.concatenate([self.pi, other.pi])
             ev_new.pi = ev_new.pi[order]
 
-        if (self.energy is None) and (other.energy is None):
-            ev_new.energy = None
-        elif (self.energy is None) or (other.energy is None):
-            self.energy = assign_value_if_none(self.energy, np.zeros_like(self.time))
-            other.energy = assign_value_if_none(other.energy,
+        if (self.pha is None) and (other.pha is None):
+            ev_new.pha = None
+        elif (self.pha is None) or (other.pha is None):
+            self.pha = assign_value_if_none(self.pha, np.zeros_like(self.time))
+            other.pha = assign_value_if_none(other.pha,
                                              np.zeros_like(other.time))
 
-        if (self.energy is not None) and (other.energy is not None):
-            ev_new.energy = np.concatenate([self.energy, other.energy])
-            ev_new.energy = ev_new.energy[order]
+        if (self.pha is not None) and (other.pha is not None):
+            ev_new.pha = np.concatenate([self.pha, other.pha])
+            ev_new.pha = ev_new.pha[order]
 
         if self.gti is None and other.gti is not None and len(self.time) > 0:
             self.gti = \
@@ -383,7 +383,7 @@ class EventList(object):
         -------
         ev: `EventList` object
         """
-        attributes = ['time', 'energy', 'ncounts', 'mjdref', 'dt',
+        attributes = ['time', 'pha', 'ncounts', 'mjdref', 'dt', 
                 'notes', 'gti', 'pi']
         data = read(filename, format_, cols=attributes)
 
@@ -405,7 +405,7 @@ class EventList(object):
                 else:
                     values.append(None)
                     
-            return EventList(time=values[0], energy=values[1], ncounts=values[2],
+            return EventList(time=values[0], pha=values[1], ncounts=values[2], 
                 mjdref=values[3], dt=values[4], notes=values[5], gti=values[6], pi=values[7])
 
         elif format_ == 'pickle':
diff --git a/stingray/modeling/scripts.py b/stingray/modeling/scripts.py
index 8aa5d058..bc66c16e 100644
--- a/stingray/modeling/scripts.py
+++ b/stingray/modeling/scripts.py
@@ -99,11 +99,11 @@ def fit_powerspectrum(ps, model, starting_pars, max_post=False, priors=None,
     Now we have to guess starting parameters. For each Lorentzian, we have
     amplitude, centroid position and fwhm, and this pattern repeats for each
     Lorentzian in the fit. The white noise level is the last parameter.
-    >>> t0 = [80, 1., 1.5, 2.5]
+    >>> t0 = [80, 1.5, 2.5]
 
     Let's also make a model to test:
     >>> model_to_test = models.PowerLaw1D() + models.Const1D()
-    >>> model_to_test.amplitude_1.fixed = True
+    >>> model_to_test.x_0_0.fixed = True
 
     We're ready for doing the fit:
     >>> parest, res = fit_powerspectrum(ps, model_to_test, t0)
diff --git a/stingray/varenergyspectrum.py b/stingray/varenergyspectrum.py
index fd46d101..8b780a84 100644
--- a/stingray/varenergyspectrum.py
+++ b/stingray/varenergyspectrum.py
@@ -22,7 +22,6 @@ def _decode_energy_specification(energy_spec):
     True
     >>> a = _decode_energy_specification((1, 4, 2, 'log'))
     >>> np.all(a == [1, 2, 4])
-    True
     """
     if not isinstance(energy_spec, tuple):
         raise ValueError("Energy specification must be a tuple")
@@ -132,8 +131,8 @@ class VarEnergySpectrum(object):
             energies1 = self.events1.pi
             energies2 = self.events2.pi
         else:
-            energies2 = self.events2.energy
-            energies1 = self.events1.energy
+            energies2 = self.events2.pha
+            energies1 = self.events1.pha
 
         gti = cross_two_gtis(self.events1.gti, self.events2.gti)
 


Test patch:
diff --git a/stingray/tests/test_crosscorrelation.py b/stingray/tests/test_crosscorrelation.py
new file mode 100644
index 00000000..fb98d232
--- /dev/null
+++ b/stingray/tests/test_crosscorrelation.py
@@ -0,0 +1,195 @@
+import numpy as np
+
+import pytest
+import warnings
+import os
+
+from stingray import Lightcurve
+from stingray.crosscorrelation import CrossCorrelation
+from stingray.exceptions import StingrayError
+
+try:
+    import matplotlib.pyplot as plt
+
+    HAS_MPL = True
+except ImportError:
+    HAS_MPL = False
+
+
+class TestCrossCorrelation(object):
+    @classmethod
+    def setup_class(cls):
+        cls.lc1 = Lightcurve([1, 2, 3, 4, 5], [2, 3, 2, 4, 1])
+        cls.lc2 = Lightcurve([1, 2, 3, 4, 5], [4, 8, 1, 9, 11])
+        # Smaller Light curve
+        cls.lc_s = Lightcurve([1, 2, 3], [5, 3, 2])
+        # lc with different time resolution
+        cls.lc_u = Lightcurve([1, 3, 5, 7, 9], [4, 8, 1, 9, 11])
+
+    def test_empty_cross_correlation(self):
+        cr = CrossCorrelation()
+        assert cr.lc1 is None
+        assert cr.lc2 is None
+        assert cr.corr is None
+        assert cr.time_shift is None
+        assert cr.time_lags is None
+        assert cr.dt is None
+        assert cr.n is None
+
+    def test_empty_cross_correlation_with_dt(self):
+        cr = CrossCorrelation()
+        with pytest.raises(StingrayError):
+            cr.cal_timeshift(dt=2.0)
+        assert np.isclose(cr.dt, 2.0)
+
+    def test_init_with_only_one_lc(self):
+        with pytest.raises(TypeError):
+            CrossCorrelation(self.lc1)
+
+    def test_init_with_invalid_lc1(self):
+        data = np.array([[2, 3, 2, 4, 1]])
+        with pytest.raises(TypeError):
+            CrossCorrelation(data, self.lc2)
+
+    def test_init_with_invalid_lc2(self):
+        data = np.array([[2, 3, 2, 4, 1]])
+        with pytest.raises(TypeError):
+            CrossCorrelation(self.lc1, data)
+
+    def test_init_with_diff_time_bin(self):
+        with pytest.raises(StingrayError):
+            CrossCorrelation(self.lc_u, self.lc2)
+
+    def test_corr_is_correct(self):
+        result = np.array([1.92, 2.16, 1.8, -14.44, 11.12])
+        lags_result = np.array([-2, -1, 0, 1, 2])
+        cr = CrossCorrelation(self.lc1, self.lc2)
+        assert np.allclose(cr.lc1, self.lc1)
+        assert np.allclose(cr.lc2, self.lc2)
+        assert np.allclose(cr.corr, result)
+        assert np.isclose(cr.dt, self.lc1.dt)
+        assert cr.n == 5
+        assert np.allclose(cr.time_lags, lags_result)
+        assert np.isclose(cr.time_shift, 2.0)
+        assert cr.mode == 'same'
+
+    def test_cross_correlation_with_unequal_lc(self):
+        result = np.array([-0.66666667, -0.33333333, -1., 0.66666667, 3.13333333])
+        lags_result = np.array([-2, -1, 0, 1, 2])
+        cr = CrossCorrelation(self.lc1, self.lc_s)
+        assert np.allclose(cr.lc1, self.lc1)
+        assert np.allclose(cr.lc2, self.lc_s)
+        assert np.allclose(cr.corr, result)
+        assert np.isclose(cr.dt, self.lc1.dt)
+        assert cr.n == 5
+        assert np.allclose(cr.time_lags, lags_result)
+        assert np.isclose(cr.time_shift, 2.0)
+        assert cr.mode == 'same'
+
+    def test_mode_with_bad_input(self):
+        with pytest.raises(TypeError):
+            CrossCorrelation(self.lc1, self.lc2, mode=123)
+
+    def test_mode_with_wrong_input(self):
+        with pytest.raises(ValueError):
+            CrossCorrelation(self.lc1, self.lc2, mode='default')
+
+    def test_full_mode_is_correct(self):
+        result = np.array([-1.76, 1.68, 1.92, 2.16, 1.8, -14.44, 11.12, -6.12, 3.64])
+        lags_result = np.array([-4, -3, -2, -1, 0, 1, 2, 3, 4])
+        cr = CrossCorrelation(self.lc1, self.lc2, mode='full')
+        assert np.allclose(cr.lc1, self.lc1)
+        assert np.allclose(cr.lc2, self.lc2)
+        assert cr.mode == 'full'
+        assert cr.n == 9
+        assert np.allclose(cr.corr, result)
+        assert np.allclose(cr.time_lags, lags_result)
+        assert np.isclose(cr.time_shift, 2.0)
+
+    def test_timeshift_with_no_corr_but_lc_assigned(self):
+        result = np.array([1.92, 2.16, 1.8, -14.44, 11.12])
+        lags_result = np.array([-2, -1, 0, 1, 2])
+        cr = CrossCorrelation()
+        cr.lc1 = self.lc1
+        cr.lc2 = self.lc2
+        cr.cal_timeshift(dt=1.0)
+        assert np.allclose(cr.lc1, self.lc1)
+        assert np.allclose(cr.lc2, self.lc2)
+        assert np.allclose(cr.corr, result)
+        assert np.isclose(cr.dt, self.lc1.dt)
+        assert cr.n == 5
+        assert np.allclose(cr.time_lags, lags_result)
+        assert np.isclose(cr.time_shift, 2.0)
+        assert cr.mode == 'same'
+
+    def test_timeshift_with_corr_and_lc_assigned(self):
+        result = np.array([1.92, 2.16, 1.8, -14.44, 11.12])
+        lags_result = np.array([-2, -1, 0, 1, 2])
+        cr = CrossCorrelation()
+        cr.lc1 = self.lc1
+        cr.lc2 = self.lc2
+        cr.corr = result
+        cr.cal_timeshift(dt=1.0)
+        assert np.allclose(cr.lc1, self.lc1)
+        assert np.allclose(cr.lc2, self.lc2)
+        assert np.allclose(cr.corr, result)
+        assert np.isclose(cr.dt, self.lc1.dt)
+        assert cr.n == 5
+        assert np.allclose(cr.time_lags, lags_result)
+        assert np.isclose(cr.time_shift, 2.0)
+        assert cr.mode == 'same'
+
+    @pytest.mark.skipif(HAS_MPL, reason='Matplotlib is already installed if condition is met')
+    def test_plot_matplotlib_not_installed(self):
+        cr = CrossCorrelation(self.lc1, self.lc2)
+        with pytest.raises(ImportError) as excinfo:
+            cr.plot()
+        message = str(excinfo.value)
+        assert "Matplotlib required for plot()" in message
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_simple_plot(self):
+        cr = CrossCorrelation(self.lc1, self.lc2)
+        cr.plot()
+        assert plt.fignum_exists(1)
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_wrong_label_type(self):
+        cr = CrossCorrelation(self.lc1, self.lc2)
+        with pytest.raises(TypeError):
+            with warnings.catch_warnings(record=True) as w:
+                cr.plot(labels=123)
+                assert "must be either a list or tuple" in str(w[0].message)
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_labels_index_error(self):
+        cr = CrossCorrelation(self.lc1, self.lc2)
+        with warnings.catch_warnings(record=True) as w:
+            cr.plot(labels='x')
+            assert "must have two labels" in str(w[0].message)
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_axis(self):
+        cr = CrossCorrelation(self.lc1, self.lc2)
+        cr.plot(axis=[0, 1, 0, 100])
+        assert plt.fignum_exists(1)
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_title(self):
+        cr = CrossCorrelation(self.lc1, self.lc2)
+        cr.plot(title="Test for Cross Correlation")
+        assert plt.fignum_exists(1)
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_default_filename(self):
+        cr = CrossCorrelation(self.lc1, self.lc2)
+        cr.plot(save=True)
+        assert os.path.isfile('corr.png')
+        os.unlink('corr.png')
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_custom_filename(self):
+        cr = CrossCorrelation(self.lc1, self.lc2)
+        cr.plot(save=True, filename='cr.png')
+        assert os.path.isfile('cr.png')
+        os.unlink('cr.png')
diff --git a/stingray/tests/test_events.py b/stingray/tests/test_events.py
index 13915309..63823177 100644
--- a/stingray/tests/test_events.py
+++ b/stingray/tests/test_events.py
@@ -28,10 +28,10 @@ class TestEvents(object):
 
     def test_inequal_length(self):
         """Check that exception is raised in case of
-        disparity in length of 'time' and 'energy'
+        disparity in length of 'time' and 'pha'
         """
         with pytest.raises(ValueError):
-            EventList(time=[1, 2, 3], energy=[10, 12])
+            EventList(time=[1, 2, 3], pha=[10, 12])    
 
     def test_to_lc(self):
         """Create a light curve from event list."""
@@ -101,17 +101,17 @@ class TestEvents(object):
         ev = EventList()
         ev.simulate_energies(self.spectrum)
 
-    def test_compare_energy(self):
+    def test_compare_pha(self):
         """Compare the simulated energy distribution to actual distribution.
         """
         fluxes = np.array(self.spectrum[1])
         ev = EventList(ncounts=1000)
         ev.simulate_energies(self.spectrum)
-        energy = [int(p) for p in ev.energy]
+        pha = [int(p) for p in ev.pha]
 
-        # Histogram energy to get shape approximation
-        gen_energy = ((np.array(energy) - 1) / 1).astype(int)
-        lc = np.bincount(energy)
+        # Histogram pha to get shape approximation
+        gen_pha = ((np.array(pha) - 1) / 1).astype(int)
+        lc = np.bincount(pha)
 
         # Remove first entry as it contains occurences of '0' element
         lc = lc[1:len(lc)]
@@ -151,7 +151,7 @@ class TestEvents(object):
         assert ev_new.time == None
         assert ev_new.gti == None
         assert ev_new.pi == None
-        assert ev_new.energy == None
+        assert ev_new.pha == None
 
         ev = EventList(time=[1, 2, 3])
         ev_other = EventList([])
@@ -168,12 +168,12 @@ class TestEvents(object):
         ev_new = ev.join(ev_other)
         assert ev_new.dt == 3
 
-    def test_join_without_energy(self):
-        ev = EventList(time=[1, 2, 3], energy=[3, 3, 3])
+    def test_join_without_pha(self):
+        ev = EventList(time=[1, 2, 3], pha=[3, 3, 3])
         ev_other = EventList(time=[4, 5])
         ev_new = ev.join(ev_other)
 
-        assert np.all(ev_new.energy == [3, 3, 3, 0, 0])
+        assert np.all(ev_new.pha == [3, 3, 3, 0, 0])
 
     def test_join_without_pi(self):
         ev = EventList(time=[1, 2, 3], pi=[3, 3, 3])
@@ -205,14 +205,14 @@ class TestEvents(object):
         """Join two overlapping event lists.
         """
         ev = EventList(time=[1, 1, 2, 3, 4], 
-                        energy=[3, 4, 7, 4, 3], gti=[[1, 2],[3, 4]])
+                        pha=[3, 4, 7, 4, 3], gti=[[1, 2],[3, 4]])
         ev_other = EventList(time=[5, 6, 6, 7, 10], 
-                            energy=[4, 3, 8, 1, 2], gti=[[6, 7]])
+                            pha=[4, 3, 8, 1, 2], gti=[[6, 7]])
         ev_new = ev.join(ev_other)
 
         assert (ev_new.time ==
                 np.array([1, 1, 2, 3, 4, 5, 6, 6, 7, 10])).all()
-        assert (ev_new.energy ==
+        assert (ev_new.pha ==
                 np.array([3, 4, 7, 4, 3, 4, 3, 8, 1, 2])).all()
         assert (ev_new.gti == 
                 np.array([[1, 2],[3, 4],[6, 7]])).all()
@@ -221,14 +221,14 @@ class TestEvents(object):
         """Join two non-overlapping event lists.
         """
         ev = EventList(time=[1, 1, 10, 6, 5], 
-                        energy=[10, 6, 3, 11, 2], gti=[[1, 3],[5, 6]])
+                        pha=[10, 6, 3, 11, 2], gti=[[1, 3],[5, 6]])
         ev_other = EventList(time=[5, 7, 6, 6, 10], 
-                            energy=[2, 3, 8, 1, 2], gti=[[5, 7],[8, 10]])
+                            pha=[2, 3, 8, 1, 2], gti=[[5, 7],[8, 10]])
         ev_new = ev.join(ev_other)
 
         assert (ev_new.time ==
                 np.array([1, 1, 5, 5, 6, 6, 6, 7, 10, 10])).all()
-        assert (ev_new.energy ==
+        assert (ev_new.pha ==
                 np.array([10, 6, 2, 2, 11, 8, 1, 3, 3, 2])).all()
         assert (ev_new.gti == np.array([[5, 6]])).all()
 
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index d50c0e66..1517f55f 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -33,7 +33,7 @@ class TestLightcurve(object):
         Demonstrate that we can create a trivial Lightcurve object.
         """
         lc = Lightcurve(self.times, self.counts)
-	
+    
     def test_irregular_time_warning(self):
         """
         Check if inputting an irregularly spaced time iterable throws out
@@ -614,3 +614,8 @@ class TestLightcurveRebin(object):
         lc = Lightcurve(times, counts, gti=gti)
         baseline = lc.baseline(10000, 0.01)
         assert np.all(lc.counts - baseline < 1)
+
+    def test_change_mjdref(self):
+        lc_new = self.lc.change_mjdref(57000)
+        assert lc_new.mjdref == 57000
+    
\ No newline at end of file
diff --git a/stingray/tests/test_varenergyspectrum.py b/stingray/tests/test_varenergyspectrum.py
index 9785fb32..996d365e 100644
--- a/stingray/tests/test_varenergyspectrum.py
+++ b/stingray/tests/test_varenergyspectrum.py
@@ -42,7 +42,7 @@ class TestVarEnergySpectrum(object):
         alltimes = np.random.uniform(tstart, tend, nphot)
         alltimes.sort()
         cls.events = EventList(alltimes,
-                               energy=np.random.uniform(0.3, 12, nphot),
+                               pha=np.random.uniform(0.3, 12, nphot),
                                gti = [[tstart, tend]])
         cls.vespec = DummyVarEnergy(cls.events, [0., 10000],
                                     (0.5, 5, 10, "lin"), [0.3, 10],
@@ -63,7 +63,7 @@ class TestVarEnergySpectrum(object):
 
     def test_ref_band_none(self):
         events = EventList([0.09, 0.21, 0.23, 0.32, 0.4, 0.54],
-                           energy=[0,0,0,0,1,1],
+                           pha=[0,0,0,0,1,1],
                            gti=[[0, 0.65]])
         vespec = DummyVarEnergy(events, [0., 10000],
                                 (0, 1, 2, "lin"),
@@ -72,7 +72,7 @@ class TestVarEnergySpectrum(object):
 
     def test_energy_spec_wrong_list_not_tuple(self):
         events = EventList([0.09, 0.21, 0.23, 0.32, 0.4, 0.54],
-                           energy=[0, 0, 0, 0, 1, 1],
+                           pha=[0, 0, 0, 0, 1, 1],
                            gti=[[0, 0.65]])
         # Test using a list instead of tuple
         # with pytest.raises(ValueError):
@@ -82,7 +82,7 @@ class TestVarEnergySpectrum(object):
 
     def test_energy_spec_wrong_str(self):
         events = EventList([0.09, 0.21, 0.23, 0.32, 0.4, 0.54],
-                           energy=[0, 0, 0, 0, 1, 1],
+                           pha=[0, 0, 0, 0, 1, 1],
                            gti=[[0, 0.65]])
         # Test using a list instead of tuple
         with pytest.raises(ValueError):
@@ -92,7 +92,7 @@ class TestVarEnergySpectrum(object):
 
     def test_construct_lightcurves(self):
         events = EventList([0.09, 0.21, 0.23, 0.32, 0.4, 0.54],
-                           energy=[0,0,0,0,1,1],
+                           pha=[0,0,0,0,1,1],
                            gti=[[0, 0.65]])
         vespec = DummyVarEnergy(events, [0., 10000],
                                 (0, 1, 2, "lin"), [0.5, 1.1],
@@ -105,7 +105,7 @@ class TestVarEnergySpectrum(object):
 
     def test_construct_lightcurves_no_exclude(self):
         events = EventList([0.09, 0.21, 0.23, 0.32, 0.4, 0.54],
-                           energy=[0,0,0,0,1,1],
+                           pha=[0,0,0,0,1,1],
                            gti=[[0, 0.65]])
 
         vespec = DummyVarEnergy(events, [0., 10000],
@@ -141,8 +141,8 @@ class TestRMSEnergySpectrum(object):
         test_ev1, test_ev2 = EventList(), EventList()
         test_ev1.simulate_times(test_lc)
         test_ev2.simulate_times(test_lc)
-        test_ev1.energy = np.random.uniform(0.3, 12, len(test_ev1.time))
-        test_ev2.energy = np.random.uniform(0.3, 12, len(test_ev2.time))
+        test_ev1.pha = np.random.uniform(0.3, 12, len(test_ev1.time))
+        test_ev2.pha = np.random.uniform(0.3, 12, len(test_ev2.time))
 
         cls.rms = RmsEnergySpectrum(test_ev1, [0., 100],
                                     (0.3, 12, 5, "lin"),
@@ -178,8 +178,8 @@ class TestLagEnergySpectrum(object):
         test_ev1, test_ev2 = EventList(), EventList()
         test_ev1.simulate_times(test_lc1)
         test_ev2.simulate_times(test_lc2)
-        test_ev1.energy = np.random.uniform(0.3, 9, len(test_ev1.time))
-        test_ev2.energy = np.random.uniform(9, 12, len(test_ev2.time))
+        test_ev1.pha = np.random.uniform(0.3, 9, len(test_ev1.time))
+        test_ev2.pha = np.random.uniform(9, 12, len(test_ev2.time))
 
         cls.lag = LagEnergySpectrum(test_ev1, [0., 0.5],
                                     (0.3, 9, 4, "lin"), [9, 12],