astrobase/lcmath.py
def fill_magseries_gaps(times, mags, errs,
                        fillgaps=0.0,
                        sigclip=3.0,
                        magsarefluxes=False,
                        filterwindow=11,
                        forcetimebin=None,
                        verbose=True):
    '''This fills in gaps in a light curve.

    This is mainly intended for use in ACF period-finding, but maybe useful
    otherwise (i.e. when we figure out ARMA stuff for LCs). The main steps here
    are:

    - normalize the light curve to zero
    - remove giant outliers
    - interpolate gaps in the light curve
      (since ACF requires evenly spaced sampling)

    From McQuillan+ 2013a (https://doi.org/10.1093/mnras/stt536):

        "The ACF calculation requires the light curves to be regularly sampled
        and normalized to zero. We divided the flux in each quarter by its
        median and subtracted unity. Gaps in the light curve longer than the
        Kepler long cadence were filled using linear interpolation with added
        white Gaussian noise. This noise level was estimated using the variance
        of the residuals following subtraction of a smoothed version of the
        flux. To smooth the flux, we applied an iterative non-linear filter
        which consists of a median filter followed by a boxcar filter, both with
        11-point windows, with iterative 3σ clipping of outliers."

    Parameters
    ----------

    times,mags,errs : np.array
        The magnitude/flux time-series and associated measurement errors to
        operate on. Non-finite elements will be removed from these arrays. At
        least 10 elements in each array are required for this function to
        operate.

    fillgaps : {'noiselevel', 'nan'} or float
        If `fillgap='noiselevel'`, fills the gaps with the noise level obtained
        via the procedure above. If `fillgaps='nan'`, fills the gaps with
        `np.nan`. Otherwise, if `fillgaps` is a float, will use that value to
        fill the gaps. The default is to fill the gaps with 0.0 (as in
        McQuillan+ 2014) to "...prevent them contributing to the ACF".

    sigclip : float or int or sequence of two floats/ints or None
        If a single float or int, a symmetric sigma-clip will be performed using
        the number provided as the sigma-multiplier to cut out from the input
        time-series.

        If a list of two ints/floats is provided, the function will perform an
        'asymmetric' sigma-clip. The first element in this list is the sigma
        value to use for fainter flux/mag values; the second element in this
        list is the sigma value to use for brighter flux/mag values. For
        example, `sigclip=[10., 3.]`, will sigclip out greater than 10-sigma
        dimmings and greater than 3-sigma brightenings. Here the meaning of
        "dimming" and "brightening" is set by *physics* (not the magnitude
        system), which is why the `magsarefluxes` kwarg must be correctly set.

        If `sigclip` is None, no sigma-clipping will be performed, and the
        time-series (with non-finite elems removed) will be passed through to
        the output.

    magsareflux : bool
        True if your "mags" are in fact fluxes, i.e. if "fainter" corresponds to
        `mags` getting smaller.

    filterwindow : int
        The number of time-series points to include in the Savitsky-Golay filter
        operation when smoothing the light curve. This should be an odd integer.

    forcetimebin : float or None
        If `forcetimebin` is a float, this value will be used to generate the
        interpolated time series, effectively binning the light curve to this
        cadence. If `forcetimebin` is None, the mode of the gaps (the forward
        difference between successive time values in `times`) in the provided
        light curve will be used as the effective cadence. NOTE: `forcetimebin`
        must be in the same units as `times`, e.g. if times are JD then
        `forcetimebin` must be in days as well

    verbose : bool
        If this is True, will indicate progress at various stages in the
        operation.

    Returns
    -------

    dict
        A dict of the following form is returned::

            {'itimes': the interpolated time values after gap-filling,
             'imags': the interpolated mag/flux values after gap-filling,
             'ierrs': the interpolated mag/flux values after gap-filling,
             'cadence': the cadence of the output mag/flux time-series}

    '''

    # remove nans
    finind = np.isfinite(times) & np.isfinite(mags) & np.isfinite(errs)
    ftimes, fmags, ferrs = times[finind], mags[finind], errs[finind]

    # remove zero errs
    nzind = np.nonzero(ferrs)
    ftimes, fmags, ferrs = ftimes[nzind], fmags[nzind], ferrs[nzind]

    # sigma-clip
    stimes, smags, serrs = sigclip_magseries(ftimes, fmags, ferrs,
                                             magsarefluxes=magsarefluxes,
                                             sigclip=sigclip)

    # normalize to zero
    if magsarefluxes:
        smags = smags / np.median(smags) - 1.0
    else:
        smags = smags - np.median(smags)

    if isinstance(fillgaps, float):

        gapfiller = fillgaps

    elif isinstance(fillgaps, str) and fillgaps == 'noiselevel':

        # figure out the gaussian noise level by subtracting a Savitsky-Golay
        # filtered version of the light curve
        smoothed = smags - savgol_filter(smags, filterwindow, 2)
        noiselevel = 1.483 * np.median(np.abs(smoothed - np.median(smoothed)))
        gapfiller = noiselevel

    elif isinstance(fillgaps, str) and fillgaps == 'nan':

        gapfiller = np.nan

    # figure out the gap size and where to interpolate. we do this by figuring
    # out the most common gap (this should be the cadence). to do this, we need
    # to calculate the mode of the gap distribution.

    # get the gaps
    gaps = np.diff(stimes)

    # just use scipy.stats.mode instead of our hacked together nonsense earlier.
    gapmoderes = scipy.stats.mode(gaps)
    gapmode = gapmoderes[0].item()

    LOGINFO('auto-cadence for mag series: %.5f' % gapmode)

    # sort the gaps
    if forcetimebin:
        LOGWARNING('forcetimebin is set, forcing cadence to %.5f' %
                   forcetimebin)
        gapmode = forcetimebin

    if gapmode == 0.0:
        LOGERROR('the smallest cadence of this light curve appears to be 0.0, '
                 'the automatic cadence finder probably failed. '
                 'try setting forcetimebin?')
        return None

    starttime, endtime = np.min(stimes), np.max(stimes)
    ntimes = int(np.ceil((endtime - starttime)/gapmode) + 1)
    if verbose:
        LOGINFO('generating new time series with %s measurements' % ntimes)

    # first, generate the full time series
    interpolated_times = np.linspace(starttime, endtime, ntimes)
    interpolated_mags = np.full_like(interpolated_times, gapfiller)
    interpolated_errs = np.full_like(interpolated_times, gapfiller)

    for ind, itime in enumerate(interpolated_times[:-1]):

        nextitime = itime + gapmode
        # find the mags between this and the next time bin
        itimeind = np.where((stimes > itime) & (stimes < nextitime))

        # if there's more than one elem in this time bin, median them
        if itimeind[0].size > 1:

            interpolated_mags[ind] = np.median(smags[itimeind[0]])
            interpolated_errs[ind] = np.median(serrs[itimeind[0]])

        # otherwise, if there's only one elem in this time bin, take it
        elif itimeind[0].size == 1:

            interpolated_mags[ind] = smags[itimeind[0]]
            interpolated_errs[ind] = serrs[itimeind[0]]


    return {'itimes':interpolated_times,
            'imags':interpolated_mags,
            'ierrs':interpolated_errs,
            'cadence':gapmode}
