def read(self, filename, format_='pickle'):
        """
        Imports LightCurve object.

        Parameters
        ----------
        filename: str
            Name of the LightCurve object to be read.

        format\_: str
            Available options are 'pickle', 'hdf5', 'ascii'

        Returns
        --------
        If format\_ is 'ascii': astropy.table is returned.
        If format\_ is 'hdf5': dictionary with key-value pairs is returned.
        If format\_ is 'pickle': class object is set.
        """

        if format_ == 'ascii' or format_ == 'hdf5':
            return io.read(filename, format_)

        elif format_ == 'pickle':
            self = io.read(filename, format_)

        else:
            utils.simon("Format not understood.")

    def _construct_energy_covar(self, energy_events, energy_covar,
                                xs_var=None):
        """Form the actual output covaraince dictionary and array."""
        self._init_energy_covar(energy_events, energy_covar)

        if not self.avg_covar:
            xs_var = dict()

        for energy in energy_covar.keys():
            lc, lc_ref = self._create_lc_and_lc_ref(energy, energy_events)

            covar = self._compute_covariance(lc, lc_ref)

            energy_covar[energy] = covar
            if not self.avg_covar:
                self.covar_error[energy] = self._calculate_covariance_error(
                                                lc, lc_ref)

            # Excess variance in ref band
            xs_var[energy] = self._calculate_excess_variance(lc_ref)

        for key, value in energy_covar.items():
            if not xs_var[key] > 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

            self.covar = np.vstack(energy_covar.items())

            self.covar_error = np.vstack(self.covar_error.items())

    def write(self, filename, format_='pickle', **kwargs):
        """
        Exports LightCurve object.

        Parameters
        ----------
        filename: str
            Name of the LightCurve object to be created.

        format\_: str
            Available options are 'pickle', 'hdf5', 'ascii'
        """

        if format_ == 'ascii':
            io.write(np.array([self.time, self.counts]).T,
                     filename, format_, fmt=["%s", "%s"])

        elif format_ == 'pickle':
            io.write(self, filename, format_)

        elif format_ == 'hdf5':
            io.write(self, filename, format_)

        else:
            utils.simon("Format not understood.")

    def _spectrum_function(self):

        rms_spec = np.zeros(len(self.energy_intervals))
        rms_spec_err = np.zeros_like(rms_spec)
        for i, eint in enumerate(self.energy_intervals):
            base_lc, ref_lc = self._construct_lightcurves(eint,
                                                          exclude=False)
            try:
                xspect = AveragedCrossspectrum(base_lc, ref_lc,
                                               segment_size=self.segment_size,
                                               norm='frac')
            except AssertionError as e:
                # Avoid "Mean count rate is <= 0. Something went wrong" assertion.
                simon("AssertionError: " + str(e))
            else:
                good = (xspect.freq >= self.freq_interval[0]) & \
                       (xspect.freq < self.freq_interval[1])
                rms_spec[i] = np.sqrt(np.sum(xspect.power[good]*xspect.df))

                # Root squared sum of errors of the spectrum
                root_sq_err_sum = np.sqrt(np.sum((xspect.power_err[good]*xspect.df)**2))
                # But the rms is the squared root. So,
                # Error propagation
                rms_spec_err[i] = 1 / (2 * rms_spec[i]) * root_sq_err_sum

        return rms_spec, rms_spec_err

