class StingrayTimeseries(StingrayObject):
    def to_astropy_timeseries(self):
        """Save the ``StingrayTimeseries`` to an ``Astropy`` timeseries.

        Array attributes (time, pi, energy, etc.) are converted
        into columns, while meta attributes (mjdref, gti, etc.)
        are saved into the ``meta`` dictionary.

        Returns
        -------
        ts : `astropy.timeseries.TimeSeries`
            A ``TimeSeries`` object with the array attributes as columns,
            and the meta attributes in the `meta` dictionary
        """
        from astropy.timeseries import TimeSeries
        from astropy.time import TimeDelta
        from astropy import units as u

        data = {}
        array_attrs = self.array_attrs()

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

        if data == {}:
            data = None

        if self.time is not None and np.size(self.time) > 0:
            times = TimeDelta(self.time * u.s)
            ts = TimeSeries(data=data, time=times)
        else:
            ts = TimeSeries()

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

        return ts

    @classmethod
    def from_astropy_timeseries(cls, ts):
        """Create a `StingrayTimeseries` from data in an Astropy TimeSeries

        The timeseries has to define at least a column called time,
        the rest of columns will form the array attributes of the
        new event list, while the attributes in table.meta will
        form the new meta attributes of the event list.

        It is strongly advisable to define such attributes and columns
        using the standard attributes of EventList: time, pi, energy, gti etc.

        Parameters
        ----------
        ts : `astropy.timeseries.TimeSeries`
            A ``TimeSeries`` object with the array attributes as columns,
            and the meta attributes in the `meta` dictionary

        Returns
        -------
        ts : `StingrayTimeseries`
            Timeseries object
        """

        time = ts["time"]
        mjdref = None
        if "mjdref" in ts.meta:
            mjdref = ts.meta["mjdref"]

        time, mjdref = interpret_times(time, mjdref)
        cls.time = np.asarray(time)

        array_attrs = ts.colnames
        for key, val in ts.meta.items():
            setattr(cls, key, val)

        for attr in array_attrs:
            if attr == "time":
                continue
            setattr(cls, attr, ts[attr])

        return cls

    def change_mjdref(self, new_mjdref):
        """Change the MJD reference time (MJDREF) of the time series

        The times of the time series will be shifted in order to be referred to
        this new MJDREF

        Parameters
        ----------
        new_mjdref : float
            New MJDREF

        Returns
        -------
        new_lc : :class:`StingrayTimeseries` object
            The new time series, shifted by MJDREF
        """
        time_shift = (self.mjdref - new_mjdref) * 86400

        ts = self.shift(time_shift)
        ts.mjdref = new_mjdref
        return ts

    def shift(self, time_shift):
        """Shift the time and the GTIs by the same amount

        Parameters
        ----------
        time_shift: float
            The time interval by which the light curve will be shifted (in
            the same units as the time array in :class:`Lightcurve`

        Returns
        -------
        ts : ``StingrayTimeseries`` object
            The new time series shifted by ``time_shift``

        """
        ts = copy.deepcopy(self)
        ts.time = np.asarray(ts.time) + time_shift
        if hasattr(ts, "gti"):
            ts.gti = np.asarray(ts.gti) + time_shift

        return ts 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 create_gti_mask(time, gtis, safe_interval=None, min_length=0,
                    return_new_gtis=False, dt=None, epsilon=0.001):
    """Create GTI mask.

    Assumes that no overlaps are present between GTIs

    Parameters
    ----------
    time : numpy.ndarray
        An array of time stamps

    gtis : ``[[g0_0, g0_1], [g1_0, g1_1], ...]``, float array-like
        The list of GTIs


    Other parameters
    ----------------
    safe_interval : float or ``[float, float]``, default None
        A safe interval to exclude at both ends (if single float) or the start
        and the end (if pair of values) of GTIs. If None, no safe interval
        is applied to data.

    min_length : float
        An optional minimum length for the GTIs to be applied. Only GTIs longer
        than ``min_length`` will be considered when creating the mask.

    return_new_gtis : bool
        If ``True```, return the list of new GTIs (if ``min_length > 0``)

    dt : float
        Time resolution of the data, i.e. the interval between time stamps

    epsilon : float
        fraction of ``dt`` that is tolerated at the borders of a GTI

    Returns
    -------
    mask : bool array
        A mask labelling all time stamps that are included in the GTIs versus those that are not.

    new_gtis : ``Nx2`` array
        An array of new GTIs created by this function
    """
    gtis = np.array(gtis, dtype=np.longdouble)
    if time.size == 0:
        raise ValueError("Passing an empty time array to create_gti_mask")
    if gtis.size == 0:
        raise ValueError("Passing an empty GTI array to create_gti_mask")

    mask = np.zeros(len(time), dtype=bool)

    if min_length > 0:
        lengths = gtis[:, 1] - gtis[:, 0]
        good = lengths >= np.max(min_length, dt)
        if np.all(~good):
            warnings.warn("No GTIs longer than "
                          "min_length {}".format(min_length))
            return mask
        gtis = gtis[good]

    if not HAS_NUMBA:
        return create_gti_mask_complete(time, gtis,
                                        safe_interval=safe_interval,
                                        min_length=min_length,
                                        return_new_gtis=return_new_gtis,
                                        dt=dt, epsilon=epsilon)

    check_gtis(gtis)

    dt = apply_function_if_none(dt, time,
                                lambda x: np.median(np.diff(x)))
    epsilon_times_dt = epsilon * dt
    gtis_new = copy.deepcopy(gtis)
    gti_mask = np.zeros(len(gtis), dtype=bool)

    if safe_interval is not None:
        if not isinstance(safe_interval, Iterable):
            safe_interval = np.array([safe_interval, safe_interval])
        # These are the gtis that will be returned (filtered!). They are only
        # modified by the safe intervals
        gtis_new[:, 0] = gtis[:, 0] + safe_interval[0]
        gtis_new[:, 1] = gtis[:, 1] - safe_interval[1]

    # These are false gtis, they contain a few boundary modifications
    # in order to simplify the calculation of the mask, but they will _not_
    # be returned.
    gtis_to_mask = copy.deepcopy(gtis_new)
    gtis_to_mask[:, 0] = gtis_new[:, 0] - epsilon_times_dt + dt / 2
    gtis_to_mask[:, 1] = gtis_new[:, 1] + epsilon_times_dt - dt / 2

    mask, gtimask = \
        create_gti_mask_jit((time - time[0]).astype(np.float64),
                            (gtis_to_mask - time[0]).astype(np.float64),
                            mask, gti_mask=gti_mask,
                            min_length=float(
            min_length - 2 * (1 + epsilon) * dt
        ))

    if return_new_gtis:
        return mask, gtis_new[gtimask]

    return mask 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 gti_border_bins(gtis, time, dt=None, epsilon=0.001):
    """Find the indices in a time array corresponding to the borders of GTIs.

    GTIs shorter than the bin time are not returned.

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

    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 of each GTI

    spectrum_stop_bins : array-like
        List of stop bins of each GTI. The elements corresponding to these bins
        should *not* be included.

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

    >>> gti_border_bins([[16., 18.]], times)
    Traceback (most recent call last):
        ...
    ValueError: Invalid time interval for the given GTIs

    >>> start_bins, stop_bins = gti_border_bins(
    ...    [[0, 5], [6, 8]], times)

    >>> np.allclose(start_bins, [0, 6])
    True
    >>> np.allclose(stop_bins, [5, 8])
    True
    >>> np.allclose(times[start_bins[0]:stop_bins[0]], [0.5, 1.5, 2.5, 3.5, 4.5])
    True
    >>> np.allclose(times[start_bins[1]:stop_bins[1]], [6.5, 7.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

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

    spectrum_start_bins = []
    spectrum_stop_bins = []

    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):
        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

        spectrum_start_bins.append(startbin)
        spectrum_stop_bins.append(stopbin)

    return np.array(spectrum_start_bins), np.array(spectrum_stop_bins) def join_gtis(gti0, gti1):
    """Union of two GTIs.

    If GTIs are mutually exclusive, it calls ``append_gtis``. Otherwise we put
    the extremes of partially overlapping GTIs on an ideal line and look at the
    number of opened and closed intervals. When the number of closed and opened
    intervals is the same, the full GTI is complete and we close it.

    In practice, we assign to each opening time of a GTI the value ``-1``, and
    the value ``1`` to each closing time; when the cumulative sum is zero, the
    GTI has ended. The timestamp after each closed GTI is the start of a new
    one.

    ::

        (cumsum)   -1   -2         -1   0   -1 -2           -1  -2  -1        0
        GTI A      |-----:----------|   :    |--:------------|   |---:--------|
        FINAL GTI  |-----:--------------|    |--:--------------------:--------|
        GTI B            |--------------|       |--------------------|

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

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

    Returns
    -------
    gti: 2-d float array
        The newly created GTI
    """

    gti0 = np.asarray(gti0)
    gti1 = np.asarray(gti1)

    # Check if independently GTIs are well behaved.
    check_gtis(gti0)
    check_gtis(gti1)

    if check_separate(gti0, gti1):
        return append_gtis(gti0, gti1)

    g0 = gti0.flatten()
    # Opening GTI: type = 1; Closing: type = -1
    g0_type = np.asarray(list(zip(-np.ones(int(len(g0) / 2)),
                                  np.ones(int(len(g0) / 2)))))
    g1 = gti1.flatten()
    g1_type = np.asarray(list(zip(-np.ones(int(len(g1) / 2)),
                                  np.ones(int(len(g1) / 2)))))

    g_all = np.append(g0, g1)
    g_type_all = np.append(g0_type, g1_type)
    order = np.argsort(g_all)
    g_all = g_all[order]
    g_type_all = g_type_all[order]

    sums = np.cumsum(g_type_all)

    # Where the cumulative sum is zero, we close the GTI
    closing_bins = sums == 0
    # The next element in the sequence is the start of the new GTI. In the case
    # of the last element, the next is the first. Numpy.roll gives this for
    # free.
    starting_bins = np.roll(closing_bins, 1)

    starting_times = g_all[starting_bins]
    closing_times = g_all[closing_bins]

    final_gti = []
    for start, stop in zip(starting_times, closing_times):
        final_gti.append([start, stop])

    return np.sort(final_gti, axis=0)     HAS_PYFFTW = True
    HAS_PYFFTW = False
HAS_NUMBA = False
    HAS_NUMBA = True
    L = len(y)
    D = sparse.csc_matrix(np.diff(np.eye(L), 2))
        W = sparse.spdiags(w, 0, L, L)
        Z = W + lam * D.dot(D.transpose())
    N = len(lc.counts)
    N = abs(N)
        H = np.lib.format.open_memmap(
            tmp, mode="w+", dtype=dtype, shape=shape
        )
        H = np.zeros(shape, dtype=dtype)
    H = _allocate_array_or_memmap(bins, np.uint64, use_memmap=use_memmap, tmp=tmp)
    H = _allocate_array_or_memmap(bins, np.uint64, use_memmap=use_memmap, tmp=tmp)
    H = _allocate_array_or_memmap(bins, np.double, use_memmap=use_memmap, tmp=tmp)
    H = _allocate_array_or_memmap(bins, np.double, use_memmap=use_memmap, tmp=tmp)
    H = _allocate_array_or_memmap(bins, np.uint64, use_memmap=use_memmap, tmp=tmp)
    class jit(object):
    class njit(object):
    class vectorize(object):
        def __init__(self, *args, **kwargs):
        def __call__(self, func):
            def wrapped_f(*args, **kwargs):
        def __init__(self, *args, **kwargs):
        def __call__(self, func):
            def wrapped_f(*args, **kwargs):
        def __init__(self, *args, **kwargs):
        def __call__(self, func):
    def generic(x, y=None):
    def prange(x):
    def show_progress(a):
    def mad(data, c=0.6745, axis=None):
def _root_squared_mean(array):
def simon(message, **kwargs):
def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
def rebin_data_log(x, y, f, y_err=None, dx=None):
def apply_function_if_none(variable, value, func):
def assign_value_if_none(value, default):
def look_for_array_in_array(array1, array2):
def is_string(s):
def is_iterable(var):
def order_list_of_arrays(data, order):
def optimal_bin_time(fftlen, tbin):
def contiguous_regions(condition):
def is_int(obj):
def get_random_state(random_state=None):
def _offset(x, off):
def offset_fit(x, y, offset_start=0):
def _als(y, lam, p, niter=10):
def baseline_als(x, y, lam=None, p=None, niter=10, return_baseline=False,
                 offset_correction=False):
def excess_variance(lc, normalization='fvar'):
def create_window(N, window_type='uniform'):
def poisson_symmetrical_errors(counts):
def standard_error(xs, mean):
def nearest_power_of_two(x):
def find_nearest(array, value):
def genDataPath(dir_path):
def check_iterables_close(iter0, iter1, **kwargs):
def compute_bin(x, bin_edges):
def _hist1d_numba_seq(H, tracks, bins, ranges):
def _allocate_array_or_memmap(shape, dtype, use_memmap=False, tmp=None):
def hist1d_numba_seq(a, bins, ranges, use_memmap=False, tmp=None):
def _hist2d_numba_seq(H, tracks, bins, ranges):
def hist2d_numba_seq(x, y, bins, ranges, use_memmap=False, tmp=None):
def _hist3d_numba_seq(H, tracks, bins, ranges):
def hist3d_numba_seq(tracks, bins, ranges, use_memmap=False, tmp=None):
def _hist2d_numba_seq_weight(H, tracks, weights, bins, ranges):
def hist2d_numba_seq_weight(x, y, weights, bins, ranges, use_memmap=False, tmp=None):
def _hist3d_numba_seq_weight(H, tracks, weights, bins, ranges):
def hist3d_numba_seq_weight(tracks, weights, bins, ranges, use_memmap=False, tmp=None):
def _index_arr(a, ix_arr):
def _index_set_arr(a, ix_arr, val):
def _histnd_numba_seq(H, tracks, bins, ranges, slice_int):
def histnd_numba_seq(tracks, bins, ranges, use_memmap=False, tmp=None):
    def histogram2d(*args, **kwargs):
    def histogram(*args, **kwargs):
    def histogram2d(*args, **kwargs):
    def histogram(*args, **kwargs):
def equal_count_energy_ranges(energies, n_ranges, emin=None, emax=None):
def sum_if_not_none_or_initialize(A, B):
def assign_if_not_finite(value, default): def lcurve_from_fits(
    fits_file,
    gtistring="GTI",
    timecolumn="TIME",
    ratecolumn=None,
    ratehdu=1,
    fracexp_limit=0.9,
    outfile=None,
    noclobber=False,
    outdir=None,
):
    """Load a lightcurve from a fits file.

    .. note ::
        FITS light curve handling is still under testing.
        Absolute times might be incorrect depending on the light curve format.

    Parameters
    ----------
    fits_file : str
        File name of the input light curve in FITS format

    Returns
    -------
    data : dict
        Dictionary containing all information needed to create a
        :class:`stingray.Lightcurve` object

    Other Parameters
    ----------------
    gtistring : str
        Name of the GTI extension in the FITS file
    timecolumn : str
        Name of the column containing times in the FITS file
    ratecolumn : str
        Name of the column containing rates in the FITS file
    ratehdu : str or int
        Name or index of the FITS extension containing the light curve
    fracexp_limit : float
        Minimum exposure fraction allowed
    noclobber : bool
        If True, do not overwrite existing files
    """
    warnings.warn(
        """WARNING! FITS light curve handling is still under testing.
        Absolute times might be incorrect."""
    )
    # TODO:
    # treat consistently TDB, UTC, TAI, etc. This requires some documentation
    # reading. For now, we assume TDB
    from astropy.io import fits as pf
    from astropy.time import Time
    import numpy as np
    from stingray.gti import create_gti_from_condition

    lchdulist = pf.open(fits_file)
    lctable = lchdulist[ratehdu].data

    # Units of header keywords
    tunit = lchdulist[ratehdu].header["TIMEUNIT"]

    try:
        mjdref = high_precision_keyword_read(
            lchdulist[ratehdu].header, "MJDREF"
        )
        mjdref = Time(mjdref, scale="tdb", format="mjd")
    except Exception:
        mjdref = None

    try:
        instr = lchdulist[ratehdu].header["INSTRUME"]
    except Exception:
        instr = "EXTERN"

    # ----------------------------------------------------------------
    # Trying to comply with all different formats of fits light curves.
    # It's a madness...
    try:
        tstart = high_precision_keyword_read(
            lchdulist[ratehdu].header, "TSTART"
        )
        tstop = high_precision_keyword_read(lchdulist[ratehdu].header, "TSTOP")
    except Exception:  # pragma: no cover
        raise (Exception("TSTART and TSTOP need to be specified"))

    # For nulccorr lcs this whould work

    timezero = high_precision_keyword_read(
        lchdulist[ratehdu].header, "TIMEZERO"
    )
    # Sometimes timezero is "from tstart", sometimes it's an absolute time.
    # This tries to detect which case is this, and always consider it
    # referred to tstart
    timezero = assign_value_if_none(timezero, 0)

    # for lcurve light curves this should instead work
    if tunit == "d":
        # TODO:
        # Check this. For now, I assume TD (JD - 2440000.5).
        # This is likely wrong
        timezero = Time(2440000.5 + timezero, scale="tdb", format="jd")
        tstart = Time(2440000.5 + tstart, scale="tdb", format="jd")
        tstop = Time(2440000.5 + tstop, scale="tdb", format="jd")
        # if None, use NuSTAR defaulf MJDREF
        mjdref = assign_value_if_none(
            mjdref,
            Time(
                np.longdouble("55197.00076601852"), scale="tdb", format="mjd"
            ),
        )

        timezero = (timezero - mjdref).to("s").value
        tstart = (tstart - mjdref).to("s").value
        tstop = (tstop - mjdref).to("s").value

    if timezero > tstart:
        timezero -= tstart

    time = np.array(lctable.field(timecolumn), dtype=np.longdouble)
    if time[-1] < tstart:
        time += timezero + tstart
    else:
        time += timezero

    try:
        dt = high_precision_keyword_read(lchdulist[ratehdu].header, "TIMEDEL")
        if tunit == "d":
            dt *= 86400
    except Exception:
        warnings.warn(
            "Assuming that TIMEDEL is the median difference between the"
            " light curve times",
            AstropyUserWarning,
        )
        # Avoid NaNs
        good = time == time
        dt = np.median(np.diff(time[good]))

    # ----------------------------------------------------------------
    if ratecolumn is None:
        for name in ["RATE", "RATE1", "COUNTS"]:
            if name in lctable.names:
                ratecolumn = name
                break
        else:  # pragma: no cover
            raise ValueError(
                "None of the accepted rate columns were found in the file")

    rate = np.array(lctable.field(ratecolumn), dtype=float)

    errorcolumn = "ERROR"
    if ratecolumn == "RATE1":
        errorcolumn = "ERROR1"

    try:
        rate_e = np.array(lctable.field(errorcolumn), dtype=np.longdouble)
    except Exception:
        rate_e = np.zeros_like(rate)

    if "RATE" in ratecolumn:
        rate *= dt
        rate_e *= dt

    try:
        fracexp = np.array(lctable.field("FRACEXP"), dtype=np.longdouble)
    except Exception:
        fracexp = np.ones_like(rate)

    good_intervals = (
        (rate == rate) * (fracexp >= fracexp_limit) * (fracexp <= 1)
    )

    rate[good_intervals] /= fracexp[good_intervals]
    rate_e[good_intervals] /= fracexp[good_intervals]

    rate[~good_intervals] = 0

    try:
        gtitable = lchdulist[gtistring].data
        gti_list = np.array(
            [
                [a, b]
                for a, b in zip(
                    gtitable.field("START"), gtitable.field("STOP")
                )
            ],
            dtype=np.longdouble,
        )
    except Exception:
        gti_list = create_gti_from_condition(time, good_intervals)

    lchdulist.close()

    res = {"time": time,
           "counts": rate,
           "err": rate_e,
           "gti": gti_list,
           "mjdref": mjdref.mjd,
           "dt": dt,
           "instr": instr,
           "header": lchdulist[ratehdu].header.tostring()}
    return res
import os
import logging
import warnings
import pickle
import numpy as np
import matplotlib.pyplot as plt
import stingray.utils as utils
from astropy.table import Table
from astropy.time import TimeDelta, Time
from astropy import units as u
from stingray.base import StingrayTimeseries
from stingray.exceptions import StingrayError
from stingray.gti import (bin_intervals_from_gtis, check_gtis, create_gti_mask,
                          cross_two_gtis, gti_border_bins, join_gtis)
from stingray.utils import (assign_value_if_none, baseline_als,
                            poisson_symmetrical_errors, simon, interpret_times)
from stingray.io import lcurve_from_fits
            from collections import Counter
            from lightkurve import LightCurve as lk
            from astropy.timeseries import TimeSeries

        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)