def fill_magseries_gaps(times, mags, errs,
                        fillgaps=0.0,
                        sigclip=3.0,
                        magsarefluxes=False,
                        filterwindow=11,
                        forcetimebin=None,
                        verbose=True):
    '''This fills in gaps in a light curve.

    This is mainly intended for use in ACF period-finding, but maybe useful
    otherwise (i.e. when we figure out ARMA stuff for LCs). The main steps here
    are:

    - normalize the light curve to zero
    - remove giant outliers
    - interpolate gaps in the light curve
      (since ACF requires evenly spaced sampling)

    From McQuillan+ 2013a (https://doi.org/10.1093/mnras/stt536):

        "The ACF calculation requires the light curves to be regularly sampled
        and normalized to zero. We divided the flux in each quarter by its
        median and subtracted unity. Gaps in the light curve longer than the
        Kepler long cadence were filled using linear interpolation with added
        white Gaussian noise. This noise level was estimated using the variance
        of the residuals following subtraction of a smoothed version of the
        flux. To smooth the flux, we applied an iterative non-linear filter
        which consists of a median filter followed by a boxcar filter, both with
        11-point windows, with iterative 3σ clipping of outliers."

    Parameters
    ----------

    times,mags,errs : np.array
        The magnitude/flux time-series and associated measurement errors to
        operate on. Non-finite elements will be removed from these arrays. At
        least 10 elements in each array are required for this function to
        operate.

    fillgaps : {'noiselevel', 'nan'} or float
        If `fillgap='noiselevel'`, fills the gaps with the noise level obtained
        via the procedure above. If `fillgaps='nan'`, fills the gaps with
        `np.nan`. Otherwise, if `fillgaps` is a float, will use that value to
        fill the gaps. The default is to fill the gaps with 0.0 (as in
        McQuillan+ 2014) to "...prevent them contributing to the ACF".

    sigclip : float or int or sequence of two floats/ints or None
        If a single float or int, a symmetric sigma-clip will be performed using
        the number provided as the sigma-multiplier to cut out from the input
        time-series.

        If a list of two ints/floats is provided, the function will perform an
        'asymmetric' sigma-clip. The first element in this list is the sigma
        value to use for fainter flux/mag values; the second element in this
        list is the sigma value to use for brighter flux/mag values. For
        example, `sigclip=[10., 3.]`, will sigclip out greater than 10-sigma
        dimmings and greater than 3-sigma brightenings. Here the meaning of
        "dimming" and "brightening" is set by *physics* (not the magnitude
        system), which is why the `magsarefluxes` kwarg must be correctly set.

        If `sigclip` is None, no sigma-clipping will be performed, and the
        time-series (with non-finite elems removed) will be passed through to
        the output.

    magsareflux : bool
        True if your "mags" are in fact fluxes, i.e. if "fainter" corresponds to
        `mags` getting smaller.

    filterwindow : int
        The number of time-series points to include in the Savitsky-Golay filter
        operation when smoothing the light curve. This should be an odd integer.

    forcetimebin : float or None
        If `forcetimebin` is a float, this value will be used to generate the
        interpolated time series, effectively binning the light curve to this
        cadence. If `forcetimebin` is None, the mode of the gaps (the forward
        difference between successive time values in `times`) in the provided
        light curve will be used as the effective cadence. NOTE: `forcetimebin`
        must be in the same units as `times`, e.g. if times are JD then
        `forcetimebin` must be in days as well

    verbose : bool
        If this is True, will indicate progress at various stages in the
        operation.

    Returns
    -------

    dict
        A dict of the following form is returned::

            {'itimes': the interpolated time values after gap-filling,
             'imags': the interpolated mag/flux values after gap-filling,
             'ierrs': the interpolated mag/flux values after gap-filling,
             'cadence': the cadence of the output mag/flux time-series}

    '''

    # remove nans
    finind = np.isfinite(times) & np.isfinite(mags) & np.isfinite(errs)
    ftimes, fmags, ferrs = times[finind], mags[finind], errs[finind]

    # remove zero errs
    nzind = np.nonzero(ferrs)
    ftimes, fmags, ferrs = ftimes[nzind], fmags[nzind], ferrs[nzind]

    # sigma-clip
    stimes, smags, serrs = sigclip_magseries(ftimes, fmags, ferrs,
                                             magsarefluxes=magsarefluxes,
                                             sigclip=sigclip)

    # normalize to zero
    if magsarefluxes:
        smags = smags / np.median(smags) - 1.0
    else:
        smags = smags - np.median(smags)

    if isinstance(fillgaps, float):

        gapfiller = fillgaps

    elif isinstance(fillgaps, str) and fillgaps == 'noiselevel':

        # figure out the gaussian noise level by subtracting a Savitsky-Golay
        # filtered version of the light curve
        smoothed = smags - savgol_filter(smags, filterwindow, 2)
        noiselevel = 1.483 * np.median(np.abs(smoothed - np.median(smoothed)))
        gapfiller = noiselevel

    elif isinstance(fillgaps, str) and fillgaps == 'nan':

        gapfiller = np.nan

    # figure out the gap size and where to interpolate. we do this by figuring
    # out the most common gap (this should be the cadence). to do this, we need
    # to calculate the mode of the gap distribution.

    # get the gaps
    gaps = np.diff(stimes)

    # just use scipy.stats.mode instead of our hacked together nonsense earlier.
    gapmoderes = scipy.stats.mode(gaps)
    gapmode = gapmoderes[0].item()

    LOGINFO('auto-cadence for mag series: %.5f' % gapmode)

    # sort the gaps
    if forcetimebin:
        LOGWARNING('forcetimebin is set, forcing cadence to %.5f' %
                   forcetimebin)
        gapmode = forcetimebin

    if gapmode == 0.0:
        LOGERROR('the smallest cadence of this light curve appears to be 0.0, '
                 'the automatic cadence finder probably failed. '
                 'try setting forcetimebin?')
        return None

    starttime, endtime = np.min(stimes), np.max(stimes)
    ntimes = int(np.ceil((endtime - starttime)/gapmode) + 1)
    if verbose:
        LOGINFO('generating new time series with %s measurements' % ntimes)

    # first, generate the full time series
    interpolated_times = np.linspace(starttime, endtime, ntimes)
    interpolated_mags = np.full_like(interpolated_times, gapfiller)
    interpolated_errs = np.full_like(interpolated_times, gapfiller)

    for ind, itime in enumerate(interpolated_times[:-1]):

        nextitime = itime + gapmode
        # find the mags between this and the next time bin
        itimeind = np.where((stimes > itime) & (stimes < nextitime))

        # if there's more than one elem in this time bin, median them
        if itimeind[0].size > 1:

            interpolated_mags[ind] = np.median(smags[itimeind[0]])
            interpolated_errs[ind] = np.median(serrs[itimeind[0]])

        # otherwise, if there's only one elem in this time bin, take it
        elif itimeind[0].size == 1:

            interpolated_mags[ind] = smags[itimeind[0]]
            interpolated_errs[ind] = serrs[itimeind[0]]


    return {'itimes':interpolated_times,
            'imags':interpolated_mags,
            'ierrs':interpolated_errs,
            'cadence':gapmode}