def read(filename, format_='pickle', **kwargs):
    """
    Return a pickled class instance.

    Parameters
    ----------
    filename: str
        name of the file to be retrieved.
    format_: str
        pickle, hdf5, ascii ...
    """

    if format_ == 'pickle':
        return _retrieve_pickle_object(filename)

    elif format_ == 'hdf5':
        if _H5PY_INSTALLED:
            return _retrieve_hdf5_object(filename)
        else:
            utils.simon('h5py not installed, cannot read an hdf5 object.')

    elif format_ == 'ascii':
        return _retrieve_ascii_object(filename, **kwargs)
    
    else:
        utils.simon('Format not understood.')

    def read(self, filename, format_='pickle'):
        """
        Read a :class:`Lightcurve` object from file. Currently supported formats are

        * pickle (not recommended for long-term storage)
        * HDF5
        * ASCII

        Parameters
        ----------
        filename: str
            Path and file name for the file to be read.

        format\_: str
            Available options are 'pickle', 'hdf5', 'ascii'

        Returns
        --------
        lc : ``astropy.table`` or ``dict`` or :class:`Lightcurve` object
            * If ``format\_`` is ``ascii``: ``astropy.table`` is returned.
            * If ``format\_`` is ``hdf5``: dictionary with key-value pairs is returned.
            * If ``format\_`` is ``pickle``: :class:`Lightcurve` object is returned.
        """

        if format_ == 'ascii' or format_ == 'hdf5':
            return io.read(filename, format_)

        elif format_ == 'pickle':
            self = io.read(filename, format_)

        else:
            utils.simon("Format not understood.")

    def write(self, filename, format_='pickle', **kwargs):
        """
        Write a :class:`Lightcurve` object to file. Currently supported formats are

        * pickle (not recommended for long-term storage)
        * HDF5
        * ASCII

        Parameters
        ----------
        filename: str
            Path and file name for the output file.

        format\_: str
            Available options are 'pickle', 'hdf5', 'ascii'
        """

        if format_ == 'ascii':
            io.write(np.array([self.time, self.counts]).T,
                     filename, format_, fmt=["%s", "%s"])

        elif format_ == 'pickle':
            io.write(self, filename, format_)

        elif format_ == 'hdf5':
            io.write(self, filename, format_)

        else:
            utils.simon("Format not understood.")

def savefig(filename, **kwargs):
    """
    Save a figure plotted by Matplotlib.

    Note : This function is supposed to be used after the ``plot``
    function. Otherwise it will save a blank image with no plot.

    Parameters
    ----------
    filename : str
        The name of the image file. Extension must be specified in the
        file name. For example filename with `.png` extension will give a
        rasterized image while `.pdf` extension will give a vectorized
        output.

    kwargs : keyword arguments
        Keyword arguments to be passed to ``savefig`` function of
        ``matplotlib.pyplot``. For example use `bbox_inches='tight'` to
        remove the undesirable whitepace around the image.
    """

    try:
        import matplotlib.pyplot as plt
    except ImportError:
        raise ImportError("Matplotlib required for savefig()")

    if not plt.fignum_exists(1):
        utils.simon("use ``plot`` function to plot the image first and "
                    "then use ``savefig`` to save the figure.")

    plt.savefig(filename, **kwargs)

def _save_hdf5_object(object, filename):
    """
    Save a class object in hdf5 format.

    Parameters
    ----------
    object: class instance
        A class object whose attributes would be saved in a dictionary format.

    filename: str
        The file name to save to
    """
    items = vars(object)
    attrs = [name for name in items]

    with h5py.File(filename, 'w') as hf:   
        for attr in attrs:
            data = items[attr]
            
            # If data is a single number, store as an attribute.
            if _isattribute(data):

                if isinstance(data, np.longdouble):
                    data = np.double(data) 
                    utils.simon("Casting data as double instead of longdouble.")
                hf.attrs[attr] = data
            
            # If data is a numpy array, create a dataset.
            else:

                if isinstance(data[0], np.longdouble):
                    data = np.double(data) 
                    utils.simon("Casting data as double instead of longdouble.")
                hf.create_dataset(attr, data=data) 

