class StingrayError(Exception):
    pass def bin_intervals_from_gtis(gtis, segment_size, time, dt=None, fraction_step=1,
                            epsilon=0.001):
    """Compute start/stop times of equal time intervals, compatible with GTIs, and map them
    to the indices of an array of time stamps.

    Used to start each FFT/PDS/cospectrum from the start of a GTI,
    and stop before the next gap in data (end of GTI).
    In this case, it is necessary to specify the time array containing the
    times of the light curve bins.
    Returns start and stop bins of the intervals to use for the PDS

    Parameters
    ----------
    gtis : 2-d float array
        List of GTIs of the form ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``

    segment_size : float
        Length of each time segment

    time : array-like
        Array of time stamps

    Other Parameters
    ----------------
    dt : float, default median(diff(time))
        Time resolution of the light curve.

    epsilon : float, default 0.001
        The tolerance, in fraction of ``dt``, for the comparisons at the borders

    fraction_step : float
        If the step is not a full ``segment_size`` but less (e.g. a moving window),
        this indicates the ratio between step step and ``segment_size`` (e.g.
        ``0.5`` means that the window shifts of half ``segment_size``)

    Returns
    -------
    spectrum_start_bins : array-like
        List of starting bins in the original time array to use in spectral
        calculations.

    spectrum_stop_bins : array-like
        List of end bins to use in the spectral calculations.

    Examples
    --------
    >>> time = np.arange(0.5, 13.5)

    >>> gtis = [[0, 5], [6, 8], [9, 10]]

    >>> segment_size = 2

    >>> start_bins, stop_bins = bin_intervals_from_gtis(gtis,segment_size,time)

    >>> np.allclose(start_bins, [0, 2, 6])
    True
    >>> np.allclose(stop_bins, [2, 4, 8])
    True
    >>> np.allclose(time[start_bins[0]:stop_bins[0]], [0.5, 1.5])
    True
    >>> np.allclose(time[start_bins[1]:stop_bins[1]], [2.5, 3.5])
    True
    """
    time = np.asarray(time)
    gtis = np.asarray(gtis)
    if dt is None:
        dt = np.median(np.diff(time))

    epsilon_times_dt = epsilon * dt
    nbin = np.rint(segment_size / dt).astype(int)

    if time[-1] < np.min(gtis) or time[0] > np.max(gtis):
        raise ValueError("Invalid time interval for the given GTIs")

    spectrum_start_bins = np.array([], dtype=int)

    gti_low = gtis[:, 0] + dt / 2 - epsilon_times_dt
    gti_up = gtis[:, 1] - dt / 2 + epsilon_times_dt

    for g0, g1 in zip(gti_low, gti_up):
        if (g1 - g0 + dt + epsilon_times_dt) < segment_size:
            continue
        startbin, stopbin = np.searchsorted(time, [g0, g1], "left")
        stopbin += 1
        if stopbin > time.size:
            stopbin = time.size

        if time[startbin] < g0:
            startbin += 1
        # Would be g[1] - dt/2, but stopbin is the end of an interval
        # so one has to add one bin
        if time[stopbin - 1] > g1:
            stopbin -= 1

        newbins = calculate_segment_bin_start(
            startbin, stopbin, nbin, fraction_step=fraction_step)
        spectrum_start_bins = \
            np.append(spectrum_start_bins,
                      newbins)
    assert len(spectrum_start_bins) > 0, \
        ("No GTIs are equal to or longer than segment_size.")
    return spectrum_start_bins, spectrum_start_bins + nbin def check_gtis(gti):
    """Check if GTIs are well-behaved.

    Check that:

    1. the shape of the GTI array is correct;
    2. no start > end
    3. no overlaps.

    Parameters
    ----------
    gti : list
        A list of GTI ``(start, stop)`` pairs extracted from the FITS file.

    Raises
    ------
    TypeError
        If GTIs are of the wrong shape
    ValueError
        If GTIs have overlapping or displaced values
    """
    if len(gti) < 1:
        raise ValueError("Empty GTIs")

    for g in gti:
        if np.size(g) != 2 or np.ndim(g) != 1:
            raise TypeError(
                "Please check the formatting of GTIs. They need to be"
                " provided as [[gti00, gti01], [gti10, gti11], ...]")

    gti = np.array(gti)
    gti_start = gti[:, 0]
    gti_end = gti[:, 1]

    # Check that GTIs are well-behaved
    if not np.all(gti_end >= gti_start):
        raise ValueError('This GTI end times must be larger than '
                         'GTI start times')

    # Check that there are no overlaps in GTIs
    if not np.all(gti_start[1:] >= gti_end[:-1]):
        raise ValueError('This GTI has overlaps')

    return def cross_two_gtis(gti0, gti1):
    """Extract the common intervals from two GTI lists *EXACTLY*.

    Parameters
    ----------
    gti0 : iterable of the form ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``
    gti1 : iterable of the form ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``
        The two lists of GTIs to be crossed.

    Returns
    -------
    gtis : ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``
        The newly created GTIs

    See Also
    --------
    cross_gtis : From multiple GTI lists, extract common intervals *EXACTLY*

    Examples
    --------
    >>> gti1 = np.array([[1, 2]])
    >>> gti2 = np.array([[1, 2]])
    >>> newgti = cross_two_gtis(gti1, gti2)
    >>> np.allclose(newgti, [[1, 2]])
    True
    >>> gti1 = np.array([[1, 4]])
    >>> gti2 = np.array([[1, 2], [2, 4]])
    >>> newgti = cross_two_gtis(gti1, gti2)
    >>> np.allclose(newgti, [[1, 4]])
    True
    """
    gti0 = join_equal_gti_boundaries(np.asarray(gti0))
    gti1 = join_equal_gti_boundaries(np.asarray(gti1))
    # Check GTIs
    check_gtis(gti0)
    check_gtis(gti1)

    gti0_start = gti0[:, 0]
    gti0_end = gti0[:, 1]
    gti1_start = gti1[:, 0]
    gti1_end = gti1[:, 1]

    # Create a list that references to the two start and end series
    gti_start = [gti0_start, gti1_start]
    gti_end = [gti0_end, gti1_end]

    # Concatenate the series, while keeping track of the correct origin of
    # each start and end time
    gti0_tag = np.array([0 for g in gti0_start], dtype=bool)
    gti1_tag = np.array([1 for g in gti1_start], dtype=bool)
    conc_start = np.concatenate((gti0_start, gti1_start))
    conc_end = np.concatenate((gti0_end, gti1_end))
    conc_tag = np.concatenate((gti0_tag, gti1_tag))

    # Put in time order
    order = np.argsort(conc_end)
    conc_start = conc_start[order]
    conc_end = conc_end[order]
    conc_tag = conc_tag[order]

    last_end = conc_start[0] - 1
    final_gti = []
    for ie, e in enumerate(conc_end):
        # Is this ending in series 0 or 1?
        this_series = int(conc_tag[ie])
        other_series = int(this_series == 0)

        # Check that this closes intervals in both series.
        # 1. Check that there is an opening in both series 0 and 1 lower than e
        try:
            st_pos = \
                np.argmax(gti_start[this_series][gti_start[this_series] < e])
            so_pos = \
                np.argmax(gti_start[other_series][gti_start[other_series] < e])
            st = gti_start[this_series][st_pos]
            so = gti_start[other_series][so_pos]

            s = np.max([st, so])
        except:  # pragma: no cover
            continue

        # If this start is inside the last interval (It can happen for equal
        # GTI start times between the two series), then skip!
        if s <= last_end:
            continue
        # 2. Check that there is no closing before e in the "other series",
        # from intervals starting either after s, or starting and ending
        # between the last closed interval and this one
        cond1 = (gti_end[other_series] > s) * (gti_end[other_series] < e)
        cond2 = gti_end[other_series][so_pos] < s
        condition = np.any(np.logical_or(cond1, cond2))
        # Well, if none of the conditions at point 2 apply, then you can
        # create the new gti!
        if not condition:
            final_gti.append([s, e])
            last_end = e

    return np.array(final_gti) def genDataPath(dir_path):
    """Generates data path to chunks.

    Parameters
    ----------
    dir_path: string
        Path to zarr datastore + Top level directory name for data

    Returns
    -------
    list
        List of path's to datastore

    Raises
    ------
    IOError
        If directory does not exist
    """
    path_list = []
    if os.path.isdir(dir_path):
        if not (os.path.isdir(os.path.join(dir_path, 'main_data/'))
                or os.path.join(dir_path, 'meta_data/')):
            raise IOError(("Directory does not exist."))

        else:
            path_list.append(os.path.join(dir_path, 'main_data/'))
            path_list.append(os.path.join(dir_path, 'meta_data/'))

            return path_list

    else:
        raise IOError(("Directory does not exist.")) def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
    """Rebin some data to an arbitrary new data resolution. Either sum
    the data points in the new bins or average them.

    Parameters
    ----------
    x: iterable
        The dependent variable with some resolution, which can vary throughout
        the time series.

    y: iterable
        The independent variable to be binned

    dx_new: float
        The new resolution of the dependent variable ``x``

    Other parameters
    ----------------
    yerr: iterable, optional
        The uncertainties of ``y``, to be propagated during binning.

    method: {``sum`` | ``average`` | ``mean``}, optional, default ``sum``
        The method to be used in binning. Either sum the samples ``y`` in
        each new bin of ``x``, or take the arithmetic mean.

    dx: float
        The old resolution (otherwise, calculated from difference between
        time bins)

    Returns
    -------
    xbin: numpy.ndarray
        The midpoints of the new bins in ``x``

    ybin: numpy.ndarray
        The binned quantity ``y``

    ybin_err: numpy.ndarray
        The uncertainties of the binned values of ``y``.

    step_size: float
        The size of the binning step

    Examples
    --------
    >>> x = np.arange(0, 100, 0.01)
    >>> y = np.ones(x.size)
    >>> yerr = np.ones(x.size)
    >>> xbin, ybin, ybinerr, step_size = rebin_data(
    ...     x, y, 4, yerr=yerr, method='sum', dx=0.01)
    >>> np.allclose(ybin, 400)
    True
    >>> np.allclose(ybinerr, 20)
    True
    >>> xbin, ybin, ybinerr, step_size = rebin_data(
    ...     x, y, 4, yerr=yerr, method='mean')
    >>> np.allclose(ybin, 1)
    True
    >>> np.allclose(ybinerr, 0.05)
    True
    """

    y = np.asarray(y)
    if yerr is None:
        yerr = np.zeros_like(y)
    else:
        yerr = np.asarray(yerr)

    if not dx:
        dx_old = np.diff(x)
    elif np.size(dx) == 1:
        dx_old = np.array([dx])
    else:
        dx_old = dx

    if np.any(dx_new < dx_old):
        raise ValueError("New frequency resolution must be larger than "
                         "old frequency resolution.")

    # left and right bin edges
    # assumes that the points given in `x` correspond to
    # the left bin edges
    xedges = np.hstack([x, x[-1]+dx_old[-1]])

    # new regularly binned resolution
    xbin = np.arange(xedges[0], xedges[-1]+dx_new, dx_new)

    output = np.zeros(xbin.shape[0] - 1, dtype=type(y[0]))
    outputerr = np.zeros(xbin.shape[0] - 1, dtype=type(yerr[0]))
    step_size = np.zeros(xbin.shape[0] - 1)

    all_x = np.searchsorted(xedges, xbin)
    min_inds = all_x[:-1]
    max_inds = all_x[1:]
    xmins = xbin[:-1]
    xmaxs = xbin[1:]
    for i, (xmin, xmax, min_ind, max_ind) in enumerate(zip(xmins, xmaxs, min_inds, max_inds)):
        filtered_y = y[min_ind:max_ind-1]
        filtered_yerr = yerr[min_ind:max_ind-1]
        output[i] = np.sum(filtered_y)
        outputerr[i] = np.sum(filtered_yerr)
        step_size[i] = max_ind - 1 - min_ind

        prev_dx = xedges[min_ind] - xedges[min_ind-1]
        prev_frac = (xedges[min_ind] - xmin)/prev_dx
        output[i] += y[min_ind-1]*prev_frac
        outputerr[i] += yerr[min_ind-1]*prev_frac
        step_size[i] += prev_frac

        if not max_ind == xedges.size:
            dx_post = xedges[max_ind] - xedges[max_ind-1]
            post_frac = (xmax-xedges[max_ind-1])/dx_post
            output[i] += y[max_ind-1]*post_frac
            outputerr[i] += yerr[max_ind-1]*post_frac
            step_size[i] += post_frac

    if method in ['mean', 'avg', 'average', 'arithmetic mean']:
        ybin = output / step_size
        ybinerr = np.sqrt(outputerr) / step_size

    elif method == "sum":
        ybin = output
        ybinerr = np.sqrt(outputerr)

    else:
        raise ValueError("Method for summing or averaging not recognized. "
                         "Please enter either 'sum' or 'mean'.")

    tseg = x[-1] - x[0] + dx_old[-1]

    if (tseg / dx_new % 1) > 0:
        ybin = ybin[:-1]
        ybinerr = ybinerr[:-1]
        step_size = step_size[:-1]

    dx_var = np.var(dx_old) / np.mean(dx_old)

    if np.size(dx_old) == 1 or dx_var < 1e-6:
        step_size = step_size[0]

    new_x0 = (x[0] - (0.5 * dx_old[0])) + (0.5 * dx_new)
    xbin = np.arange(ybin.shape[0]) * dx_new + new_x0

    return xbin, ybin, ybinerr, step_size def rebin_data_log(x, y, f, y_err=None, dx=None):
    """Logarithmic re-bin of some data. Particularly useful for the power
    spectrum.

    The new dependent variable depends on the previous dependent variable
    modified by a factor f:

    .. math::

        d\\nu_j = d\\nu_{j-1} (1+f)

    Parameters
    ----------
    x: iterable
        The dependent variable with some resolution ``dx_old = x[1]-x[0]``

    y: iterable
        The independent variable to be binned

    f: float
        The factor of increase of each bin wrt the previous one.

    Other Parameters
    ----------------
    yerr: iterable, optional
        The uncertainties of ``y`` to be propagated during binning.

    method: {``sum`` | ``average`` | ``mean``}, optional, default ``sum``
        The method to be used in binning. Either sum the samples ``y`` in
        each new bin of ``x`` or take the arithmetic mean.

    dx: float, optional
        The binning step of the initial ``x``

    Returns
    -------
    xbin: numpy.ndarray
        The midpoints of the new bins in ``x``

    ybin: numpy.ndarray
        The binned quantity ``y``

    ybin_err: numpy.ndarray
        The uncertainties of the binned values of ``y``

    step_size: float
        The size of the binning step
    """

    dx_init = apply_function_if_none(dx, np.diff(x), np.median)
    x = np.asarray(x)
    y = np.asarray(y)
    y_err = np.asarray(apply_function_if_none(y_err, y, np.zeros_like))

    if x.shape[0] != y.shape[0]:
        raise ValueError("x and y must be of the same length!")
    if y.shape[0] != y_err.shape[0]:
        raise ValueError("y and y_err must be of the same length!")

    minx = x[0] * 0.5  # frequency to start from
    maxx = x[-1]  # maximum frequency to end
    binx = [minx, minx + dx_init]  # first
    dx = dx_init  # the frequency resolution of the first bin

    # until we reach the maximum frequency, increase the width of each
    # frequency bin by f
    while binx[-1] <= maxx:
        binx.append(binx[-1] + dx * (1.0 + f))
        dx = binx[-1] - binx[-2]

    binx = np.asarray(binx)

    real = y.real
    real_err = y_err.real
    # compute the mean of the ys that fall into each new frequency bin.
    # we cast to np.double due to scipy's bad handling of longdoubles
    biny, bin_edges, binno = scipy.stats.binned_statistic(
        x.astype(np.double), real.astype(np.double),
        statistic="mean", bins=binx)

    biny_err, bin_edges, binno = scipy.stats.binned_statistic(
        x.astype(np.double), real_err.astype(np.double),
        statistic=_root_squared_mean, bins=binx)

    if np.iscomplexobj(y):
        imag = y.imag
        biny_imag, bin_edges, binno = scipy.stats.binned_statistic(
            x.astype(np.double), imag.astype(np.double),
            statistic="mean", bins=binx)

        biny = biny + 1j * biny_imag

    if np.iscomplexobj(y_err):
        imag_err = y_err.imag

        biny_err_imag, bin_edges, binno = scipy.stats.binned_statistic(
            x.astype(np.double), imag_err.astype(np.double),
            statistic=_root_squared_mean, bins=binx)
        biny_err = biny_err + 1j * biny_err_imag

    # compute the number of powers in each frequency bin
    nsamples = np.array([len(binno[np.where(binno == i)[0]])
                         for i in range(1, np.max(binno) + 1, 1)])

    return binx, biny, biny_err, nsamples def simon(message, **kwargs):
    """The Statistical Interpretation MONitor.

    A warning system designed to always remind the user that Simon
    is watching him/her.

    Parameters
    ----------
    message : string
        The message that is thrown

    kwargs : dict
        The rest of the arguments that are passed to ``warnings.warn``
    """

    warnings.warn("SIMON says: {0}".format(message), **kwargs) class StingrayObject(object):
    """This base class defines some general-purpose utilities.

    The main purpose is to have a consistent mechanism for:

    + round-tripping to and from Astropy Tables and other dataframes

    + round-tripping to files in different formats

    The idea is that any object inheriting :class:`StingrayObject` should,
    just by defining an attribute called ``main_array_attr``, be able to perform
    the operations above, with no additional effort.

    ``main_array_attr`` is, e.g. ``time`` for :class:`EventList` and
    :class:`Lightcurve`, ``freq`` for :class:`Crossspectrum`, ``energy`` for
    :class:`VarEnergySpectrum`, and so on. It is the array with wich all other
    attributes are compared: if they are of the same shape, they get saved as
    columns of the table/dataframe, otherwise as metadata.
    """

    def __init__(cls, *args, **kwargs):
        if not hasattr(cls, "main_array_attr"):
            raise RuntimeError(
                "A StingrayObject needs to have the main_array_attr attribute specified"
            )

    def array_attrs(self):
        """List the names of the array attributes of the Stingray Object.

        By array attributes, we mean the ones with the same size and shape as
        ``main_array_attr`` (e.g. ``time`` in ``EventList``)
        """

        main_attr = getattr(self, getattr(self, "main_array_attr"))
        if main_attr is None:
            return []

        return [
            attr
            for attr in dir(self)
            if (
                isinstance(getattr(self, attr), Iterable)
                and np.shape(getattr(self, attr)) == np.shape(main_attr)
            )
        ]

    def meta_attrs(self):
        """List the names of the meta attributes of the Stingray Object.

        By array attributes, we mean the ones with a different size and shape
        than ``main_array_attr`` (e.g. ``time`` in ``EventList``)
        """
        array_attrs = self.array_attrs()
        return [
            attr
            for attr in dir(self)
            if (
                attr not in array_attrs
                and not attr.startswith("_")
                # Use new assignment expression (PEP 572). I'm testing that
                # self.attribute is not callable, and assigning its value to
                # the variable attr_value for further checks
                and not callable(attr_value := getattr(self, attr))
                # a way to avoid EventLists, Lightcurves, etc.
                and not hasattr(attr_value, "meta_attrs")
            )
        ]

    def get_meta_dict(self):
        """Give a dictionary with all non-None meta attrs of the object."""
        meta_attrs = self.meta_attrs()
        meta_dict = {}
        for key in meta_attrs:
            val = getattr(self, key)
            if val is not None:
                meta_dict[key] = val
        return meta_dict

    def to_astropy_table(self):
        """Create an Astropy Table from a ``StingrayObject``

        Array attributes (e.g. ``time``, ``pi``, ``energy``, etc. for
        ``EventList``) are converted into columns, while meta attributes
        (``mjdref``, ``gti``, etc.) are saved into the ``meta`` dictionary.
        """
        data = {}
        array_attrs = self.array_attrs()

        for attr in array_attrs:
            data[attr] = np.asarray(getattr(self, attr))

        ts = Table(data)

        ts.meta.update(self.get_meta_dict())

        return ts

    @classmethod
    def from_astropy_table(cls, ts):
        """Create a Stingray Object object from data in an Astropy Table.

        The table MUST contain at least a column named like the
        ``main_array_attr``.
        The rest of columns will form the array attributes of the
        new object, while the attributes in ds.attrs will
        form the new meta attributes of the object.

        It is strongly advisable to define such attributes and columns
        using the standard attributes of the wanted StingrayObject (e.g.
        ``time``, ``pi``, etc. for ``EventList``)
        """
        cls = cls()

        if len(ts) == 0:
            # return an empty object
            return cls

        array_attrs = ts.colnames

        # Set the main attribute first
        mainarray = np.array(ts[cls.main_array_attr])
        setattr(cls, cls.main_array_attr, mainarray)

        for attr in array_attrs:
            if attr == cls.main_array_attr:
                continue
            setattr(cls, attr.lower(), np.array(ts[attr]))

        for key, val in ts.meta.items():
            setattr(cls, key.lower(), val)

        return cls

    def to_xarray(self):
        """Create an ``xarray`` Dataset from a `StingrayObject`.

        Array attributes (e.g. ``time``, ``pi``, ``energy``, etc. for
        ``EventList``) are converted into columns, while meta attributes
        (``mjdref``, ``gti``, etc.) are saved into the ``ds.attrs`` dictionary.
        """
        from xarray import Dataset

        data = {}
        array_attrs = self.array_attrs()

        for attr in array_attrs:
            data[attr] = np.asarray(getattr(self, attr))

        ts = Dataset(data)

        ts.attrs.update(self.get_meta_dict())

        return ts

    @classmethod
    def from_xarray(cls, ts):
        """Create a `StingrayObject` from data in an xarray Dataset.

        The dataset MUST contain at least a column named like the
        ``main_array_attr``.
        The rest of columns will form the array attributes of the
        new object, while the attributes in ds.attrs will
        form the new meta attributes of the object.

        It is strongly advisable to define such attributes and columns
        using the standard attributes of the wanted StingrayObject (e.g.
        ``time``, ``pi``, etc. for ``EventList``)
        """
        cls = cls()

        if len(ts[cls.main_array_attr]) == 0:
            # return an empty object
            return cls

        array_attrs = ts.coords

        # Set the main attribute first
        mainarray = np.array(ts[cls.main_array_attr])
        setattr(cls, cls.main_array_attr, mainarray)

        for attr in array_attrs:
            if attr == cls.main_array_attr:
                continue
            setattr(cls, attr, np.array(ts[attr]))

        for key, val in ts.attrs.items():
            if key not in array_attrs:
                setattr(cls, key, val)

        return cls

    def to_pandas(self):
        """Create a pandas ``DataFrame`` from a :class:`StingrayObject`.

        Array attributes (e.g. ``time``, ``pi``, ``energy``, etc. for
        ``EventList``) are converted into columns, while meta attributes
        (``mjdref``, ``gti``, etc.) are saved into the ``ds.attrs`` dictionary.
        """
        from pandas import DataFrame

        data = {}
        array_attrs = self.array_attrs()

        for attr in array_attrs:
            data[attr] = np.asarray(getattr(self, attr))

        ts = DataFrame(data)

        ts.attrs.update(self.get_meta_dict())

        return ts

    @classmethod
    def from_pandas(cls, ts):
        """Create an `StingrayObject` object from data in a pandas DataFrame.

        The dataframe MUST contain at least a column named like the
        ``main_array_attr``.
        The rest of columns will form the array attributes of the
        new object, while the attributes in ds.attrs will
        form the new meta attributes of the object.

        It is strongly advisable to define such attributes and columns
        using the standard attributes of the wanted StingrayObject (e.g.
        ``time``, ``pi``, etc. for ``EventList``)

        """
        cls = cls()

        if len(ts) == 0:
            # return an empty object
            return cls

        array_attrs = ts.columns

        # Set the main attribute first
        mainarray = np.array(ts[cls.main_array_attr])
        setattr(cls, cls.main_array_attr, mainarray)

        for attr in array_attrs:
            if attr == cls.main_array_attr:
                continue
            setattr(cls, attr, np.array(ts[attr]))

        for key, val in ts.attrs.items():
            if key not in array_attrs:
                setattr(cls, key, val)

        return cls

    @classmethod
    def read(cls, filename, fmt=None, format_=None):
        r"""Generic reader for :class`StingrayObject`

        Currently supported formats are

        * pickle (not recommended for long-term storage)
        * any other formats compatible with the writers in
          :class:`astropy.table.Table` (ascii.ecsv, hdf5, etc.)

        Files that need the :class:`astropy.table.Table` interface MUST contain
        at least a column named like the ``main_array_attr``.
        The default ascii format is enhanced CSV (ECSV). Data formats
        supporting the serialization of metadata (such as ECSV and HDF5) can
        contain all attributes such as ``mission``, ``gti``, etc with
        no significant loss of information. Other file formats might lose part
        of the metadata, so must be used with care.

        ..note::

            Complex values can be dealt with out-of-the-box in some formats
            like HDF5 or FITS, not in others (e.g. all ASCII formats).
            With these formats, and in any case when fmt is ``None``, complex
            values should be stored as two columns of real numbers, whose names
            are of the format <variablename>.real and <variablename>.imag

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

        fmt: str
            Available options are 'pickle', 'hea', and any `Table`-supported
            format such as 'hdf5', 'ascii.ecsv', etc.

        Returns
        -------
        obj: :class:`StingrayObject` object
            The object reconstructed from file
        """
        if fmt is None and format_ is not None:
            warnings.warn(
                "The format_ keyword for read and write is deprecated. Use fmt instead",
                DeprecationWarning,
            )
            fmt = format_

        if fmt is None:
            pass
        elif fmt.lower() == "pickle":
            with open(filename, "rb") as fobj:
                return pickle.load(fobj)
        elif fmt.lower() == "ascii":
            fmt = "ascii.ecsv"

        ts = Table.read(filename, format=fmt)

        # For specific formats, and in any case when the format is not
        # specified, make sure that complex values are treated correctly.
        if fmt is None or "ascii" in fmt:
            for col in ts.colnames:
                if not (
                    (is_real := col.endswith(".real"))
                    or (is_imag := col.endswith(".imag"))
                ):
                    continue

                new_value = ts[col]

                if is_imag:
                    new_value = new_value * 1j

                # Make sure it's complex, even if we find the real part first
                new_value = new_value + 0.0j

                col_strip = col.replace(".real", "").replace(".imag", "")

                if col_strip not in ts.colnames:
                    # If the column without ".real" or ".imag" doesn't exist,
                    # define it, and make sure it's complex-valued
                    ts[col_strip] = new_value
                else:
                    # If it does exist, sum the new value to it.
                    ts[col_strip] += new_value

                ts.remove_column(col)

        return cls.from_astropy_table(ts)

    def write(self, filename, fmt=None, format_=None):
        """Generic writer for :class`StingrayObject`

        Currently supported formats are

        * pickle (not recommended for long-term storage)
        * any other formats compatible with the writers in
          :class:`astropy.table.Table` (ascii.ecsv, hdf5, etc.)

        ..note::

            Complex values can be dealt with out-of-the-box in some formats
            like HDF5 or FITS, not in others (e.g. all ASCII formats).
            With these formats, and in any case when fmt is ``None``, complex
            values will be stored as two columns of real numbers, whose names
            are of the format <variablename>.real and <variablename>.imag

        Parameters
        ----------
        filename: str
            Name and path of the file to save the object list to.

        fmt: str
            The file format to store the data in.
            Available options are ``pickle``, ``hdf5``, ``ascii``, ``fits``
        """
        if fmt is None and format_ is not None:
            warnings.warn(
                "The format_ keyword for read and write is deprecated. Use fmt instead",
                DeprecationWarning,
            )
            fmt = format_
        if fmt is None:
            pass
        elif fmt.lower() == "pickle":
            with open(filename, "wb") as fobj:
                pickle.dump(self, fobj)
            return
        elif fmt.lower() == "ascii":
            fmt = "ascii.ecsv"

        ts = self.to_astropy_table()
        if fmt is None or "ascii" in fmt:
            for col in ts.colnames:
                if np.iscomplex(ts[col][0]):
                    ts[f"{col}.real"] = ts[col].real
                    ts[f"{col}.imag"] = ts[col].imag
                    ts.remove_column(col)

        try:
            ts.write(filename, format=fmt, overwrite=True, serialize_meta=True)
        except TypeError:
            ts.write(filename, format=fmt, overwrite=True) class EventList(StingrayTimeseries):
    """
    Basic class for event list data. Event lists generally correspond to individual events (e.g. photons)
    recorded by the detector, and their associated properties. For X-ray data where this type commonly occurs,
    events are time stamps of when a photon arrived in the detector, and (optionally) the photon energy associated
    with the event.

    Parameters
    ----------
    time: iterable
        A list or array of time stamps

    Other Parameters
    ----------------
    dt: float
        The time resolution of the events. Only relevant when using events
        to produce light curves with similar bin time.

    energy: iterable
        A list of array of photon energy values in keV

    mjdref : float
        The MJD used as a reference for the time array.

    ncounts: int
        Number of desired data points in event list.

    gtis: ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``
        Good Time Intervals

    pi : integer, numpy.ndarray
        PI channels

    notes : str
        Any useful annotations

    high_precision : bool
        Change the precision of self.time to float128. Useful while dealing with fast pulsars.

    mission : str
        Mission that recorded the data (e.g. NICER)

    instr : str
        Instrument onboard the mission

    header : str
        The full header of the original FITS file, if relevant

    detector_id : iterable
        The detector that recorded each photon (if the instrument has more than
        one, e.g. XMM/EPIC-pn)

    timeref : str
        The time reference, as recorded in the FITS file (e.g. SOLARSYSTEM)

    timesys : str
        The time system, as recorded in the FITS file (e.g. TDB)

    ephem : str
        The JPL ephemeris used to barycenter the data, if any (e.g. DE430)

    **other_kw :
        Used internally. Any other keyword arguments will be ignored

    Attributes
    ----------
    time: numpy.ndarray
        The array of event arrival times, in seconds from the reference
        MJD defined in ``mjdref``

    energy: numpy.ndarray
        The array of photon energy values

    ncounts: int
        The number of data points in the event list

    dt: float
        The time resolution of the events. Only relevant when using events
        to produce light curves with similar bin time.

    mjdref : float
        The MJD used as a reference for the time array.

    gtis: ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``
        Good Time Intervals

    pi : integer, numpy.ndarray
        PI channels

    high_precision : bool
        Change the precision of self.time to float128. Useful while dealing with fast pulsars.

    mission : str
        Mission that recorded the data (e.g. NICER)

    instr : str
        Instrument onboard the mission

    detector_id : iterable
        The detector that recoded each photon, if relevant (e.g. XMM, Chandra)

    header : str
        The full header of the original FITS file, if relevant

    """
    main_array_attr = "time"
    def __init__(self, time=None, energy=None, ncounts=None, mjdref=0, dt=0,
                 notes="", gti=None, pi=None, high_precision=False,
                 mission=None, instr=None, header=None, detector_id=None,
                 ephem=None, timeref=None, timesys=None,
                 **other_kw):
        StingrayObject.__init__(self)

        self.energy = None if energy is None else np.asarray(energy)
        self.notes = notes
        self.dt = dt
        self.mjdref = mjdref
        self.gti = np.asarray(gti) if gti is not None else None
        self.pi = None if pi is None else np.asarray(pi)
        self.ncounts = ncounts
        self.mission = mission
        self.instr = instr
        self.detector_id = detector_id
        self.header = header
        self.ephem = ephem
        self.timeref = timeref
        self.timesys = timesys

        if other_kw != {}:
            warnings.warn(f"Unrecognized keywords: {list(other_kw.keys())}")

        if time is not None:
            time, mjdref = interpret_times(time, mjdref)
            if not high_precision:
                self.time = np.asarray(time)
            else:
                self.time = np.asarray(time, dtype=np.longdouble)
            self.ncounts = self.time.size
        else:
            self.time = None

        if (self.time is not None) and (self.energy is not None):
            if self.time.size != self.energy.size:
                raise ValueError('Lengths of time and energy must be equal.')

    def to_lc(self, dt, tstart=None, tseg=None):
        """
        Convert event list to a :class:`stingray.Lightcurve` object.

        Parameters
        ----------
        dt: float
            Binning time of the light curve

        Other Parameters
        ----------------
        tstart : float
            Start time of the light curve

        tseg: float
            Total duration of light curve

        Returns
        -------
        lc: :class:`stingray.Lightcurve` object
        """
        if tstart is None and self.gti is not None:
            tstart = self.gti[0][0]
            tseg = self.gti[-1][1] - tstart

        return Lightcurve.make_lightcurve(self.time, dt, tstart=tstart,
                                          gti=self.gti, tseg=tseg,
                                          mjdref=self.mjdref)

    def to_lc_iter(self, dt, segment_size=None):
        """Convert event list to a generator of Lightcurves.

        Parameters
        ----------
        dt: float
            Binning time of the light curves

        Other parameters
        ----------------
        segment_size : float, default None
            Optional segment size. If None, use the GTI boundaries

        Returns
        -------
        lc_gen: `generator`
            Generates one :class:`stingray.Lightcurve` object for each GTI or segment
        """

        segment_iter = generate_indices_of_boundaries(
            self.time, self.gti, segment_size=segment_size, dt=0)

        for st, end, idx_st, idx_end in segment_iter:
            tseg = end - st

            lc = Lightcurve.make_lightcurve(self.time[idx_st:idx_end + 1], dt,
                                            tstart=st,
                                            gti=np.asarray([[st, end]]),
                                            tseg=tseg,
                                            mjdref=self.mjdref, use_hist=True)
            yield lc

    def to_lc_list(self, dt, segment_size=None):
        """Convert event list to a list of Lightcurves.

        Parameters
        ----------
        dt: float
            Binning time of the light curves

        Other parameters
        ----------------
        segment_size : float, default None
            Optional segment size. If None, use the GTI boundaries

        Returns
        -------
        lc_list: `List`
            List containig one :class:`stingray.Lightcurve` object for each GTI or segment
        """
        return list(self.to_lc_iter(dt, segment_size))

    @staticmethod
    def from_lc(lc):
        """
        Create an :class:`EventList` from a :class:`stingray.Lightcurve` object. Note that all
        events in a given time bin will have the same time stamp.

        Parameters
        ----------
        lc: :class:`stingray.Lightcurve` object
            Light curve to use for creation of the event list.

        Returns
        -------
        ev: :class:`EventList` object
            The resulting list of photon arrival times generated from the light curve.
        """

        # Multiply times by number of counts
        times = [[i] * int(j) for i, j in zip(lc.time, lc.counts)]
        # Concatenate all lists
        times = [i for j in times for i in j]

        return EventList(time=times, gti=lc.gti)

    def simulate_times(self, lc, use_spline=False, bin_time=None):
        """Simulate times from an input light curve.

        Randomly simulate photon arrival times to an :class:`EventList` from a
        :class:`stingray.Lightcurve` object, using the inverse CDF method.

        ..note::
            Preferably use model light curves containing **no Poisson noise**,
            as this method will intrinsically add Poisson noise to them.

        Parameters
        ----------
        lc: :class:`stingray.Lightcurve` object

        Other Parameters
        ----------------
        use_spline : bool
            Approximate the light curve with a spline to avoid binning effects

        bin_time : float default None
            Ignored and deprecated, maintained for backwards compatibility.

        Returns
        -------
        times : array-like
            Simulated photon arrival times
        """
        # Need import here, or there will be a circular import
        from .simulator.base import simulate_times
        if bin_time is not None:
            warnings.warn("Bin time will be ignored in simulate_times",
                          DeprecationWarning)

        self.time = simulate_times(lc, use_spline=use_spline)
        self.gti = lc.gti
        self.ncounts = len(self.time)

    def simulate_energies(self, spectrum, use_spline=False):
        """
        Assign (simulate) energies to event list from a spectrum.

        Parameters
        ----------
        spectrum: 2-d array or list [energies, spectrum]
            Energies versus corresponding fluxes. The 2-d array or list must
            have energies across the first dimension and fluxes across the
            second one. If the dimension of the energies is the same as
            spectrum, they are interpreted as bin centers.
            If it is longer by one, they are interpreted as proper bin edges
            (similarly to the bins of `np.histogram`).
            Note that for non-uniformly binned spectra, it is advisable to pass
            the exact edges.
        """
        from .simulator.base import simulate_with_inverse_cdf

        if self.ncounts is None:
            simon("Either set time values or explicity provide counts.")
            return

        if isinstance(spectrum, list) or isinstance(spectrum, np.ndarray):

            energy = np.asarray(spectrum)[0]
            fluxes = np.asarray(spectrum)[1]

            if not isinstance(energy, np.ndarray):
                raise IndexError("Spectrum must be a 2-d array or list")

        else:
            raise TypeError("Spectrum must be a 2-d array or list")

        if energy.size == fluxes.size:
            de = energy[1] - energy[0]
            energy = np.concatenate([energy - de / 2, [energy[-1] + de / 2]])

        self.energy = simulate_with_inverse_cdf(
            fluxes, self.ncounts, edges=energy, sorted=False, interp_kind="linear")

    def sort(self, inplace=False):
        """Sort the event list in time.

        Other parameters
        ----------------
        inplace : bool, default False
            Sort in place. If False, return a new event list.

        Returns
        -------
        eventlist : `EventList`
            The sorted event list. If ``inplace=True``, it will be a shallow copy
            of ``self``.

        Examples
        --------
        >>> events = EventList(time=[0, 2, 1], energy=[0.3, 2, 0.5], pi=[3, 20, 5])
        >>> e1 = events.sort()
        >>> np.allclose(e1.time, [0, 1, 2])
        True
        >>> np.allclose(e1.energy, [0.3, 0.5, 2])
        True
        >>> np.allclose(e1.pi, [3, 5, 20])
        True

        But the original event list has not been altered (``inplace=False`` by
        default):
        >>> np.allclose(events.time, [0, 2, 1])
        True

        Let's do it in place instead
        >>> e2 = events.sort(inplace=True)
        >>> np.allclose(e2.time, [0, 1, 2])
        True

        In this case, the original event list has been altered.
        >>> np.allclose(events.time, [0, 1, 2])
        True

        """
        order = np.argsort(self.time)
        return self.apply_mask(order, inplace=inplace)

    def join(self, other):
        """
        Join two :class:`EventList` objects into one.

        If both are empty, an empty :class:`EventList` is returned.

        GTIs are crossed if the event lists are over a common time interval,
        and appended otherwise.

        ``pi`` and ``pha`` remain ``None`` if they are ``None`` in both. Otherwise, 0 is used
        as a default value for the :class:`EventList` where they were None.

        Parameters
        ----------
        other : :class:`EventList` object
            The other :class:`EventList` object which is supposed to be joined with.

        Returns
        -------
        `ev_new` : :class:`EventList` object
            The resulting :class:`EventList` object.
        """

        ev_new = EventList()

        if self.dt != other.dt:
            simon("The time resolution is different."
                  " Using the rougher by default")
            ev_new.dt = np.max([self.dt, other.dt])

        if self.time is None and other.time is None:
            return ev_new

        if (self.time is None):
            simon("One of the event lists you are concatenating is empty.")
            self.time = np.asarray([])

        elif (other.time is None):
            simon("One of the event lists you are concatenating is empty.")
            other.time = np.asarray([])

        # Tolerance for MJDREF:1 microsecond
        if not np.isclose(self.mjdref, other.mjdref, atol=1e-6 / 86400):
            other = other.change_mjdref(self.mjdref)

        ev_new.time = np.concatenate([self.time, other.time])
        order = np.argsort(ev_new.time)
        ev_new.time = ev_new.time[order]

        if (self.pi is None) and (other.pi is None):
            ev_new.pi = None
        elif (self.pi is None) or (other.pi is None):
            self.pi = assign_value_if_none(self.pi, np.zeros_like(self.time))
            other.pi = assign_value_if_none(other.pi,
                                            np.zeros_like(other.time))

        if (self.pi is not None) and (other.pi is not None):
            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,
                                                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.gti is None and other.gti is not None and len(self.time) > 0:
            self.gti = \
                assign_value_if_none(
                    self.gti, np.asarray([[self.time[0] - self.dt / 2,
                                           self.time[-1] + self.dt / 2]]))
        if other.gti is None and self.gti is not None and len(other.time) > 0:
            other.gti = \
                assign_value_if_none(
                    other.gti, np.asarray([[other.time[0] - other.dt / 2,
                                            other.time[-1] + other.dt / 2]]))

        if (self.gti is None) and (other.gti is None):
            ev_new.gti = None

        elif (self.gti is not None) and (other.gti is not None):
            if check_separate(self.gti, other.gti):
                ev_new.gti = append_gtis(self.gti, other.gti)
                simon('GTIs in these two event lists do not overlap at all.'
                      'Merging instead of returning an overlap.')
            else:
                ev_new.gti = cross_gtis([self.gti, other.gti])

        for attr in ['mission', 'instr']:
            if getattr(self, attr) != getattr(other, attr):
                setattr(ev_new, attr, getattr(self, attr) + ',' + getattr(other, attr))
            else:
                setattr(ev_new, attr, getattr(self, attr))

        ev_new.mjdref = self.mjdref

        return ev_new

    @classmethod
    def read(cls, filename, fmt=None, format_=None, **kwargs):
        r"""Read a :class:`EventList` object from file.

        Currently supported formats are

        * pickle (not recommended for long-term storage)
        * hea or ogip : FITS Event files from (well, some) HEASARC-supported missions.
        * any other formats compatible with the writers in
          :class:`astropy.table.Table` (ascii.ecsv, hdf5, etc.)

        Files that need the :class:`astropy.table.Table` interface MUST contain
        at least a ``time`` column. Other recognized columns are ``energy`` and
        ``pi``.
        The default ascii format is enhanced CSV (ECSV). Data formats
        supporting the serialization of metadata (such as ECSV and HDF5) can
        contain all eventlist attributes such as ``mission``, ``gti``, etc with
        no significant loss of information. Other file formats might lose part
        of the metadata, so must be used with care.

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

        fmt: str
            Available options are 'pickle', 'hea', and any `Table`-supported
            format such as 'hdf5', 'ascii.ecsv', etc.

        Other parameters
        ----------------
        kwargs : dict
            Any further keyword arguments to be passed to `load_events_and_gtis`
            for reading in event lists in OGIP/HEASOFT format

        Returns
        -------
        ev: :class:`EventList` object
            The :class:`EventList` object reconstructed from file
        """
        if fmt is None and format_ is not None:
            warnings.warn(
                "The format_ keyword for read and write is deprecated. "
                "Use fmt instead", DeprecationWarning)
            fmt = format_

        if fmt.lower() in ('hea', 'ogip'):
            evtdata = load_events_and_gtis(filename, **kwargs)

            evt = EventList(time=evtdata.ev_list,
                            gti=evtdata.gti_list,
                            pi=evtdata.pi_list,
                            energy=evtdata.energy_list,
                            mjdref=evtdata.mjdref,
                            instr=evtdata.instr,
                            mission=evtdata.mission,
                            header=evtdata.header,
                            detector_id=evtdata.detector_id,
                            ephem=evtdata.ephem,
                            timeref=evtdata.timeref,
                            timesys=evtdata.timesys)
            if 'additional_columns' in kwargs:
                for key in evtdata.additional_data:
                    if not hasattr(evt, key.lower()):
                        setattr(evt, key.lower(), evtdata.additional_data[key])
            return evt

        return super().read(filename=filename, fmt=fmt)

    def filter_energy_range(self, energy_range, inplace=False, use_pi=False):
        """Filter the event list from a given energy range.

        Parameters
        ----------
        energy_range: [float, float]
            Energy range in keV, or in PI channel (if ``use_pi`` is True)

        Other Parameters
        ----------------
        inplace : bool, default False
            Do the change in place (modify current event list). Otherwise, copy
            to a new event list.
        use_pi : bool, default False
            Use PI channel instead of energy in keV

        Examples
        --------
        >>> events = EventList(time=[0, 1, 2], energy=[0.3, 0.5, 2], pi=[3, 5, 20])
        >>> e1 = events.filter_energy_range([0, 1])
        >>> np.allclose(e1.time, [0, 1])
        True
        >>> np.allclose(events.time, [0, 1, 2])
        True
        >>> e2 = events.filter_energy_range([0, 10], use_pi=True, inplace=True)
        >>> np.allclose(e2.time, [0, 1])
        True
        >>> np.allclose(events.time, [0, 1])
        True

        """
        if use_pi:
            energies = self.pi
        else:
            energies = self.energy
        mask = (energies >= energy_range[0]) & (energies < energy_range[1])

        return self.apply_mask(mask, inplace=inplace)

    def apply_mask(self, mask, inplace=False):
        """Apply a mask to all array attributes of the event list

        Parameters
        ----------
        mask : array of ``bool``
            The mask. Has to be of the same length as ``self.time``

        Other parameters
        ----------------
        inplace : bool
            If True, overwrite the current event list. Otherwise, return a new one.

        Examples
        --------
        >>> evt = EventList(time=[0, 1, 2], mission="nustar")
        >>> evt.bubuattr = [222, 111, 333]
        >>> newev0 = evt.apply_mask([True, True, False], inplace=False);
        >>> newev1 = evt.apply_mask([True, True, False], inplace=True);
        >>> newev0.mission == "nustar"
        True
        >>> np.allclose(newev0.time, [0, 1])
        True
        >>> np.allclose(newev0.bubuattr, [222, 111])
        True
        >>> np.allclose(newev1.time, [0, 1])
        True
        >>> evt is newev1
        True
        """
        array_attrs = self.array_attrs()

        if inplace:
            new_ev = self
        else:
            new_ev = EventList()
            for attr in self.meta_attrs():
                setattr(new_ev, attr, copy.deepcopy(getattr(self, attr)))

        for attr in array_attrs:
            if hasattr(self, attr) and getattr(self, attr) is not None:
                setattr(new_ev, attr, copy.deepcopy(np.asarray(getattr(self, attr))[mask]))
        return new_ev

    def apply_deadtime(self, deadtime, inplace=False, **kwargs):
        """Apply deadtime filter to this event list.

        Additional arguments in ``kwargs`` are passed to `get_deadtime_mask`

        Parameters
        ----------
        deadtime : float
            Value of dead time to apply to data
        inplace : bool, default False
            If True, apply the deadtime to the current event list. Otherwise,
            return a new event list.

        Returns
        -------
        new_event_list : `EventList` object
            Filtered event list. if `inplace` is True, this is the input object
            filtered for deadtime, otherwise this is a new object.
        additional_output : object
            Only returned if `return_all` is True. See `get_deadtime_mask` for
            more details.

        Examples
        --------
        >>> events = np.array([1, 1.05, 1.07, 1.08, 1.1, 2, 2.2, 3, 3.1, 3.2])
        >>> events = EventList(events)
        >>> events.pi=np.array([1, 2, 2, 2, 2, 1, 1, 1, 2, 1])
        >>> events.energy=np.array([1, 2, 2, 2, 2, 1, 1, 1, 2, 1])
        >>> events.mjdref = 10
        >>> filt_events, retval = events.apply_deadtime(0.11, inplace=False,
        ...                                             verbose=False,
        ...                                             return_all=True)
        >>> filt_events is events
        False
        >>> expected = np.array([1, 2, 2.2, 3, 3.2])
        >>> np.allclose(filt_events.time, expected)
        True
        >>> np.allclose(filt_events.pi, 1)
        True
        >>> np.allclose(filt_events.energy, 1)
        True
        >>> np.allclose(events.pi, 1)
        False
        >>> filt_events = events.apply_deadtime(0.11, inplace=True,
        ...                                     verbose=False)
        >>> filt_events is events
        True
        """
        local_retall = kwargs.pop('return_all', False)

        mask, retall = get_deadtime_mask(self.time, deadtime,
                                         return_all=True,
                                         **kwargs)

        new_ev = self.apply_mask(mask, inplace=inplace)

        if local_retall:
            new_ev = [new_ev, retall]

        return new_ev class Lightcurve(StingrayTimeseries):
    """
    Make a light curve object from an array of time stamps and an
    array of counts.

    Parameters
    ----------
    time: Iterable, `:class:astropy.time.Time`, or `:class:astropy.units.Quantity` object
        A list or array of time stamps for a light curve. Must be a type that
        can be cast to `:class:np.array` or `:class:List` of floats, or that
        has a `value` attribute that does (e.g. a
        `:class:astropy.units.Quantity` or `:class:astropy.time.Time` object).

    counts: iterable, optional, default ``None``
        A list or array of the counts in each bin corresponding to the
        bins defined in `time` (note: use ``input_counts=False`` to
        input the count range, i.e. counts/second, otherwise use
        counts/bin).

    err: iterable, optional, default ``None``
        A list or array of the uncertainties in each bin corresponding to
        the bins defined in ``time`` (note: use ``input_counts=False`` to
        input the count rage, i.e. counts/second, otherwise use
        counts/bin). If ``None``, we assume the data is poisson distributed
        and calculate the error from the average of the lower and upper
        1-sigma confidence intervals for the Poissonian distribution with
        mean equal to ``counts``.

    input_counts: bool, optional, default True
        If True, the code assumes that the input data in ``counts``
        is in units of counts/bin. If False, it assumes the data
        in ``counts`` is in counts/second.

    gti: 2-d float array, default ``None``
        ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``
        Good Time Intervals. They are *not* applied to the data by default.
        They will be used by other methods to have an indication of the
        "safe" time intervals to use during analysis.

    err_dist: str, optional, default ``None``
        Statistical distribution used to calculate the
        uncertainties and other statistical values appropriately.
        Default makes no assumptions and keep errors equal to zero.

    mjdref: float
        MJD reference (useful in most high-energy mission data)

    skip_checks: bool
        If True, the user specifies that data are already sorted and contain no
        infinite or nan points. Use at your own risk

    low_memory: bool
        If True, all the lazily evaluated attribute (e.g., countrate and
        countrate_err if input_counts is True) will _not_ be stored in memory,
        but calculated every time they are requested.

    mission : str
        Mission that recorded the data (e.g. NICER)

    instr : str
        Instrument onboard the mission

    header : str
        The full header of the original FITS file, if relevant

    **other_kw :
        Used internally. Any other keyword arguments will be ignored

    Attributes
    ----------
    time: numpy.ndarray
        The array of midpoints of time bins.

    bin_lo: numpy.ndarray
        The array of lower time stamp of time bins.

    bin_hi: numpy.ndarray
        The array of higher time stamp of time bins.

    counts: numpy.ndarray
        The counts per bin corresponding to the bins in ``time``.

    counts_err: numpy.ndarray
        The uncertainties corresponding to ``counts``

    countrate: numpy.ndarray
        The counts per second in each of the bins defined in ``time``.

    countrate_err: numpy.ndarray
        The uncertainties corresponding to ``countrate``

    meanrate: float
        The mean count rate of the light curve.

    meancounts: float
        The mean counts of the light curve.

    n: int
        The number of data points in the light curve.

    dt: float
        The time resolution of the light curve.

    mjdref: float
        MJD reference date (``tstart`` / 86400 gives the date in MJD at the
        start of the observation)

    tseg: float
        The total duration of the light curve.

    tstart: float
        The start time of the light curve.

    gti: 2-d float array
        ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``
        Good Time Intervals. They indicate the "safe" time intervals
        to be used during the analysis of the light curve.

    err_dist: string
        Statistic of the Lightcurve, it is used to calculate the
        uncertainties and other statistical values appropriately.
        It propagates to Spectrum classes.

    mission : str
        Mission that recorded the data (e.g. NICER)

    instr : str
        Instrument onboard the mission

    detector_id : iterable
        The detector that recoded each photon, if relevant (e.g. XMM, Chandra)

    header : str
        The full header of the original FITS file, if relevant

    """
    main_array_attr = "time"

    def __init__(self, time, counts, err=None, input_counts=True,
                 gti=None, err_dist='poisson', mjdref=0, dt=None,
                 skip_checks=False, low_memory=False, mission=None,
                 instr=None, header=None, **other_kw):

        StingrayTimeseries.__init__(self)

        if other_kw != {}:
            warnings.warn(f"Unrecognized keywords: {list(other_kw.keys())}")

        time, mjdref = interpret_times(time, mjdref=mjdref)

        time = np.asarray(time)
        counts = np.asarray(counts)

        if err is not None:
            err = np.asarray(err)

        if not skip_checks:
            time, counts, err = self.initial_optional_checks(time, counts, err)

        if time.size != counts.size:
            raise StingrayError("time and counts array are not "
                                "of the same length!")

        if time.size <= 1:
            raise StingrayError("A single or no data points can not create "
                                "a lightcurve!")

        if err_dist.lower() not in valid_statistics:
            # err_dist set can be increased with other statistics
            raise StingrayError("Statistic not recognized."
                                "Please select one of these: ",
                                "{}".format(valid_statistics))
        elif not err_dist.lower() == 'poisson':
            simon("Stingray only uses poisson err_dist at the moment. "
                  "All analysis in the light curve will assume Poisson "
                  "errors. "
                  "Sorry for the inconvenience.")

        if err is not None:
            err = np.asarray(err)
            if not skip_checks and not np.all(np.isfinite(err)):
                raise ValueError("There are inf or NaN values in "
                                 "your err array")

        self.mjdref = mjdref
        self._time = time

        if dt is None:
            logging.info("Computing the bin time ``dt``. This can take "
                         "time. If you know the bin time, please specify it"
                         " at light curve creation")
            dt = np.median(np.diff(self._time))

        self.dt = dt
        self.err_dist = err_dist

        self.tstart = self._time[0] - 0.5 * self.dt
        self.tseg = self._time[-1] - self._time[0] + self.dt

        self._gti = None
        if gti is not None:
            self._gti = np.asarray(gti)

        self._mask = None
        self._counts = None
        self._counts_err = None
        self._countrate = None
        self._countrate_err = None
        self._meanrate = None
        self._meancounts = None
        self._bin_lo = None
        self._bin_hi = None
        self._n = None
        self.mission = mission
        self.instr = instr
        self.header = header

        self.input_counts = input_counts
        self.low_memory = low_memory
        if input_counts:
            self._counts = np.asarray(counts)
            self._counts_err = err
        else:
            self._countrate = np.asarray(counts)
            self._countrate_err = err

        if not skip_checks:
            self.check_lightcurve()
        if os.name == 'nt':
            warnings.warn(
                "On Windows, the size of an integer is 32 bits. "
                "To avoid integer overflow, I'm converting the input array to float")
            counts = counts.astype(float)

    @property
    def time(self):
        return self._time

    @time.setter
    def time(self, value):
        value = np.asarray(value)
        if not value.shape == self.time.shape:
            raise ValueError('Can only assign new times of the same shape as '
                             'the original array')
        self._time = value
        self._bin_lo = None
        self._bin_hi = None

    @property
    def gti(self):
        if self._gti is None:
            self._gti = \
                np.asarray([[self.tstart, self.tstart + self.tseg]])
        return self._gti

    @gti.setter
    def gti(self, value):
        value = np.asarray(value)
        self._gti = value
        self._mask = None

    @property
    def mask(self):
        if self._mask is None:
            self._mask = create_gti_mask(self.time, self.gti, dt=self.dt)
        return self._mask

    @property
    def n(self):
        if self._n is None:
            self._n = self.counts.shape[0]
        return self._n

    @property
    def meanrate(self):
        if self._meanrate is None:
            self._meanrate = np.mean(self.countrate[self.mask])
        return self._meanrate

    @property
    def meancounts(self):
        if self._meancounts is None:
            self._meancounts = np.mean(self.counts[self.mask])
        return self._meancounts

    @property
    def counts(self):
        counts = self._counts
        if self._counts is None:
            counts = self._countrate * self.dt
            # If not in low-memory regime, cache the values
            if not self.low_memory or self.input_counts:
                self._counts = counts

        return counts

    @counts.setter
    def counts(self, value):
        value = np.asarray(value)
        if not value.shape == self.counts.shape:
            raise ValueError('Can only assign new counts array of the same '
                             'shape as the original array')
        self._counts = value
        self._countrate = None
        self._meancounts = None
        self._meancountrate = None
        self.input_counts = True

    @property
    def counts_err(self):
        counts_err = self._counts_err
        if counts_err is None and self._countrate_err is not None:
            counts_err = self._countrate_err * self.dt
        elif counts_err is None:
            if self.err_dist.lower() == 'poisson':
                counts_err = poisson_symmetrical_errors(self.counts)
            else:
                counts_err = np.zeros_like(self.counts)

        # If not in low-memory regime, cache the values ONLY if they have
        # been changed!
        if self._counts_err is not counts_err:
            if not self.low_memory or self.input_counts:
                self._counts_err = counts_err

        return counts_err

    @counts_err.setter
    def counts_err(self, value):
        value = np.asarray(value)
        if not value.shape == self.counts.shape:
            raise ValueError('Can only assign new error array of the same '
                             'shape as the original array')
        self._counts_err = value
        self._countrate_err = None

    @property
    def countrate(self):
        countrate = self._countrate
        if countrate is None:
            countrate = self._counts / self.dt
            # If not in low-memory regime, cache the values
            if not self.low_memory or not self.input_counts:
                self._countrate = countrate

        return countrate

    @countrate.setter
    def countrate(self, value):
        value = np.asarray(value)
        if not value.shape == self.countrate.shape:
            raise ValueError('Can only assign new countrate array of the same '
                             'shape as the original array')
        self._countrate = value
        self._counts = None
        self._meancounts = None
        self._meancountrate = None
        self.input_counts = False

    @property
    def countrate_err(self):
        countrate_err = self._countrate_err
        if countrate_err is None and self._counts_err is not None:
            countrate_err = self._counts_err / self.dt
        elif countrate_err is None:
            countrate_err = 0

        # If not in low-memory regime, cache the values ONLY if they have
        # been changed!
        if countrate_err is not self._countrate_err:
            if not self.low_memory or not self.input_counts:
                self._countrate_err = countrate_err

        return countrate_err

    @countrate_err.setter
    def countrate_err(self, value):
        value = np.asarray(value)
        if not value.shape == self.countrate.shape:
            raise ValueError('Can only assign new error array of the same '
                             'shape as the original array')
        self._countrate_err = value
        self._counts_err = None

    @property
    def bin_lo(self):
        if self._bin_lo is None:
            self._bin_lo = self.time - 0.5 * self.dt
        return self._bin_lo

    @property
    def bin_hi(self):
        if self._bin_hi is None:
            self._bin_hi = self.time + 0.5 * self.dt
        return self._bin_hi

    def initial_optional_checks(self, time, counts, err):
        logging.info("Checking if light curve is well behaved. This "
                     "can take time, so if you are sure it is already "
                     "sorted, specify skip_checks=True at light curve "
                     "creation.")

        if not np.all(np.isfinite(time)):
            raise ValueError("There are inf or NaN values in "
                             "your time array!")

        if not np.all(np.isfinite(counts)):
            raise ValueError("There are inf or NaN values in "
                             "your counts array!")

        logging.info("Checking if light curve is sorted.")
        dt_array = np.diff(time)
        unsorted = np.any(dt_array < 0)

        if unsorted:
            logging.warning("The light curve is unsorted. Now, sorting...")
            order = np.argsort(time)
            time = time[order]
            counts = counts[order]
            if err is not None:
                err = err[order]
        return time, counts, err

    def check_lightcurve(self):
        """Make various checks on the lightcurve.

        It can be slow, use it if you are not sure about your
        input data.
        """
        # Issue a warning if the input time iterable isn't regularly spaced,
        # i.e. the bin sizes aren't equal throughout.

        check_gtis(self.gti)

        idxs = np.searchsorted(self.time, self.gti)
        uneven = False
        for idx in range(idxs.shape[0]):
            istart, istop = idxs[idx, 0], min(idxs[idx, 1], self.time.size - 1)

            local_diff = np.diff(self.time[istart:istop])
            if np.any(~np.isclose(local_diff, self.dt)):
                uneven = True

                break
        if uneven:
            simon("Bin sizes in input time array aren't equal throughout! "
                  "This could cause problems with Fourier transforms. "
                  "Please make the input time evenly sampled.")

    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:
            warnings.warn("MJDref is different in the two light curves")
            other = other.change_mjdref(self.mjdref)

        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, skip_checks=True,
                            dt=self.dt)

        return lc_new

    def __add__(self, other):
        """
        Add the counts of two light curves element by element, assuming the light curves
        have the same time array.

        This magic method adds two :class:`Lightcurve` objects having the same time
        array such that the corresponding counts arrays get summed up.

        GTIs are crossed, so that only common intervals are saved.

        Examples
        --------
        >>> time = [5, 10, 15]
        >>> count1 = [300, 100, 400]
        >>> count2 = [600, 1200, 800]
        >>> gti1 = [[0, 20]]
        >>> gti2 = [[0, 25]]
        >>> lc1 = Lightcurve(time, count1, gti=gti1, dt=5)
        >>> lc2 = Lightcurve(time, count2, gti=gti2, dt=5)
        >>> lc = lc1 + lc2
        >>> np.allclose(lc.counts, [ 900, 1300, 1200])
        True
        """

        return self._operation_with_other_lc(other, np.add)

    def __sub__(self, other):
        """
        Subtract the counts/flux of one light curve from the counts/flux of another
        light curve element by element, assuming the ``time`` arrays of the light curves
        match exactly.

        This magic method takes two :class:`Lightcurve` objects having the same
        ``time`` array and subtracts the ``counts`` of one :class:`Lightcurve` with
        that of another, while also updating ``countrate``, ``counts_err`` and ``countrate_err``
        correctly.

        GTIs are crossed, so that only common intervals are saved.

        Examples
        --------
        >>> time = [10, 20, 30]
        >>> count1 = [600, 1200, 800]
        >>> count2 = [300, 100, 400]
        >>> gti1 = [[0, 35]]
        >>> gti2 = [[5, 40]]
        >>> lc1 = Lightcurve(time, count1, gti=gti1, dt=10)
        >>> lc2 = Lightcurve(time, count2, gti=gti2, dt=10)
        >>> lc = lc1 - lc2
        >>> np.allclose(lc.counts, [ 300, 1100,  400])
        True
        """

        return self._operation_with_other_lc(other, np.subtract)

    def __neg__(self):
        """
        Implement the behavior of negation of the light curve objects.

        The negation operator ``-`` is supposed to invert the sign of the count
        values of a light curve object.

        Examples
        --------
        >>> time = [1, 2, 3]
        >>> count1 = [100, 200, 300]
        >>> count2 = [200, 300, 400]
        >>> lc1 = Lightcurve(time, count1)
        >>> lc2 = Lightcurve(time, count2)
        >>> lc_new = -lc1 + lc2
        >>> np.allclose(lc_new.counts, [100, 100, 100])
        True
        """
        lc_new = Lightcurve(self.time, -1 * self.counts,
                            err=self.counts_err, gti=self.gti,
                            mjdref=self.mjdref, skip_checks=True,
                            dt=self.dt)

        return lc_new

    def __len__(self):
        """
        Return the number of time bins of a light curve.

        This method implements overrides the ``len`` function for a :class:`Lightcurve`
        object and returns the length of the ``time`` array (which should be equal to the
        length of the ``counts`` and ``countrate`` arrays).

        Examples
        --------
        >>> time = [1, 2, 3]
        >>> count = [100, 200, 300]
        >>> lc = Lightcurve(time, count, dt=1)
        >>> len(lc)
        3
        """
        return self.n

    def __getitem__(self, index):
        """
        Return the corresponding count value at the index or a new :class:`Lightcurve`
        object upon slicing.

        This method adds functionality to retrieve the count value at
        a particular index. This also can be used for slicing and generating
        a new :class:`Lightcurve` object. GTIs are recalculated based on the new light
        curve segment

        If the slice object is of kind ``start:stop:step``, GTIs are also sliced,
        and rewritten as ``zip(time - self.dt /2, time + self.dt / 2)``

        Parameters
        ----------
        index : int or slice instance
            Index value of the time array or a slice object.

        Examples
        --------
        >>> time = [1, 2, 3, 4, 5, 6, 7, 8, 9]
        >>> count = [11, 22, 33, 44, 55, 66, 77, 88, 99]
        >>> lc = Lightcurve(time, count, dt=1)
        >>> np.isclose(lc[2], 33)
        True
        >>> np.allclose(lc[:2].counts, [11, 22])
        True
        """
        if isinstance(index, (int, np.integer)):
            return self.counts[index]
        elif isinstance(index, slice):
            start = assign_value_if_none(index.start, 0)
            stop = assign_value_if_none(index.stop, len(self.counts))
            step = assign_value_if_none(index.step, 1)

            new_counts = self.counts[start:stop:step]
            new_time = self.time[start:stop:step]

            new_gti = [[self.time[start] - 0.5 * self.dt,
                        self.time[stop - 1] + 0.5 * self.dt]]
            new_gti = np.asarray(new_gti)
            if step > 1:
                new_gt1 = np.array(list(zip(new_time - self.dt / 2,
                                            new_time + self.dt / 2)))
                new_gti = cross_two_gtis(new_gti, new_gt1)
            new_gti = cross_two_gtis(self.gti, new_gti)

            lc = Lightcurve(new_time, new_counts, mjdref=self.mjdref,
                            gti=new_gti, dt=self.dt, skip_checks=True,
                            err_dist=self.err_dist)
            if self._counts_err is not None:
                lc._counts_err = self._counts_err[start:stop:step]
            return lc
        else:
            raise IndexError("The index must be either an integer or a slice "
                             "object !")

    def __eq__(self, other_lc):
        """
        Compares two :class:`Lightcurve` objects.

        Light curves are equal only if their counts as well as times at which those counts occur equal.

        Examples
        --------
        >>> time = [1, 2, 3]
        >>> count1 = [100, 200, 300]
        >>> count2 = [100, 200, 300]
        >>> lc1 = Lightcurve(time, count1, dt=1)
        >>> lc2 = Lightcurve(time, count2, dt=1)
        >>> lc1 == lc2
        True
        """
        if not isinstance(other_lc, Lightcurve):
            raise ValueError('Lightcurve can only be compared with a Lightcurve Object')
        if (np.allclose(self.time, other_lc.time) and
                np.allclose(self.counts, other_lc.counts)):
            return True
        return False

    def baseline(self, lam, p, niter=10, offset_correction=False):
        """Calculate the baseline of the light curve, accounting for GTIs.

        Parameters
        ----------
        lam : float
            "smoothness" parameter. Larger values make the baseline stiffer
            Typically ``1e2 < lam < 1e9``
        p : float
            "asymmetry" parameter. Smaller values make the baseline more
            "horizontal". Typically ``0.001 < p < 0.1``, but not necessary.

        Other parameters
        ----------------
        offset_correction : bool, default False
            by default, this method does not align to the running mean of the
            light curve, but it goes below the light curve. Setting align to
            True, an additional step is done to shift the baseline so that it
            is shifted to the middle of the light curve noise distribution.


        Returns
        -------
        baseline : numpy.ndarray
            An array with the baseline of the light curve
        """
        baseline = np.zeros_like(self.time)
        for g in self.gti:
            good = create_gti_mask(self.time, [g], dt=self.dt)
            _, baseline[good] = \
                baseline_als(self.time[good], self.counts[good], lam, p,
                             niter, offset_correction=offset_correction,
                             return_baseline=True)

        return baseline

    @staticmethod
    def make_lightcurve(toa, dt, tseg=None, tstart=None, gti=None, mjdref=0,
                        use_hist=False):
        """
        Make a light curve out of photon arrival times, with a given time resolution ``dt``.
        Note that ``dt`` should be larger than the native time resolution of the instrument
        that has taken the data.

        Parameters
        ----------
        toa: iterable
            list of photon arrival times

        dt: float
            time resolution of the light curve (the bin width)

        tseg: float, optional, default ``None``
            The total duration of the light curve.
            If this is ``None``, then the total duration of the light curve will
            be the interval between the arrival between the first and the last
            photon in ``toa``.

                **Note**: If ``tseg`` is not divisible by ``dt`` (i.e. if ``tseg``/``dt`` is
                not an integer number), then the last fractional bin will be
                dropped!

        tstart: float, optional, default ``None``
            The start time of the light curve.
            If this is ``None``, the arrival time of the first photon will be used
            as the start time of the light curve.

        gti: 2-d float array
            ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``
            Good Time Intervals

        use_hist : bool
            Use ``np.histogram`` instead of ``np.bincounts``. Might be advantageous
            for very short datasets.

        Returns
        -------
        lc: :class:`Lightcurve` object
            A :class:`Lightcurve` object with the binned light curve
        """
        toa, mjdref = interpret_times(toa, mjdref=mjdref)

        toa = np.sort(np.asarray(toa))
        # tstart is an optional parameter to set a starting time for
        # the light curve in case this does not coincide with the first photon
        if tstart is None:
            # if tstart is not set, assume light curve starts with first photon
            tstart = toa[0]

        # compute the number of bins in the light curve
        # for cases where tseg/dt is not integer.
        # TODO: check that this is always consistent and that we
        # are not throwing away good events.

        if tseg is None:
            tseg = toa[-1] - tstart

        logging.info("make_lightcurve: tseg: " + str(tseg))

        timebin = int(tseg / dt)
        # If we are missing the next bin by just 1%, let's round up:
        if tseg / dt - timebin >= 0.99:
            timebin += 1

        logging.info("make_lightcurve: timebin:  " + str(timebin))

        tend = tstart + timebin * dt
        good = (tstart <= toa) & (toa < tend)
        if not use_hist:
            binned_toas = ((toa[good] - tstart) // dt).astype(np.int64)
            counts = \
                np.bincount(binned_toas, minlength=timebin)
            time = tstart + np.arange(0.5, 0.5 + len(counts)) * dt
        else:
            histbins = np.arange(tstart, tend + dt, dt)
            counts, histbins = np.histogram(toa[good], bins=histbins)
            time = histbins[:-1] + 0.5 * dt

        return Lightcurve(time, counts, gti=gti, mjdref=mjdref, dt=dt,
                          skip_checks=True, err_dist='poisson')

    def rebin(self, dt_new=None, f=None, method='sum'):
        """
        Rebin the light curve to a new time resolution. While the new
        resolution need not be an integer multiple of the previous time
        resolution, be aware that if it is not, the last bin will be cut
        off by the fraction left over by the integer division.

        Parameters
        ----------
        dt_new: float
            The new time resolution of the light curve. Must be larger than
            the time resolution of the old light curve!

        method: {``sum`` | ``mean`` | ``average``}, optional, default ``sum``
            This keyword argument sets whether the counts in the new bins
            should be summed or averaged.

        Other Parameters
        ----------------
        f: float
            the rebin factor. If specified, it substitutes ``dt_new`` with
            ``f*self.dt``

        Returns
        -------
        lc_new: :class:`Lightcurve` object
            The :class:`Lightcurve` object with the new, binned light curve.
        """

        if f is None and dt_new is None:
            raise ValueError('You need to specify at least one between f and '
                             'dt_new')
        elif f is not None:
            dt_new = f * self.dt

        if dt_new < self.dt:
            raise ValueError("New time resolution must be larger than "
                             "old time resolution!")

        bin_time, bin_counts, bin_err = [], [], []
        gti_new = []
        for g in self.gti:
            if g[1] - g[0] < dt_new:
                continue
            else:
                # find start and end of GTI segment in data
                start_ind = self.time.searchsorted(g[0])
                end_ind = self.time.searchsorted(g[1])

                t_temp = self.time[start_ind:end_ind]
                c_temp = self.counts[start_ind:end_ind]
                e_temp = self.counts_err[start_ind:end_ind]

                bin_t, bin_c, bin_e, _ = \
                    utils.rebin_data(t_temp, c_temp, dt_new,
                                     yerr=e_temp, method=method)

                bin_time.extend(bin_t)
                bin_counts.extend(bin_c)
                bin_err.extend(bin_e)
                gti_new.append(g)

        if len(gti_new) == 0:
            raise ValueError("No valid GTIs after rebin.")

        lc_new = Lightcurve(bin_time, bin_counts, err=bin_err,
                            mjdref=self.mjdref, dt=dt_new, gti=gti_new,
                            skip_checks=True)
        return lc_new

    def join(self, other, skip_checks=False):
        """
        Join two lightcurves into a single object.

        The new :class:`Lightcurve` object will contain time stamps from both the
        objects. The ``counts`` and ``countrate`` attributes in the resulting object
        will contain the union of the non-overlapping parts of the two individual objects,
        or the average in case of overlapping ``time`` arrays of both :class:`Lightcurve` objects.

        Good Time Intervals are also joined.

        Note : Ideally, the ``time`` array of both lightcurves should not overlap.

        Parameters
        ----------
        other : :class:`Lightcurve` object
            The other :class:`Lightcurve` object which is supposed to be joined with.
        skip_checks: bool
            If True, the user specifies that data are already sorted and
            contain no infinite or nan points. Use at your own risk.

        Returns
        -------
        lc_new : :class:`Lightcurve` object
            The resulting :class:`Lightcurve` object.

        Examples
        --------
        >>> time1 = [5, 10, 15]
        >>> count1 = [300, 100, 400]
        >>> time2 = [20, 25, 30]
        >>> count2 = [600, 1200, 800]
        >>> lc1 = Lightcurve(time1, count1, dt=5)
        >>> lc2 = Lightcurve(time2, count2, dt=5)
        >>> lc = lc1.join(lc2)
        >>> lc.time
        array([ 5, 10, 15, 20, 25, 30])
        >>> np.allclose(lc.counts, [ 300,  100,  400,  600, 1200,  800])
        True
        """
        if self.mjdref != other.mjdref:
            warnings.warn("MJDref is different in the two light curves")
            other = other.change_mjdref(self.mjdref)

        if self.dt != other.dt:
            utils.simon("The two light curves have different bin widths.")

        if(self.tstart < other.tstart):
            first_lc = self
            second_lc = other
        else:
            first_lc = other
            second_lc = self

        if len(np.intersect1d(self.time, other.time) > 0):

            utils.simon("The two light curves have overlapping time ranges. "
                        "In the common time range, the resulting count will "
                        "be the average of the counts in the two light "
                        "curves. If you wish to sum, use `lc_sum = lc1 + "
                        "lc2`.")
            valid_err = False

            if self.err_dist.lower() != other.err_dist.lower():
                simon("Lightcurves have different statistics!"
                      "We are setting the errors to zero.")

            elif self.err_dist.lower() in valid_statistics:
                valid_err = True
            # More conditions can be implemented for other statistics
            else:
                raise StingrayError("Statistics not recognized."
                                    " Please use one of these: "
                                    "{}".format(valid_statistics))

            from collections import Counter
            counts = Counter()
            counts_err = Counter()

            for i, time in enumerate(first_lc.time):
                counts[time] = first_lc.counts[i]
                counts_err[time] = first_lc.counts_err[i]

            for i, time in enumerate(second_lc.time):

                if counts.get(time) is not None:  # Common time
                    counts[time] = (counts[time] + second_lc.counts[i]) / 2
                    counts_err[time] = \
                        np.sqrt(((counts_err[time] ** 2) +
                                 (second_lc.counts_err[i] ** 2)) / 2)

                else:
                    counts[time] = second_lc.counts[i]
                    counts_err[time] = second_lc.counts_err[i]

            new_time = list(counts.keys())
            new_counts = list(counts.values())
            if(valid_err):
                new_counts_err = list(counts_err.values())
            else:
                new_counts_err = np.zeros_like(new_counts)

            del[counts, counts_err]

        else:

            new_time = np.concatenate([first_lc.time, second_lc.time])
            new_counts = np.concatenate([first_lc.counts, second_lc.counts])
            new_counts_err = \
                np.concatenate([first_lc.counts_err, second_lc.counts_err])

        new_time = np.asarray(new_time)
        new_counts = np.asarray(new_counts)
        new_counts_err = np.asarray(new_counts_err)
        gti = join_gtis(self.gti, other.gti)

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

        return lc_new

    def truncate(self, start=0, stop=None, method="index"):
        """
        Truncate a :class:`Lightcurve` object.

        This method takes a ``start`` and a ``stop`` point (either as indices,
        or as times in the same unit as those in the ``time`` attribute, and truncates
        all bins before ``start`` and after ``stop``, then returns a new :class:`Lightcurve`
        object with the truncated light curve.

        Parameters
        ----------
        start : int, default 0
            Index (or time stamp) of the starting point of the truncation. If no value is set
            for the start point, then all points from the first element in the ``time`` array
            are taken into account.

        stop : int, default ``None``
            Index (or time stamp) of the ending point (exclusive) of the truncation. If no
            value of stop is set, then points including the last point in
            the counts array are taken in count.

        method : {``index`` | ``time``}, optional, default ``index``
            Type of the start and stop values. If set to ``index`` then
            the values are treated as indices of the counts array, or
            if set to ``time``, the values are treated as actual time values.

        Returns
        -------
        lc_new: :class:`Lightcurve` object
            The :class:`Lightcurve` object with truncated time and counts
            arrays.

        Examples
        --------
        >>> time = [1, 2, 3, 4, 5, 6, 7, 8, 9]
        >>> count = [10, 20, 30, 40, 50, 60, 70, 80, 90]
        >>> lc = Lightcurve(time, count, dt=1)
        >>> lc_new = lc.truncate(start=2, stop=8)
        >>> np.allclose(lc_new.counts, [30, 40, 50, 60, 70, 80])
        True
        >>> lc_new.time
        array([3, 4, 5, 6, 7, 8])
        >>> # Truncation can also be done by time values
        >>> lc_new = lc.truncate(start=6, method='time')
        >>> lc_new.time
        array([6, 7, 8, 9])
        >>> np.allclose(lc_new.counts, [60, 70, 80, 90])
        True
        """

        if not isinstance(method, str):
            raise TypeError("method key word argument is not "
                            "a string !")

        if method.lower() not in ['index', 'time']:
            raise ValueError("Unknown method type " + method + ".")

        if method.lower() == 'index':
            return self._truncate_by_index(start, stop)
        else:
            return self._truncate_by_time(start, stop)

    def _truncate_by_index(self, start, stop):
        """Private method for truncation using index values."""
        time_new = self.time[start:stop]
        counts_new = self.counts[start:stop]
        counts_err_new = self.counts_err[start:stop]
        gti = \
            cross_two_gtis(self.gti,
                           np.asarray([[self.time[start] - 0.5 * self.dt,
                                        time_new[-1] + 0.5 * self.dt]]))

        return Lightcurve(time_new, counts_new, err=counts_err_new, gti=gti,
                          dt=self.dt, skip_checks=True)

    def _truncate_by_time(self, start, stop):
        """Helper method for truncation using time values.

        Parameters
        ----------
        start : float
            start time for new light curve; all time bins before this time will be discarded

        stop : float
            stop time for new light curve; all time bins after this point will be discarded

        Returns
        -------
            new_lc : Lightcurve
                A new :class:`Lightcurve` object with the truncated time bins

        """

        if stop is not None:
            if start > stop:
                raise ValueError("start time must be less than stop time!")

        if not start == 0:
            start = self.time.searchsorted(start)

        if stop is not None:
            stop = self.time.searchsorted(stop)

        return self._truncate_by_index(start, stop)

    def split(self, min_gap, min_points=1):
        """
        For data with gaps, it can sometimes be useful to be able to split
        the light curve into separate, evenly sampled objects along those
        data gaps. This method allows to do this: it finds data gaps of a
        specified minimum size, and produces a list of new `Lightcurve`
        objects for each contiguous segment.

        Parameters
        ----------
        min_gap : float
            The length of a data gap, in the same units as the `time` attribute
            of the `Lightcurve` object. Any smaller gaps will be ignored, any
            larger gaps will be identified and used to split the light curve.

        min_points : int, default 1
            The minimum number of data points in each light curve. Light
            curves with fewer data points will be ignored.

        Returns
        -------
        lc_split : iterable of `Lightcurve` objects
            The list of all contiguous light curves

        Examples
        --------
        >>> time = np.array([1, 2, 3, 6, 7, 8, 11, 12, 13])
        >>> counts = np.random.rand(time.shape[0])
        >>> lc = Lightcurve(time, counts, dt=1, skip_checks=True)
        >>> split_lc = lc.split(1.5)

        """

        # calculate the difference between time bins
        tdiff = np.diff(self.time)
        # find all distances between time bins that are larger than `min_gap`
        gap_idx = np.where(tdiff >= min_gap)[0]

        # tolerance for the newly created GTIs: Note that this seems to work
        # with a tolerance of 2, but not if I substitute 10. I don't know why
        epsilon = np.min(tdiff) / 2.0

        # calculate new GTIs
        gti_start = np.hstack([self.time[0] - epsilon, self.time[gap_idx + 1] - epsilon])
        gti_stop = np.hstack([self.time[gap_idx] + epsilon, self.time[-1] + epsilon])

        gti = np.vstack([gti_start, gti_stop]).T
        if hasattr(self, 'gti') and self.gti is not None:
            gti = cross_two_gtis(self.gti, gti)

        lc_split = self.split_by_gti(gti, min_points=min_points)
        return lc_split

    def sort(self, reverse=False):
        """
        Sort a Lightcurve object by time.

        A Lightcurve can be sorted in either increasing or decreasing order
        using this method. The time array gets sorted and the counts array is
        changed accordingly.

        Parameters
        ----------
        reverse : boolean, default False
            If True then the object is sorted in reverse order.

        Examples
        --------
        >>> time = [2, 1, 3]
        >>> count = [200, 100, 300]
        >>> lc = Lightcurve(time, count, dt=1)
        >>> lc_new = lc.sort()
        >>> lc_new.time
        array([1, 2, 3])
        >>> np.allclose(lc_new.counts, [100, 200, 300])
        True

        Returns
        -------
        lc_new: :class:`Lightcurve` object
            The :class:`Lightcurve` object with sorted time and counts
            arrays.
        """
        new_time, new_counts, new_counts_err = \
            zip(*sorted(zip(self.time, self.counts, self.counts_err),
                        reverse=reverse))
        new_time = np.asarray(new_time)
        new_counts = np.asarray(new_counts)
        new_counts_err = np.asarray(new_counts_err)

        new_lc = Lightcurve(new_time, new_counts, err=new_counts_err,
                            gti=self.gti, dt=self.dt, mjdref=self.mjdref,
                            skip_checks=True)

        return new_lc

    def sort_counts(self, reverse=False):
        """
        Sort a :class:`Lightcurve` object in accordance with its counts array.

        A :class:`Lightcurve` can be sorted in either increasing or decreasing order
        using this method. The counts array gets sorted and the time array is
        changed accordingly.

        Parameters
        ----------
        reverse : boolean, default ``False``
            If ``True`` then the object is sorted in reverse order.

        Returns
        -------
        lc_new: :class:`Lightcurve` object
            The :class:`Lightcurve` object with sorted ``time`` and ``counts``
            arrays.

        Examples
        --------
        >>> time = [1, 2, 3]
        >>> count = [200, 100, 300]
        >>> lc = Lightcurve(time, count, dt=1)
        >>> lc_new = lc.sort_counts()
        >>> lc_new.time
        array([2, 1, 3])
        >>> np.allclose(lc_new.counts, [100, 200, 300])
        True
        """

        new_counts, new_time, new_counts_err = \
            zip(*sorted(zip(self.counts, self.time, self.counts_err),
                        reverse=reverse))

        new_lc = Lightcurve(new_time, new_counts, err=new_counts_err,
                            gti=self.gti, dt=self.dt, mjdref=self.mjdref,
                            skip_checks=True)

        return new_lc

    def estimate_chunk_length(self, *args, **kwargs):
        """Deprecated alias of estimate_segment_size.
        """
        warnings.warn("This function was renamed to estimate_segment_size", DeprecationWarning)
        return self.estimate_segment_size(*args, **kwargs)

    def estimate_segment_size(self, min_total_counts=100, min_time_bins=100):
        """Estimate a reasonable segment length for chunk-by-chunk analysis.

        Choose a reasonable length for time segments, given a minimum number of total
        counts in the segment, and a minimum number of time bins in the segment.

        The user specifies a condition on the total counts in each segment and
        the minimum number of time bins.

        Other Parameters
        ----------------
        min_total_counts : int
            Minimum number of counts for each chunk
        min_time_bins : int
            Minimum number of time bins

        Returns
        -------
        segment_size : float
            The length of the light curve chunks that satisfies the conditions

        Examples
        --------
        >>> import numpy as np
        >>> time = np.arange(150)
        >>> count = np.zeros_like(time) + 3
        >>> lc = Lightcurve(time, count, dt=1)
        >>> lc.estimate_segment_size(min_total_counts=10, min_time_bins=3)
        4.0
        >>> lc.estimate_segment_size(min_total_counts=10, min_time_bins=5)
        5.0
        >>> count[2:4] = 1
        >>> lc = Lightcurve(time, count, dt=1)
        >>> lc.estimate_segment_size(min_total_counts=3, min_time_bins=1)
        3.0
        >>> # A slightly more complex example
        >>> dt=0.2
        >>> time = np.arange(0, 1000, dt)
        >>> counts = np.random.poisson(100, size=len(time))
        >>> lc = Lightcurve(time, counts, dt=dt)
        >>> lc.estimate_segment_size(100, 2)
        0.4
        >>> min_total_bins = 40
        >>> lc.estimate_segment_size(100, 40)
        8.0
        """

        rough_estimate = np.ceil(min_total_counts / self.meancounts) * self.dt

        segment_size = np.max([rough_estimate, min_time_bins * self.dt])

        keep_searching = True
        while keep_searching:
            start_times, stop_times, results = \
                self.analyze_lc_chunks(segment_size, np.sum)
            mincounts = np.min(results)
            if mincounts >= min_total_counts:
                keep_searching = False
            else:
                segment_size += self.dt

        return segment_size

    def analyze_lc_chunks(self, segment_size, func, fraction_step=1, **kwargs):
        """Analyze segments of the light curve with any function.

        Parameters
        ----------
        segment_size : float
            Length in seconds of the light curve segments
        func : function
            Function accepting a :class:`Lightcurve` object as single argument, plus
            possible additional keyword arguments, and returning a number or a
            tuple - e.g., ``(result, error)`` where both ``result`` and ``error`` are
            numbers.

        Other parameters
        ----------------
        fraction_step : float
            If the step is not a full ``segment_size`` but less (e.g. a moving window),
            this indicates the ratio between step step and ``segment_size`` (e.g.
            0.5 means that the window shifts of half ``segment_size``)
        kwargs : keyword arguments
            These additional keyword arguments, if present, they will be passed
            to ``func``

        Returns
        -------
        start_times : array
            Lower time boundaries of all time segments.
        stop_times : array
            upper time boundaries of all segments.
        result : array of N elements
            The result of ``func`` for each segment of the light curve

        Examples
        --------
        >>> import numpy as np
        >>> time = np.arange(0, 10, 0.1)
        >>> counts = np.zeros_like(time) + 10
        >>> lc = Lightcurve(time, counts, dt=0.1)
        >>> # Define a function that calculates the mean
        >>> mean_func = lambda x: np.mean(x)
        >>> # Calculate the mean in segments of 5 seconds
        >>> start, stop, res = lc.analyze_lc_chunks(5, mean_func)
        >>> len(res) == 2
        True
        >>> np.allclose(res, 10)
        True
        """
        start, stop = bin_intervals_from_gtis(self.gti, segment_size,
                                              self.time,
                                              fraction_step=fraction_step,
                                              dt=self.dt)
        start_times = self.time[start] - 0.5 * self.dt

        # Remember that stop is one element above the last element, because
        # it's defined to be used in intervals start:stop
        stop_times = self.time[stop - 1] + self.dt * 1.5

        results = []
        for i, (st, sp) in enumerate(zip(start, stop)):
            lc_filt = self[st:sp]
            res = func(lc_filt, **kwargs)
            results.append(res)

        results = np.array(results)

        if len(results.shape) == 2:
            results = [results[:, i] for i in range(results.shape[1])]
        return start_times, stop_times, results

    def to_lightkurve(self):
        """
        Returns a `lightkurve.LightCurve` object.
        This feature requires `Lightkurve
        <https://docs.lightkurve.org/index.html/>`_ to be installed
        (e.g. ``pip install lightkurve``).  An `ImportError` will
        be raised if this package is not available.

        Returns
        -------
        lightcurve : `lightkurve.LightCurve`
            A lightkurve LightCurve object.
        """
        try:
            from lightkurve import LightCurve as lk
        except ImportError:
            raise ImportError("You need to install Lightkurve to use "
                              "the Lightcurve.to_lightkurve() method.")
        time = Time(self.time / 86400 + self.mjdref, format="mjd", scale="utc")
        return lk(time=time, flux=self.counts, flux_err=self.counts_err)

    @staticmethod
    def from_lightkurve(lk, skip_checks=True):
        """
        Creates a new `Lightcurve` from a `lightkurve.LightCurve`.

        Parameters
        ----------
        lk : `lightkurve.LightCurve`
            A lightkurve LightCurve object
        skip_checks: bool
            If True, the user specifies that data are already sorted and contain no
            infinite or nan points. Use at your own risk.
        """

        return Lightcurve(time=lk.time, counts=lk.flux,
                          err=lk.flux_err, input_counts=False,
                          skip_checks=skip_checks)

    def to_astropy_timeseries(self):
        return self._to_astropy_object(kind="timeseries")

    def to_astropy_table(self):
        return self._to_astropy_object(kind="table")

    def _to_astropy_object(self, kind="table"):
        data = {}

        for attr in ['_counts', '_counts_err', '_countrate', '_countrate_err',
                     '_bin_lo', '_bin_hi']:
            if hasattr(self, attr) and getattr(self, attr) is not None:
                data[attr.lstrip('_')] = np.asarray(getattr(self, attr))

        if kind.lower() == "table":
            data['time'] = self.time
            ts = Table(data)
        elif kind.lower() == "timeseries":
            from astropy.timeseries import TimeSeries
            ts = TimeSeries(data=data, time=TimeDelta(self.time * u.s))
        else:  # pragma: no cover
            raise ValueError("Invalid kind (accepted: table or timeseries)")

        for attr in ['_gti', 'mjdref', '_meancounts', '_meancountrate',
                     'instr', 'mission', 'dt', 'err_dist']:
            if hasattr(self, attr) and getattr(self, attr) is not None:
                ts.meta[attr.lstrip('_')] = getattr(self, attr)

        return ts

    @staticmethod
    def from_astropy_timeseries(ts, **kwargs):
        return Lightcurve._from_astropy_object(ts, **kwargs)

    @staticmethod
    def from_astropy_table(ts, **kwargs):
        return Lightcurve._from_astropy_object(ts, **kwargs)

    @staticmethod
    def _from_astropy_object(
            ts, err_dist='poisson', skip_checks=True):

        if hasattr(ts, 'time'):
            time = ts.time
        else:
            time = ts['time']

        kwargs = ts.meta
        err = None
        input_counts = True

        if "counts_err" in ts.colnames:
            err = ts["counts_err"]
        elif "countrate_err" in ts.colnames:
            err = ts["countrate_err"]

        if 'counts' in ts.colnames:
            counts = ts['counts']
        elif 'countrate' in ts.colnames:
            counts = ts['countrate']
            input_counts = False
        else:
            raise ValueError('Input timeseries must contain at least a '
                             '`counts` or a `countrate` column')

        kwargs.update({'time': time, 'counts': counts, 'err': err,
                       'input_counts': input_counts,
                       'skip_checks': skip_checks})
        if 'err_dist' not in kwargs:
            kwargs['err_dist'] = err_dist

        lc = Lightcurve(**kwargs)

        return lc

    def plot(self, witherrors=False, labels=None, axis=None, title=None,
             marker='-', save=False, filename=None):
        """
        Plot the light curve using ``matplotlib``.

        Plot the light curve object on a graph ``self.time`` on x-axis and
        ``self.counts`` on y-axis with ``self.counts_err`` optionally
        as error bars.

        Parameters
        ----------
        witherrors: boolean, default False
            Whether to plot the Lightcurve with errorbars or not

        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 the ``matplotlib`` figure. For example
            it can be a list like ``[xmin, xmax, ymin, ymax]`` or any other
            acceptable argument for the``matplotlib.pyplot.axis()`` method.

        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``.
        """

        fig = plt.figure()
        if witherrors:
            fig = plt.errorbar(self.time, self.counts, yerr=self.counts_err,
                               fmt=marker)
        else:
            fig = plt.plot(self.time, self.counts, 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('out.png')
            else:
                plt.savefig(filename)
        else:
            plt.show(block=False)

    @classmethod
    def read(cls, filename, fmt=None, format_=None, err_dist='gauss',
             skip_checks=False):
        """
        Read a :class:`Lightcurve` object from file.

        Currently supported formats are

        * pickle (not recommended for long-term storage)
        * hea : FITS Light curves from HEASARC-supported missions.
        * any other formats compatible with the writers in
          :class:`astropy.table.Table` (ascii.ecsv, hdf5, etc.)

        Files that need the :class:`astropy.table.Table` interface MUST contain
        at least a ``time`` column and a ``counts`` or ``countrate`` column.
        The default ascii format is enhanced CSV (ECSV). Data formats
        supporting the serialization of metadata (such as ECSV and HDF5) can
        contain all lightcurve attributes such as ``dt``, ``gti``, etc with
        no significant loss of information. Other file formats might lose part
        of the metadata, so must be used with care.

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

        fmt: str
            Available options are 'pickle', 'hea', and any `Table`-supported
            format such as 'hdf5', 'ascii.ecsv', etc.

        Other parameters
        ----------------

        err_dist: str, default='gauss'
            Default error distribution if not specified in the file (e.g. for
            ASCII files). The default is 'gauss' just because it is likely
            that people using ASCII light curves will want to specify Gaussian
            error bars, if any.
        skip_checks : bool
            See :class:`Lightcurve` documentation

        Returns
        --------
        lc : :class:`Lightcurve` object
        """
        if fmt is None and format_ is not None:
            warnings.warn(
                "The format_ keyword for read and write is deprecated. "
                "Use fmt instead", DeprecationWarning)
            fmt = format_

        if fmt.lower() in ('hea', 'ogip'):
            data = lcurve_from_fits(filename)
            data.update({'err_dist': err_dist, 'skip_checks': skip_checks})
            return Lightcurve(**data)

        return super().read(filename=filename, fmt=fmt)

    def split_by_gti(self, gti=None, min_points=2):
        """
        Split the current :class:`Lightcurve` object into a list of :class:`Lightcurve` objects, one
        for each continuous GTI segment as defined in the ``gti`` attribute.

        Parameters
        ----------
        min_points : int, default 1
            The minimum number of data points in each light curve. Light
            curves with fewer data points will be ignored.

        Returns
        -------
        list_of_lcs : list
            A list of :class:`Lightcurve` objects, one for each GTI segment
        """

        if gti is None:
            gti = self.gti

        list_of_lcs = []

        start_bins, stop_bins = gti_border_bins(gti, self.time, self.dt)
        for i in range(len(start_bins)):
            start = start_bins[i]
            stop = stop_bins[i]

            if np.isclose(stop - start, 1):
                logging.warning("Segment with a single time bin! Ignoring this segment!")
                continue
            if (stop - start) < min_points:
                continue

            # Note: GTIs are consistent with default in this case!
            new_lc = Lightcurve(self.time[start:stop], self.counts[start:stop],
                                err=self.counts_err[start:stop],
                                mjdref=self.mjdref, gti=[gti[i]],
                                dt=self.dt, err_dist=self.err_dist,
                                skip_checks=True)
            list_of_lcs.append(new_lc)

        return list_of_lcs

    def apply_gtis(self):
        """
        Apply GTIs to a light curve. Filters the ``time``, ``counts``,
        ``countrate``, ``counts_err`` and ``countrate_err`` arrays for all bins
        that fall into Good Time Intervals and recalculates mean countrate
        and the number of bins.
        """
        check_gtis(self.gti)

        good = self.mask

        # nota bene: We set the private properties, otherwise we'll get a
        # ValueError from changing the shape of the arrays.
        self._time = self.time[good]
        self._counts = self.counts[good]
        if self._counts_err is not None:
            self._counts_err = self._counts_err[good]
        self._countrate = None
        self._countrate_err = None
        self._mask = None

        self._meanrate = None
        self._meancounts = None
        self._n = None
        self.tseg = np.max(self.gti) - np.min(self.gti)
        self.tstart = self.time - 0.5 * self.dt def show_progress(a):
    return a def avg_cs_from_iterables(
    flux_iterable1,
    flux_iterable2,
    dt,
    norm="frac",
    use_common_mean=True,
    silent=False,
    fullspec=False,
    power_type="all",
    return_auxil=False
):
    """Calculate the average cross spectrum from an iterable of light curves

    Parameters
    ----------
    flux_iterable1 : `iterable` of `np.array`s or of tuples (`np.array`, `np.array`)
        Iterable providing either equal-length series of count measurements, or
        of tuples (fluxes, errors). They must all be of the same length.
    flux_iterable2 : `iterable` of `np.array`s or of tuples (`np.array`, `np.array`)
        Same as ``flux_iterable1``, for the reference channel
    dt : float
        Time resolution of the light curves used to produce periodograms

    Other Parameters
    ----------------
    norm : str, default "frac"
        The normalization of the periodogram. "abs" is absolute rms, "frac" is
        fractional rms, "leahy" is Leahy+83 normalization, and "none" is the
        unnormalized periodogram
    use_common_mean : bool, default True
        The mean of the light curve can be estimated in each interval, or on
        the full light curve. This gives different results (Alston+2013).
        Here we assume the mean is calculated on the full light curve, but
        the user can set ``use_common_mean`` to False to calculate it on a
        per-segment basis.
    fullspec : bool, default False
        Return the full periodogram, including negative frequencies
    silent : bool, default False
        Silence the progress bars
    power_type : str, default 'all'
        If 'all', give complex powers. If 'abs', the absolute value; if 'real',
        the real part
    return_auxil : bool, default False
        Return the auxiliary unnormalized PDSs from the two separate channels

    Returns
    -------
    results : :class:`astropy.table.Table`
        Table containing the following columns:
        freq : `np.array`
            The frequencies
        power : `np.array`
            The normalized cross spectral powers
        unnorm_power : `np.array`
            The unnormalized cross spectral power
        unnorm_pds1 : `np.array`
            The unnormalized auxiliary PDS from channel 1. Only returned if ``return_auxil`` is ``True``
        unnorm_pds2 : `np.array`
            The unnormalized auxiliary PDS from channel 2. Only returned if ``return_auxil`` is ``True``

        And a number of other useful diagnostics in the metadata, including
        all attributes needed to allocate Crossspectrum objects, such as all
        the input arguments of this function (``dt``, ``segment_size``), and,
        e.g.
        n : int
            the number of bins in the light curves used in each segment
        m : int
            the number of averaged periodograms
        mean : float
            the mean flux (geometrical average of the mean fluxes in the two
            channels)

    """

    local_show_progress = show_progress
    if silent:
        def local_show_progress(a):
            return a

    # Initialize stuff
    cross = unnorm_cross = unnorm_pds1 = unnorm_pds2 = pds1 = pds2 = None
    n_ave = 0

    sum_of_photons1 = sum_of_photons2 = 0
    common_variance1 = common_variance2 = common_variance = None

    for flux1, flux2 in local_show_progress(zip(flux_iterable1, flux_iterable2)):
        if flux1 is None or flux2 is None or \
                np.all(flux1 == 0) or np.all(flux2 == 0):
            continue

        # Does the flux iterable return the uncertainty?
        # If so, define the variances
        variance1 = variance2 = None
        if isinstance(flux1, tuple):
            flux1, err1 = flux1
            variance1 = np.mean(err1) ** 2
        if isinstance(flux2, tuple):
            flux2, err2 = flux2
            variance2 = np.mean(err2) ** 2

        # Only use the variance if both flux iterables define it.
        if variance1 is None or variance2 is None:
            variance1 = variance2 = None
        else:
            common_variance1 = sum_if_not_none_or_initialize(common_variance1, variance1)
            common_variance2 = sum_if_not_none_or_initialize(common_variance2, variance2)

        n_bin = flux1.size

        # At the first loop, we define the frequency array and the range of positive
        # frequency bins (after the first loop, cross will not be None nymore)
        if cross is None:
            freq = fftfreq(n_bin, dt)
            fgt0 = positive_fft_bins(n_bin)

        # Calculate the FFTs
        ft1 = fft(flux1)
        ft2 = fft(flux2)

        # Calculate the sum of each light curve, to calculate the mean
        n_ph1 = flux1.sum()
        n_ph2 = flux2.sum()
        n_ph = np.sqrt(n_ph1 * n_ph2)

        # Calculate the unnormalized cross spectrum
        unnorm_power = ft1 * ft2.conj()
        unnorm_pd1 = unnorm_pd2 = 0

        # If requested, calculate the auxiliary PDSs
        if return_auxil:
            unnorm_pd1 = (ft1 * ft1.conj()).real
            unnorm_pd2 = (ft2 * ft2.conj()).real

        # Accumulate the sum to calculate the total mean of the lc
        sum_of_photons1 += n_ph1
        sum_of_photons2 += n_ph2

        # Take only positive frequencies unless the user wants the full spectrum
        if not fullspec:
            unnorm_power = unnorm_power[fgt0]
            if return_auxil:
                unnorm_pd1 = unnorm_pd1[fgt0]
                unnorm_pd2 = unnorm_pd2[fgt0]

        cs_seg = unnorm_power
        p1_seg = unnorm_pd1
        p2_seg = unnorm_pd2

        # If normalization has to be done interval by interval, do it here.
        if not use_common_mean:
            mean1 = n_ph1 / n_bin
            mean2 = n_ph2 / n_bin
            mean = n_ph / n_bin
            variance = None

            if variance1 is not None:
                variance = np.sqrt(variance1 * variance2)

            cs_seg = normalize_periodograms(
                unnorm_power, dt, n_bin, mean, n_ph=n_ph, norm=norm, power_type=power_type, variance=variance
            )
            p1_seg = normalize_periodograms(
                unnorm_pd1, dt, n_bin, mean1, n_ph=n_ph1, norm=norm, power_type=power_type, variance=variance1
            )
            p2_seg = normalize_periodograms(
                unnorm_pd2, dt, n_bin, mean2, n_ph=n_ph2, norm=norm, power_type=power_type, variance=variance2
            )

        # Initialize or accumulate final averaged spectra
        cross = sum_if_not_none_or_initialize(cross, cs_seg)
        unnorm_cross = sum_if_not_none_or_initialize(unnorm_cross, unnorm_power)

        if return_auxil:
            unnorm_pds1 = sum_if_not_none_or_initialize(unnorm_pds1, unnorm_pd1)
            unnorm_pds2 = sum_if_not_none_or_initialize(unnorm_pds2, unnorm_pd2)
            pds1 = sum_if_not_none_or_initialize(pds1, p1_seg)
            pds2 = sum_if_not_none_or_initialize(pds2, p2_seg)

        n_ave += 1

    # If no valid intervals were found, return only `None`s
    if cross is None:
        return None

    # Calculate the mean number of photons per chunk
    n_ph1 = sum_of_photons1 / n_ave
    n_ph2 = sum_of_photons2 / n_ave
    n_ph = np.sqrt(n_ph1 * n_ph2)

    # Calculate the common mean number of photons per bin
    common_mean1 = n_ph1 / n_bin
    common_mean2 = n_ph2 / n_bin
    common_mean = n_ph / n_bin

    if common_variance1 is not None:
        # Note: the variances we summed were means, not sums. Hence M, not M*N
        common_variance1 /= n_ave
        common_variance2 /= n_ave
        common_variance = np.sqrt(common_variance1 * common_variance2)

    # Transform the sums into averages
    cross /= n_ave
    unnorm_cross /= n_ave
    if return_auxil:
        unnorm_pds1 /= n_ave
        unnorm_pds2 /= n_ave

    # Finally, normalize the cross spectrum (only if not already done on an
    # interval-to-interval basis)
    if use_common_mean:
        cross = normalize_periodograms(
            unnorm_cross,
            dt,
            n_bin,
            common_mean,
            n_ph=n_ph,
            norm=norm,
            variance=common_variance,
            power_type=power_type,
        )
        if return_auxil:
            pds1 = normalize_periodograms(
                unnorm_pds1,
                dt,
                n_bin,
                common_mean1,
                n_ph=n_ph1,
                norm=norm,
                variance=common_variance1,
                power_type=power_type,
            )
            pds2 = normalize_periodograms(
                unnorm_pds2,
                dt,
                n_bin,
                common_mean2,
                n_ph=n_ph2,
                norm=norm,
                variance=common_variance2,
                power_type=power_type,
            )
    # If the user does not want negative frequencies, don't give them
    if not fullspec:
        freq = freq[fgt0]

    results = Table()
    results["freq"] = freq
    results["power"] = cross
    results["unnorm_power"] = unnorm_cross
    results.meta.update({"n": n_bin, "m": n_ave, "dt": dt,
                         "norm": norm,
                         "df": 1 / (dt * n_bin),
                         "segment_size": dt * n_bin,
                         "nphots": n_ph,
                         "nphots1": n_ph1, "nphots2": n_ph2,
                         "countrate1": common_mean1 / dt,
                         "countrate2": common_mean2 / dt,
                         "mean": common_mean,
                         "mean1": common_mean1,
                         "mean2": common_mean2,
                         "power_type": power_type,
                         "fullspec": fullspec,
                         "variance": common_variance,
                         "variance1": common_variance1,
                         "variance2": common_variance2})

    if return_auxil:
        results["pds1"] = pds1
        results["pds2"] = pds2
        results["unnorm_pds1"] = unnorm_pds1
        results["unnorm_pds2"] = unnorm_pds2

    return results def error_on_averaged_cross_spectrum(cross_power, seg_power, ref_power, n_ave, seg_power_noise, ref_power_noise, common_ref="False"):
    """Error on cross spectral quantities, From Ingram 2019.

    Note: this is only valid for a very large number of averaged powers.
    Beware if n_ave < 50 or so.

    Parameters
    ----------
    cross_power : complex `np.array`
        cross spectrum
    seg_power : float `np.array`
        sub-band periodogram
    ref_power : float `np.array`
        reference-band periodogram
    seg_power_noise : float
        Poisson noise level of the sub-band periodogram
    ref_power_noise : float
        Poisson noise level of the reference-band periodogram
    n_ave : int
        number of intervals that have been averaged to obtain the input spectra

    Other Parameters
    ----------------
    common_ref : bool, default False
        Are data in the sub-band also included in the reference band?

    Returns
    -------
    dRe : float `np.array`
        Error on the real part of the cross spectrum
    dIm : float `np.array`
        Error on the imaginary part of the cross spectrum
    dphi : float `np.array`
        Error on the angle (or phase lag)
    dG : float `np.array`
        Error on the modulus of the cross spectrum

    """
    if n_ave < 30:
        warnings.warn("n_ave is below 30. Please note that the error bars "
                      "on the quantities derived from the cross spectrum "
                      "are only reliable for a large number of averaged powers.")
    two_n_ave = 2 * n_ave
    if common_ref:
        Gsq = (cross_power * np.conj(cross_power)).real
        bsq = bias_term(seg_power, ref_power, seg_power_noise, ref_power_noise, n_ave)
        frac = (Gsq - bsq) / (ref_power - ref_power_noise)
        power_over_2n = ref_power / two_n_ave

        # Eq. 18
        dRe = dIm = dG = np.sqrt(power_over_2n * (seg_power - frac))
        # Eq. 19
        dphi = np.sqrt(power_over_2n * (seg_power / (Gsq - bsq) -
                       1 / (ref_power - ref_power_noise)))

    else:
        PrPs = ref_power * seg_power
        dRe = np.sqrt((PrPs + cross_power.real ** 2 - cross_power.imag ** 2) / two_n_ave)
        dIm = np.sqrt((PrPs - cross_power.real ** 2 + cross_power.imag ** 2) / two_n_ave)
        gsq = raw_coherence(cross_power, seg_power, ref_power,
                            seg_power_noise, ref_power_noise, n_ave)
        dphi = np.sqrt((1 - gsq) / (2 * gsq * n_ave))
        dG = np.sqrt(PrPs / n_ave)

    return dRe, dIm, dphi, dG def avg_cs_from_events(
    times1,
    times2,
    gti,
    segment_size,
    dt,
    norm="frac",
    use_common_mean=True,
    fullspec=False,
    silent=False,
    power_type="all",
    fluxes1=None,
    fluxes2=None,
    errors1=None,
    errors2=None,
    return_auxil=False,
):
    """Calculate the average cross spectrum from a list of event times or a light curve.

    If the input is a light curve, the time arrays need to be uniformly sampled
    inside GTIs (they can have gaps outside), and the fluxes need to be passed
    through the ``fluxes1`` and ``fluxes2`` arrays.
    Otherwise, times are interpeted as photon arrival times

    Parameters
    ----------
    times1 : float `np.array`
        Array of times in the sub-band
    times2 : float `np.array`
        Array of times in the reference band
    gti : [[gti00, gti01], [gti10, gti11], ...]
        common good time intervals
    segment_size : float
        length of segments
    dt : float
        Time resolution of the light curves used to produce periodograms

    Other Parameters
    ----------------
    norm : str, default "frac"
        The normalization of the periodogram. "abs" is absolute rms, "frac" is
        fractional rms, "leahy" is Leahy+83 normalization, and "none" is the
        unnormalized periodogram
    use_common_mean : bool, default True
        The mean of the light curve can be estimated in each interval, or on
        the full light curve. This gives different results (Alston+2013).
        Here we assume the mean is calculated on the full light curve, but
        the user can set ``use_common_mean`` to False to calculate it on a
        per-segment basis.
    fullspec : bool, default False
        Return the full periodogram, including negative frequencies
    silent : bool, default False
        Silence the progress bars
    power_type : str, default 'all'
        If 'all', give complex powers. If 'abs', the absolute value; if 'real',
        the real part
    fluxes1 : float `np.array`, default None
        Array of fluxes or counts per bin for channel 1
    fluxes2 : float `np.array`, default None
        Array of fluxes or counts per bin for channel 2
    errors1 : float `np.array`, default None
        Array of errors on the fluxes on channel 1
    errors2 : float `np.array`, default None
        Array of errors on the fluxes on channel 2

    Returns
    -------
    freq : `np.array`
        The periodogram frequencies
    pds : `np.array`
        The normalized periodogram powers
    n_bin : int
        the number of bins in the light curves used in each segment
    n_ave : int
        the number of averaged periodograms
    """
    if segment_size is None:
        segment_size = gti.max() - gti.min()
    n_bin = np.rint(segment_size / dt).astype(int)
    # adjust dt
    dt = segment_size / n_bin

    flux_iterable1 = get_flux_iterable_from_segments(
        times1, gti, segment_size, n_bin, fluxes=fluxes1, errors=errors1
    )
    flux_iterable2 = get_flux_iterable_from_segments(
        times2, gti, segment_size, n_bin, fluxes=fluxes2, errors=errors2
    )

    is_events = np.all([val is None for val in (fluxes1, fluxes2, errors1, errors2)])

    if (is_events
            and silent
            and use_common_mean
            and power_type == "all"
            and not fullspec
            and not return_auxil):
        results = avg_cs_from_iterables_quick(
            flux_iterable1,
            flux_iterable2,
            dt,
            norm=norm
        )

    else:
        results = avg_cs_from_iterables(
            flux_iterable1,
            flux_iterable2,
            dt,
            norm=norm,
            use_common_mean=use_common_mean,
            silent=silent,
            fullspec=fullspec,
            power_type=power_type,
            return_auxil=return_auxil,
        )
    if results is not None:
        results.meta["gti"] = gti
    return results def poisson_level(norm="frac", meanrate=None, n_ph=None, backrate=0):
    """Poisson (white)-noise level in a periodogram of pure counting noise.

    For Leahy normalization, this is:
    .. math::
        P = 2

    For the fractional r.m.s. normalization, this is
    .. math::
        P = \frac{2}{\mu}
    where :math:`\mu` is the average count rate

    For the absolute r.m.s. normalization, this is
    .. math::
        P = 2 \mu

    Finally, for the unnormalized periodogram, this is
    .. math::
        P = N_{ph}

    Parameters
    ----------
    norm : str, default "frac"
        Normalization of the periodogram. One of ["abs", "frac", "leahy", "none"]

    Other Parameters
    ----------------
    meanrate : float, default None
        Mean count rate in counts/s. Needed for r.m.s. norms (``abs`` and ``frac``)
    n_ph : float, default None
        Total number of counts in the light curve. Needed if ``norm=="none"``
    backrate : float, default 0
        Background count rate in counts/s. Optional for fractional r.m.s. norm

    Raises
    ------
    ValueError
        If the inputs are incompatible with the required normalization.

    Returns
    -------
    power_noise : float
        The Poisson noise level in the wanted normalization.

    Examples
    --------
    >>> poisson_level(norm="leahy")
    2.0
    >>> poisson_level(norm="abs", meanrate=10.)
    20.0
    >>> poisson_level(norm="frac", meanrate=10.)
    0.2
    >>> poisson_level(norm="none", n_ph=10)
    10.0
    >>> poisson_level(norm="asdfwrqfasdh3r", meanrate=10.)
    Traceback (most recent call last):
    ...
    ValueError: Unknown value for norm: asdfwrqfasdh3r...
    >>> poisson_level(norm="none", meanrate=10)
    Traceback (most recent call last):
    ...
    ValueError: Bad input parameters for norm none...
    >>> poisson_level(norm="abs", n_ph=10)
    Traceback (most recent call last):
    ...
    ValueError: Bad input parameters for norm abs...
    """
    # Various ways the parameters are wrong.
    # We want the noise in rms norm, but don't specify the mean rate.
    bad_input = norm.lower() in ["abs", "frac"] and meanrate is None
    # We want the noise in unnormalized powers, without giving n_ph.
    bad_input = bad_input or (norm.lower() == "none" and n_ph is None)

    if bad_input:
        raise ValueError(f"Bad input parameters for norm {norm}: n_ph={n_ph}, meanrate={meanrate}")

    if norm == "abs":
        return 2. * meanrate
    if norm == "frac":
        return 2. / (meanrate - backrate)**2 * meanrate
    if norm == "leahy":
        return 2.0
    if norm == "none":
        return float(n_ph)

    raise ValueError(f"Unknown value for norm: {norm}") def positive_fft_bins(n_bin, include_zero=False):
def poisson_level(norm="frac", meanrate=None, n_ph=None, backrate=0):
def normalize_frac(unnorm_power, dt, n_bin, mean_flux, background_flux=0):
def normalize_abs(unnorm_power, dt, n_bin):
def normalize_leahy_from_variance(unnorm_power, variance, n_bin):
def normalize_leahy_poisson(unnorm_power, n_ph):
def normalize_periodograms(unnorm_power, dt, n_bin, mean_flux=None, n_ph=None, variance=None, background_flux=0., norm="frac", power_type="all"):
def bias_term(power1, power2, power1_noise, power2_noise, n_ave, intrinsic_coherence=1.0):
def raw_coherence(cross_power, power1, power2, power1_noise, power2_noise, n_ave, intrinsic_coherence=1):
def _estimate_intrinsic_coherence_single(cross_power, power1, power2, power1_noise, power2_noise, n_ave):
def estimate_intrinsic_coherence(cross_power, power1, power2, power1_noise, power2_noise, n_ave):
def error_on_averaged_cross_spectrum(cross_power, seg_power, ref_power, n_ave, seg_power_noise, ref_power_noise, common_ref="False"):
def cross_to_covariance(cross_power, ref_power, ref_power_noise, delta_nu):
def _which_segment_idx_fun(binned=False, dt=None):
        def fun(*args, **kwargs):
def get_average_ctrate(times, gti, segment_size, counts=None):
def get_flux_iterable_from_segments(times, gti, segment_size, n_bin=None, fluxes=None, errors=None):
def avg_pds_from_iterable(flux_iterable, dt, norm="frac", use_common_mean=True, silent=False):
        def local_show_progress(a):
def avg_cs_from_iterables_quick(
    flux_iterable1,
    flux_iterable2,
    dt,
    norm="frac"
):
def avg_cs_from_iterables(
    flux_iterable1,
    flux_iterable2,
    dt,
    norm="frac",
    use_common_mean=True,
    silent=False,
    fullspec=False,
    power_type="all",
    return_auxil=False
):
        def local_show_progress(a):
def avg_pds_from_events(
    times,
    gti,
    segment_size,
    dt,
    norm="frac",
    use_common_mean=True,
    silent=False,
    fluxes=None,
    errors=None,
):
def avg_cs_from_events(
    times1,
    times2,
    gti,
    segment_size,
    dt,
    norm="frac",
    use_common_mean=True,
    fullspec=False,
    silent=False,
    power_type="all",
    fluxes1=None,
    fluxes2=None,
    errors1=None,
    errors2=None,
    return_auxil=False,
): def get_flux_iterable_from_segments(times, gti, segment_size, n_bin=None, fluxes=None, errors=None):
    """Get fluxes from different segments of the observation.

    If ``fluxes`` is ``None``, the input times are interpreted as events, and
    they are split into many binned series of length ``segment_size`` with
    ``n_bin`` bins.

    If ``fluxes`` is an array, the number of events corresponding to each time
    bin is taken from ``fluxes``

    Therefore, at least one of either ``n_bin`` and ``fluxes`` needs to be
    specified.

    Parameters
    ----------
    times : float `np.array`
        Array of times
    gti : [[gti00, gti01], [gti10, gti11], ...]
        good time intervals
    segment_size : float
        length of segments

    Other parameters
    ----------------
    n_bin : int, default None
        Number of bins to divide the ``segment_size`` in
    fluxes : float `np.array`, default None
        Array of fluxes.
    errors : float `np.array`, default None
        Array of error bars corresponding to the flux values above.

    Yields
    ------
    flux : `np.array`
        Array of fluxes
    err : `np.array`
        (optional) if ``errors`` is None, an array of errors in the segment

    """
    if fluxes is None and n_bin is None:
        raise ValueError(
            "At least one between fluxes (if light curve) and n_bin (if events) has to be set"
        )

    dt = None
    binned = fluxes is not None
    if binned:
        dt = np.median(np.diff(times[:100]))

    fun = _which_segment_idx_fun(binned, dt)

    for s, e, idx0, idx1 in fun(times, gti, segment_size):
        if idx1 - idx0 < 2:
            yield None
            continue
        if not binned:
            event_times = times[idx0:idx1]
            # astype here serves to avoid integer rounding issues in Windows,
            # where long is a 32-bit integer.
            cts = histogram((event_times - s).astype(float), bins=n_bin,
                            range=[0, segment_size]).astype(float)
            cts = np.array(cts)
        else:
            cts = fluxes[idx0:idx1].astype(float)
            if errors is not None:
                cts = cts, errors[idx0:idx1]

        yield cts
import copy
import warnings
import numpy as np
import scipy
import scipy.optimize
import scipy.stats
import matplotlib.pyplot as plt
from collections.abc import Iterable, Iterator, Generator
from astropy import log
from stingray.exceptions import StingrayError
from stingray.gti import bin_intervals_from_gtis, check_gtis, cross_two_gtis
from stingray.largememory import createChunkedSpectra, saveData, HAS_ZARR
from stingray.utils import genDataPath, rebin_data, rebin_data_log, simon
from .base import StingrayObject
from .events import EventList
from .lightcurve import Lightcurve
from .utils import show_progress
from .fourier import avg_cs_from_iterables, error_on_averaged_cross_spectrum
from .fourier import avg_cs_from_events, poisson_level
from .fourier import fftfreq, fft, normalize_periodograms, raw_coherence
from .fourier import get_flux_iterable_from_segments
from scipy.special import factorial

        if self.lc2.mjdref != self.lc1.mjdref:
            raise ValueError("MJDref is different in the two light curves")

        # Then check that GTIs make sense
        if self.gti is None:
            self.gti = cross_two_gtis(lc1.gti, lc2.gti)

        check_gtis(self.gti)

        if self.gti.shape[0] != 1:
            raise TypeError("Non-averaged Cross Spectra need a single Good Time Interval")

        lc1 = lc1.split_by_gti()[0]
        lc2 = lc2.split_by_gti()[0]

        # total number of photons is the sum of the
        # counts in the light curve
        self.meancounts1 = lc1.meancounts
        self.meancounts2 = lc2.meancounts
        self.nphots1 = np.float64(np.sum(lc1.counts))
        self.nphots2 = np.float64(np.sum(lc2.counts))

        self.err_dist = "poisson"
        if lc1.err_dist == "poisson":
            self.variance1 = lc1.meancounts
        else:
            self.variance1 = np.mean(lc1.counts_err) ** 2
            self.err_dist = "gauss"

        if lc2.err_dist == "poisson":
            self.variance2 = lc2.meancounts
        else:
            self.variance2 = np.mean(lc2.counts_err) ** 2
            self.err_dist = "gauss"

        if lc1.n != lc2.n:
            raise StingrayError("Light curves do not have same number of time bins per segment.")

        # If dt differs slightly, its propagated error must not be more than
        # 1/100th of the bin
        if not np.isclose(lc1.dt, lc2.dt, rtol=0.01 * lc1.dt / lc1.tseg):
            raise StingrayError("Light curves do not have same time binning dt.")

        # In case a small difference exists, ignore it
        lc1.dt = lc2.dt

        self.dt = lc1.dt
        self.n = lc1.n

        # the frequency resolution
        self.df = 1.0 / lc1.tseg

        # the number of averaged periodograms in the final output
        # This should *always* be 1 here
        self.m = 1

        # make the actual Fourier transform and compute cross spectrum
        self.freq, self.unnorm_power = self._fourier_cross(lc1, lc2, fullspec)

        # If co-spectrum is desired, normalize here. Otherwise, get raw back
        # with the imaginary part still intact.
        self.power = self._normalize_crossspectrum(self.unnorm_power)

        if lc1.err_dist.lower() != lc2.err_dist.lower():
            simon(
                "Your lightcurves have different statistics."
                "The errors in the Crossspectrum will be incorrect."
            )
        elif lc1.err_dist.lower() != "poisson":
            simon(
                "Looks like your lightcurve statistic is not poisson."
                "The errors in the Powerspectrum will be incorrect."
            )

        if self.__class__.__name__ in ["Powerspectrum", "AveragedPowerspectrum"]:
            self.power_err = self.power / np.sqrt(self.m)
        elif self.__class__.__name__ in ["Crossspectrum", "AveragedCrossspectrum"]:
            # This is clearly a wild approximation.
            simon(
                "Errorbars on cross spectra are not thoroughly tested. "
                "Please report any inconsistencies."
            )
            unnorm_power_err = np.sqrt(2) / np.sqrt(self.m)  # Leahy-like
            unnorm_power_err /= 2 / np.sqrt(self.nphots1 * self.nphots2)
            unnorm_power_err += np.zeros_like(self.power)

            self.power_err = self._normalize_crossspectrum(unnorm_power_err)
        else:
            self.power_err = np.zeros(len(self.power))

    def _fourier_cross(self, lc1, lc2, fullspec=False):
        """
        Fourier transform the two light curves, then compute the cross spectrum.
        Computed as CS = lc1 x lc2* (where lc2 is the one that gets
        complex-conjugated). The user has the option to either get just the
        positive frequencies or the full spectrum.

        Parameters
        ----------
        lc1: :class:`stingray.Lightcurve` object
            One light curve to be Fourier transformed. Ths is the band of
            interest or channel of interest.

        lc2: :class:`stingray.Lightcurve` object
            Another light curve to be Fourier transformed.
            This is the reference band.

        fullspec: boolean. Default is False.
            If True, return the whole array of frequencies, or only positive frequencies (False).

        Returns
        -------
        fr: numpy.ndarray
            The squared absolute value of the Fourier amplitudes

        """
        fourier_1 = fft(lc1.counts)  # do Fourier transform 1
        fourier_2 = fft(lc2.counts)  # do Fourier transform 2


        freqs = fftfreq(lc1.n, lc1.dt)