This fills in gaps in a light curve.

    This is mainly intended for use in ACF period-finding, but maybe useful
    otherwise (i.e. when we figure out ARMA stuff for LCs). The main steps here
    are:

    - normalize the light curve to zero
    - remove giant outliers
    - interpolate gaps in the light curve
      (since ACF requires evenly spaced sampling)

    From McQuillan+ 2013a (https://doi.org/10.1093/mnras/stt536):

        "The ACF calculation requires the light curves to be regularly sampled
        and normalized to zero. We divided the flux in each quarter by its
        median and subtracted unity. Gaps in the light curve longer than the
        Kepler long cadence were filled using linear interpolation with added
        white Gaussian noise. This noise level was estimated using the variance
        of the residuals following subtraction of a smoothed version of the
        flux. To smooth the flux, we applied an iterative non-linear filter
        which consists of a median filter followed by a boxcar filter, both with
        11-point windows, with iterative 3σ clipping of outliers."

    Parameters
    ----------

    times,mags,errs : np.array
        The magnitude/flux time-series and associated measurement errors to
        operate on. Non-finite elements will be removed from these arrays. At
        least 10 elements in each array are required for this function to
        operate.

    fillgaps : {'noiselevel', 'nan'} or float
        If `fillgap='noiselevel'`, fills the gaps with the noise level obtained
        via the procedure above. If `fillgaps='nan'`, fills the gaps with
        `np.nan`. Otherwise, if `fillgaps` is a float, will use that value to
        fill the gaps. The default is to fill the gaps with 0.0 (as in
        McQuillan+ 2014) to "...prevent them contributing to the ACF".

    sigclip : float or int or sequence of two floats/ints or None
        If a single float or int, a symmetric sigma-clip will be performed using
        the number provided as the sigma-multiplier to cut out from the input
        time-series.

        If a list of two ints/floats is provided, the function will perform an
        'asymmetric' sigma-clip. The first element in this list is the sigma
        value to use for fainter flux/mag values; the second element in this
        list is the sigma value to use for brighter flux/mag values. For
        example, `sigclip=[10., 3.]`, will sigclip out greater than 10-sigma
        dimmings and greater than 3-sigma brightenings. Here the meaning of
        "dimming" and "brightening" is set by *physics* (not the magnitude
        system), which is why the `magsarefluxes` kwarg must be correctly set.

        If `sigclip` is None, no sigma-clipping will be performed, and the
        time-series (with non-finite elems removed) will be passed through to
        the output.

    magsareflux : bool
        True if your "mags" are in fact fluxes, i.e. if "fainter" corresponds to
        `mags` getting smaller.

    filterwindow : int
        The number of time-series points to include in the Savitsky-Golay filter
        operation when smoothing the light curve. This should be an odd integer.

    forcetimebin : float or None
        If `forcetimebin` is a float, this value will be used to generate the
        interpolated time series, effectively binning the light curve to this
        cadence. If `forcetimebin` is None, the mode of the gaps (the forward
        difference between successive time values in `times`) in the provided
        light curve will be used as the effective cadence. NOTE: `forcetimebin`
        must be in the same units as `times`, e.g. if times are JD then
        `forcetimebin` must be in days as well

    verbose : bool
        If this is True, will indicate progress at various stages in the
        operation.

    Returns
    -------

    dict
        A dict of the following form is returned::

            {'itimes': the interpolated time values after gap-filling,
             'imags': the interpolated mag/flux values after gap-filling,
             'ierrs': the interpolated mag/flux values after gap-filling,
             'cadence': the cadence of the output mag/flux time-series}