def write(input_, filename, format_='pickle', **kwargs):
    """
    Pickle a class instance. For parameters depending on
    `format_`, see individual function definitions.

    Parameters
    ----------
    object: a class instance
    filename: str
        name of the file to be created.
    format_: str
        pickle, hdf5, ascii ...
    """

    if format_ == 'pickle':
        _save_pickle_object(input_, filename)

    elif format_ == 'hdf5':
        if _H5PY_INSTALLED:
            _save_hdf5_object(input_, filename)

        else:
            utils.simon('h5py not installed, using pickle instead to save object.')
            _save_pickle_object(input_, filename.split('.')[0]+'.pickle')

    elif format_ == 'ascii':
        _save_ascii_object(input_, filename, **kwargs)

    else:
        utils.simon('Format not understood.')

    def _construct_covar(self):
        """
        Helper method to construct the covariance attribute and fill it with values.
        """
        self.avg_covar = True

        start_time = self.lcs[0].time[0]

        covar = np.zeros(len(self.lcs))
        covar_err = np.zeros(len(self.lcs))
        xs_var = np.zeros(len(self.lcs))

        for i in range(len(self.lcs)):
            lc = self.lcs[i]

            if np.size(self.ref_band_lcs) == 1:
                lc_ref = self.ref_band_lcs
            else:
                lc_ref = self.ref_band_lcs[i]

            tstart = start_time
            tend = start_time + self.segment_size
            cv = 0.0
            cv_err = 0.0
            xs = 0.0

            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)

                lc_seg = lc.truncate(start=start_ind, stop=end_ind)
                lc_ref_seg = lc_ref.truncate(start=start_ind, stop=end_ind)

                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!")

                tstart += self.segment_size
                tend += self.segment_size


            covar[i] = cv/self.nbins
            covar_err[i] = cv_err/self.nbins
            xs_var[i] = xs/self.nbins

        self.unnorm_covar = covar
        energy_covar = covar / xs_var**0.5

        self.covar = energy_covar

        self.covar_error = covar_err

        return

    def _operation_with_other_lc(self, other, operation):
        """
        Helper method to codify an operation of one light curve with another (e.g. add, subtract, ...).
        Takes into account the GTIs correctly, and returns a new :class:`Lightcurve` object.

        Parameters
        ----------
        other : :class:`Lightcurve` object
            A second light curve object

        operation : function
            An operation between the :class:`Lightcurve` object calling this method, and ``other``,
            operating on the ``counts`` attribute in each :class:`Lightcurve` object

        Returns
        -------
        lc_new : Lightcurve object
            The new light curve calculated in ``operation``
        """
        if self.mjdref != other.mjdref:
            raise ValueError("MJDref is different in the two light curves")

        common_gti = cross_two_gtis(self.gti, other.gti)
        mask_self = create_gti_mask(self.time, common_gti, dt=self.dt)
        mask_other = create_gti_mask(other.time, common_gti, dt=other.dt)

        # ValueError is raised by Numpy while asserting np.equal over arrays
        # with different dimensions.
        try:
            diff = np.abs((self.time[mask_self] - other.time[mask_other]))
            assert np.all(diff < self.dt / 100)
        except (ValueError, AssertionError):
            raise ValueError("GTI-filtered time arrays of both light curves "
                             "must be of same dimension and equal.")

        new_time = self.time[mask_self]
        new_counts = operation(self.counts[mask_self],
                               other.counts[mask_other])

        if self.err_dist.lower() != other.err_dist.lower():
            simon("Lightcurves have different statistics!"
                  "We are setting the errors to zero to avoid complications.")
            new_counts_err = np.zeros_like(new_counts)
        elif self.err_dist.lower() in valid_statistics:
                new_counts_err = \
                    np.sqrt(np.add(self.counts_err[mask_self]**2,
                                   other.counts_err[mask_other]**2))
            # More conditions can be implemented for other statistics
        else:
            raise StingrayError("Statistics not recognized."
                                " Please use one of these: "
                                "{}".format(valid_statistics))

        lc_new = Lightcurve(new_time, new_counts,
                            err=new_counts_err, gti=common_gti,
                            mjdref=self.mjdref)

        return lc_new

