<commit_before># -*- coding: utf-8 -*-
"""
Created on Sun Feb  7 13:43:01 2016

@author: fergal

A series of metrics to quantify the noise in a lightcurve:

Includes:
x sgCdpp
x Marshall's noise estimate
o An FT based estimate of 6 hour artifact strength.
o A per thruster firing estimate of 6 hour artifact strength.

$Id$
$URL$
"""

__version__ = "$Id$"
__URL__ = "$URL$"



from scipy.signal import savgol_filter
import matplotlib.pyplot as mp
import numpy as np
import fft

keplerLongCadence_s = 1765.4679
keplerLongCadence_days = keplerLongCadence_s / float(86400)


def computeRollTweakAmplitude(y, nHarmonics = 3, tweakPeriod_days = .25, \
    expTime_days=None, plot=False):
    """Compute strength of roll tweak artifact in K2 data with an FT approach.

    Compute FT of lightcurve

    Optional Inputs:
    -----------------
    plot
        Show a diagnostic plot

    Returns:
    --------
    float indicating strength of correction. A value of 1 means the
    amplitude of the tweak is approx equal to the strength of all other
    signals in the FT.
    """
    if expTime_days is None:
        expTime_days = keplerLongCadence_days

    #computes FT with frequencies in cycles per days
    ft = fft.computeFft(y, expTime_days)

      #Thruster firings every 6 hours
    artifactFreq_cd = 1/tweakPeriod_days  #cycles per day

    if plot:
        mp.clf()
        mp.plot(ft[:,0], 1e6*ft[:,1], 'b-')

    metric = 0
    nPtsForMed = 50
    for i in range(1, nHarmonics+1):
        wh = np.argmin( np.fabs(ft[:,0] - i*artifactFreq_cd))

        med = np.median(ft[wh-nPtsForMed:wh+nPtsForMed, 1])
        metric += ft[wh, 1] / med

        if plot:
            mp.axvline(i*artifactFreq_cd, color='m')
    return metric / float(nHarmonics)


def computeSgCdpp_ppm(y, transitDuration_cadences=13, plot=False):
    """Estimates 6hr CDPP using <NAME> Cleve's Savitzy-Golay technique

    An interesting estimate of the noise in a lightcurve is the scatter
    after all long term trends have been removed. This is the kernel of
    the idea behind the Combined Differential Photometric Precision (CDPP)
    metric used in classic Kepler. <NAME> devised a much simpler
    algorithm for computing CDPP using a Savitzy-Golay detrending, which
    he called Savitzy-Golay CDPP, or SG-CDPP. We implement his algorithm
    here.


    Inputs:
    ----------
    y
        (1d numpy array) normalised flux to calculate noise from. Flux
        should have a mean of zero and be in units of fractional amplitude.
        Note: Bad data in input will skew result. Some filtering of
        outliers is performed, but Nan's or Infs will not be caught.

    Optional Inputs:
    -----------------
    transitDuration_cadences
        (int) Adjust the assumed transit width, in cadences. Default is
        13, which corresponds to a 6.5 hour transit in K2
    plot
        Show a diagnostic plot

    Returns:
    ------------
    Estimated noise in parts per million.

    Notes:
    -------------
    Taken from
    svn+ssh://murzim/repo/so/trunk/Develop/jvc/common/compute_SG_noise.m
    by <NAME>
    """

    #These 3 values were chosen for the original algorithm, and we don't
    #change them here.
    window = 101
    polyorder=2
    noiseNorm = 1.40
    #Name change for consistency with original algorithm
    cadencesPerTransit = transitDuration_cadences

    if cadencesPerTransit < 4:
        raise ValueError("Cadences per transit must be >= 4")

    trend = savgol_filter(y, window_length=window, polyorder=polyorder)
    detrend = y-trend

    filtered = np.ones(cadencesPerTransit)/float(cadencesPerTransit)
    smoothed = np.convolve(detrend, filtered, mode='same')

    if plot:
        mp.clf()
        mp.plot(y, 'ko')
        mp.plot(trend, 'r-')
        mp.plot(smoothed, 'g.')

    sgCdpp_ppm = noiseNorm*robustStd(smoothed, 1)*1e6
    return sgCdpp_ppm


def estimateScatterWithMarshallMethod(flux, plot=False):
    """Estimate the typical scatter in a lightcurve.

    Uses the same method as Marshall (Mullally et al 2016 submitted)

    Inputs:
    ----------
    flux
        (np 1d array). Flux to measure scatter of. Need not have
        zero mean.<commit_msg>Fixed bug where computeSgCdpp will crash with strange error message if the number of input points in less than the window. Now it gives a saner error message.<commit_after># -*- coding: utf-8 -*-
"""
Created on Sun Feb  7 13:43:01 2016

@author: fergal