def classical_pvalue(power, nspec):
    """
    Compute the probability of detecting the current power under
    the assumption that there is no periodic oscillation in the data.

    This computes the single-trial p-value that the power was
    observed under the null hypothesis that there is no signal in
    the data.

    Important: the underlying assumptions that make this calculation valid
    are:
    (1) the powers in the power spectrum follow a chi-square distribution
    (2) the power spectrum is normalized according to Leahy (1984), such
    that the powers have a mean of 2 and a variance of 4
    (3) there is only white noise in the light curve. That is, there is no
    aperiodic variability that would change the overall shape of the power
    spectrum.

    Also note that the p-value is for a *single trial*, i.e. the power currently
    being tested. If more than one power or more than one power spectrum are
    being tested, the resulting p-value must be corrected for the number
    of trials (Bonferroni correction).

    Mathematical formulation in Groth, 1975.
    Original implementation in IDL by Anna L. Watts.

    Parameters
    ----------
    power :  float
        The squared Fourier amplitude of a spectrum to be evaluated

    nspec : int
        The number of spectra or frequency bins averaged in `power`.
        This matters because averaging spectra or frequency bins increases
        the signal-to-noise ratio, i.e. makes the statistical distributions
        of the noise narrower, such that a smaller power might be very
        significant in averaged spectra even though it would not be in a single
        power spectrum.

    """

    assert np.isfinite(power), "power must be a finite floating point number!"
    assert power > 0.0, "power must be a positive real number!"
    assert np.isfinite(nspec), "nspec must be a finite integer number"
    assert nspec >= 1, "nspec must be larger or equal to 1"
    assert np.isclose(nspec % 1, 0), "nspec must be an integer number!"

    ## If the power is really big, it's safe to say it's significant,
    ## and the p-value will be nearly zero
    if power*nspec > 30000:
        simon("Probability of no signal too miniscule to calculate.")
        return 0.0

    else:
        pval = _pavnosigfun(power, nspec)
        return pval

    def coherence(self):
        """
        Compute an averaged Coherence function of cross spectrum by computing
        coherence function of each segment and averaging them. The return type
        is a tuple with first element as the coherence function and the second
        element as the corresponding uncertainty[1] associated with it.

        Note : The uncertainty in coherence function is strictly valid for
               Gaussian statistics only.

        Returns
        -------
        tuple : tuple of np.ndarray
            Tuple of coherence function and uncertainty.

        References
        ----------
        .. [1] http://iopscience.iop.org/article/10.1086/310430/pdf

        """
        if self.m < 50:
            utils.simon("Number of segments used in averaging is "
                        "significantly low. The result might not follow the "
                        "expected statistical distributions.")

        # Calculate average coherence
        unnorm_power_avg = np.zeros_like(self.cs_all[0].unnorm_power)
        for cs in self.cs_all:
            unnorm_power_avg += cs.unnorm_power

        unnorm_power_avg /= self.m
        num = np.abs(unnorm_power_avg)**2

        # this computes the averaged power spectrum, but using the
        # cross spectrum code to avoid circular imports
        aps1 = AveragedCrossspectrum(self.lc1, self.lc1,
                                     segment_size=self.segment_size)
        aps2 = AveragedCrossspectrum(self.lc2, self.lc2,
                                     segment_size=self.segment_size)

        unnorm_powers_avg_1 = np.zeros_like(aps1.cs_all[0].unnorm_power)
        for ps in aps1.cs_all:
            unnorm_powers_avg_1 += ps.unnorm_power
        unnorm_powers_avg_1 /= aps1.m

        unnorm_powers_avg_2 = np.zeros_like(aps2.cs_all[0].unnorm_power)
        for ps in aps2.cs_all:
            unnorm_powers_avg_2 += ps.unnorm_power
        unnorm_powers_avg_2 /= aps2.m

        coh = num / (unnorm_powers_avg_1 * unnorm_powers_avg_2)

        # Calculate uncertainty
        uncertainty = (2**0.5 * coh * (1 - coh)) / (np.abs(coh) * self.m**0.5)

        return (coh, uncertainty)

    def coherence(self):
        """
        Compute an averaged Coherence function of cross spectrum by computing
        coherence function of each segment and averaging them. The return type
        is a tuple with first element as the coherence function and the second
        element as the corresponding uncertainty[1] associated with it.

        Note : The uncertainty in coherence function is strictly valid for
               Gaussian statistics only.

        Returns
        -------
        tuple : tuple of np.ndarray
            Tuple of coherence function and uncertainty.

        References
        ----------
        .. [1] http://iopscience.iop.org/article/10.1086/310430/pdf

        """
        if np.any(self.m < 50):
            simon("Number of segments used in averaging is "
                  "significantly low. The result might not follow the "
                  "expected statistical distributions.")

        # Calculate average coherence
        unnorm_power_avg = self.unnorm_power

        num = np.absolute(unnorm_power_avg) ** 2

        # The normalization was 'none'!
        unnorm_powers_avg_1 = self.pds1.power.real
        unnorm_powers_avg_2 = self.pds2.power.real

        coh = num / (unnorm_powers_avg_1 * unnorm_powers_avg_2)

        # Calculate uncertainty
        uncertainty = (2 ** 0.5 * coh * (1 - coh)) / (
        np.abs(coh) * self.m ** 0.5)

        return (coh, uncertainty)

def read(filename, format_='pickle', **kwargs):
    """
    Return a saved class instance.

    Parameters
    ----------
    filename: str
        The name of the file to be retrieved.

    format_: str
        The format used to store file. Supported formats are
        pickle, hdf5, ascii or fits.

    Returns
    -------
    data : {``object`` | ``astropy.table`` | ``dict``}

        * If ``format_`` is ``pickle``, an object is returned.
        * If ``format_`` is ``ascii``, `astropy.table` object is returned.
        * If ``format_`` is ``hdf5`` or 'fits``, a dictionary object is returned.
    """

    if format_ == 'pickle':
        return _retrieve_pickle_object(filename)

    elif format_ == 'hdf5':
        if _H5PY_INSTALLED:
            return _retrieve_hdf5_object(filename)
        else:
            utils.simon('h5py not installed, cannot read an'
                        'hdf5 object.')

    elif format_ == 'ascii':
        return _retrieve_ascii_object(filename, **kwargs)

    elif format_ == 'fits':
        return _retrieve_fits_object(filename, **kwargs)

    else:
        utils.simon('Format not understood.')

    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()]

        self.event_list = event_list

        self.event_list_T = event_list.T

        self._init_special_vars()

        if ref_band_interest is None:
            ref_band_interest = (self.min_energy, self.max_energy)

        assert type(ref_band_interest) in (list, tuple), "Ref Band interest " \
                                                         "should be either " \
                                                         "tuple or list."

        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

        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

        self.std = std

    def _operation_with_other_lc(self, other, operation):
        if self.mjdref != other.mjdref:
            raise ValueError("MJDref is different in the two light curves")

        common_gti = cross_two_gtis(self.gti, other.gti)
        mask_self = create_gti_mask(self.time, common_gti)
        mask_other = create_gti_mask(other.time, common_gti)

        # ValueError is raised by Numpy while asserting np.equal over arrays
        # with different dimensions.
        try:
            assert np.all(np.equal(self.time[mask_self],
                                   other.time[mask_other]))
        except (ValueError, AssertionError):
            raise ValueError("GTI-filtered time arrays of both light curves "
                             "must be of same dimension and equal.")

        new_time = self.time[mask_self]
        new_counts = operation(self.counts[mask_self],
                               other.counts[mask_other])

        if self.err_dist.lower() != other.err_dist.lower():
            simon("Lightcurves have different statistics!"
                  "We are setting the errors to zero to avoid complications.")
            new_counts_err = np.zeros_like(new_counts)
        elif self.err_dist.lower() in valid_statistics:
                new_counts_err = np.sqrt(np.add(self.counts_err[mask_self]**2,
                                                other.counts_err[mask_other]**2))
            # More conditions can be implemented for other statistics
        else:
            raise StingrayError("Statistics not recognized."
                                " Please use one of these: "
                                "{}".format(valid_statistics))

        lc_new = Lightcurve(new_time, new_counts,
                            err=new_counts_err, gti=common_gti,
                            mjdref=self.mjdref)

        return lc_new

def write(input_, filename, format_='pickle', **kwargs):
    """
    Pickle a class instance. For parameters depending on
    ``format_``, see individual function definitions.

    Parameters
    ----------
    object: a class instance
        The object to be stored

    filename: str
        The name of the file to be created

    format_: str
        The format in which to store file. Formats supported
        are ``pickle``, ``hdf5``, ``ascii`` or ``fits``
    """

    if format_ == 'pickle':
        _save_pickle_object(input_, filename)

    elif format_ == 'hdf5':
        if _H5PY_INSTALLED:
            _save_hdf5_object(input_, filename)
        else:
            utils.simon('h5py not installed, using pickle instead'
                        'to save object.')
            _save_pickle_object(input_, filename.split('.')[0] +
                                '.pickle')

    elif format_ == 'ascii':
        _save_ascii_object(input_, filename, **kwargs)

    elif format_ == 'fits':
        _save_fits_object(input_, filename, **kwargs)

    else:
        utils.simon('Format not understood.')

def read(filename, format_='pickle', **kwargs):
    """
    Return a pickled class instance.

    Parameters
    ----------
    filename: str
        The name of the file to be retrieved.

    format_: str
        The format used to store file. Supported formats are
        pickle, hdf5, ascii or fits.
    
    Returns
    -------
    If format_ is 'pickle', a class object is returned.
    If format_ is 'ascii', astropy.table object is returned.
    If format_ is 'hdf5' or 'fits', a dictionary object is returned.
    """

    if format_ == 'pickle':
        return _retrieve_pickle_object(filename)

    elif format_ == 'hdf5':
        if _H5PY_INSTALLED:
            return _retrieve_hdf5_object(filename)
        else:
            utils.simon('h5py not installed, cannot read an' \
                'hdf5 object.')

    elif format_ == 'ascii':
        return _retrieve_ascii_object(filename, **kwargs)

    elif format_ == 'fits':
        return _retrieve_fits_object(filename, **kwargs)
    
    else:
        utils.simon('Format not understood.')