A series of metrics to quantify the noise in a lightcurve:

Includes:
x sgCdpp
x Marshall's noise estimate
o An FT based estimate of 6 hour artifact strength.
o A per thruster firing estimate of 6 hour artifact strength.

$Id$
$URL$
"""

__version__ = "$Id$"
__URL__ = "$URL$"



from scipy.signal import savgol_filter
import matplotlib.pyplot as mp
import numpy as np
import fft

keplerLongCadence_s = 1765.4679
keplerLongCadence_days = keplerLongCadence_s / float(86400)


def computeRollTweakAmplitude(y, nHarmonics = 3, tweakPeriod_days = .25, \
    expTime_days=None, plot=False):
    """Compute strength of roll tweak artifact in K2 data with an FT approach.

    Compute FT of lightcurve

    Optional Inputs:
    -----------------
    plot
        Show a diagnostic plot

    Returns:
    --------
    float indicating strength of correction. A value of 1 means the
    amplitude of the tweak is approx equal to the strength of all other
    signals in the FT.
    """
    if expTime_days is None:
        expTime_days = keplerLongCadence_days

    #computes FT with frequencies in cycles per days
    ft = fft.computeFft(y, expTime_days)

      #Thruster firings every 6 hours
    artifactFreq_cd = 1/tweakPeriod_days  #cycles per day

    if plot:
        mp.clf()
        mp.plot(ft[:,0], 1e6*ft[:,1], 'b-')

    metric = 0
    nPtsForMed = 50
    for i in range(1, nHarmonics+1):
        wh = np.argmin( np.fabs(ft[:,0] - i*artifactFreq_cd))

        med = np.median(ft[wh-nPtsForMed:wh+nPtsForMed, 1])
        metric += ft[wh, 1] / med

        if plot:
            mp.axvline(i*artifactFreq_cd, color='m')
    return metric / float(nHarmonics)


def computeSgCdpp_ppm(y, transitDuration_cadences=13, plot=False):
    """Estimates 6hr CDPP using Jeff Van Cleve's Savitzy-Golay technique

    An interesting estimate of the noise in a lightcurve is the scatter
    after all long term trends have been removed. This is the kernel of
    the idea behind the Combined Differential Photometric Precision (CDPP)
    metric used in classic Kepler. <NAME> devised a much simpler
    algorithm for computing CDPP using a Savitzy-Golay detrending, which
    he called Savitzy-Golay CDPP, or SG-CDPP. We implement his algorithm
    here.


    Inputs:
    ----------
    y
        (1d numpy array) normalised flux to calculate noise from. Flux
        should have a mean of zero and be in units of fractional amplitude.
        Note: Bad data in input will skew result. Some filtering of
        outliers is performed, but Nan's or Infs will not be caught.

    Optional Inputs:
    -----------------
    transitDuration_cadences
        (int) Adjust the assumed transit width, in cadences. Default is
        13, which corresponds to a 6.5 hour transit in K2
    plot
        Show a diagnostic plot

    Returns:
    ------------
    Estimated noise in parts per million.

    Notes:
    -------------
    Taken from
    svn+ssh://murzim/repo/so/trunk/Develop/jvc/common/compute_SG_noise.m
    by <NAME>
    """

    #These 3 values were chosen for the original algorithm, and we don't
    #change them here.
    window = 101
    polyorder=2
    noiseNorm = 1.40
    #Name change for consistency with original algorithm
    cadencesPerTransit = transitDuration_cadences

    if cadencesPerTransit < 4:
        raise ValueError("Cadences per transit must be >= 4")

    if len(y) < window:
        raise ValueError("Can't compute CDPP for timeseries with fewer points than defined window (%i points)" %(window))

    trend = savgol_filter(y, window_length=window, polyorder=polyorder)
    detrend = y-trend

    filtered = np.ones(cadencesPerTransit)/float(cadencesPerTransit)
    smoothed = np.convolve(detrend, filtered, mode='same')

    if plot:
        mp.clf()
        mp.plot(y, 'ko')
        mp.plot(trend, 'r-')
        mp.plot(smoothed, 'g.')

    sgCdpp_ppm = noiseNorm*robustStd(smoothed, 1)*1e6
    return sgCdpp_ppm


def estimateScatterWithMarshallMethod(flux, plot=False):
    """Estimate the typical scatter in a lightcurve.

    Uses the same method as Marshall (Mullally et al 2016 submitted)

    Inputs:
    ----------
    flux
        (np 1d array). Flux to measure scatter of. Need not have
        zero mean.
{'fraction_non_alphanumeric': 0.0929295255642561, 'fraction_numerical': 0.016812528788576694, 'mean_word_length': 2.663011387600169, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6560011', 'n_tokens_mistral': 2902, 'n_tokens_neox': 2655, 'n_words': 969}