Source code for astrobase.periodbase.kbls

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# kbls.py - Waqas Bhatti (wbhatti@astro.princeton.edu) - Jan 2017

'''
Contains the Kovacs, et al. (2002) Box-Least-squared-Search period-search
algorithm implementation for periodbase.

'''

#############
## LOGGING ##
#############

import logging
from astrobase import log_sub, log_fmt, log_date_fmt

DEBUG = False
if DEBUG:
    level = logging.DEBUG
else:
    level = logging.INFO
LOGGER = logging.getLogger(__name__)
logging.basicConfig(
    level=level,
    style=log_sub,
    format=log_fmt,
    datefmt=log_date_fmt,
)

LOGDEBUG = LOGGER.debug
LOGINFO = LOGGER.info
LOGWARNING = LOGGER.warning
LOGERROR = LOGGER.error
LOGEXCEPTION = LOGGER.exception


#############
## IMPORTS ##
#############

from math import fmod
from multiprocessing import Pool, cpu_count

from numpy import (
    nan as npnan, arange as nparange, ones as npones, array as nparray,
    isfinite as npisfinite, argmax as npargmax, floor as npfloor,
    linspace as nplinspace, digitize as npdigitize, where as npwhere,
    abs as npabs, min as npmin, full_like as npfull_like, median as npmedian,
    std as npstd, sqrt as npsqrt, ceil as npceil, argsort as npargsort,
    concatenate as npconcatenate, ndarray as npndarray, inf as npinf
)

###################
## LOCAL IMPORTS ##
###################

from pyeebls import eebls

from ..lcmath import sigclip_magseries, phase_magseries_with_errs
from ..lcfit.nonphysical import savgol_fit_magseries
from ..lcfit.transits import traptransit_fit_magseries

from .utils import resort_by_time

############
## CONFIG ##
############

NCPUS = cpu_count()


######################################
## BLS (Kovacs, Zucker, Mazeh 2002) ##
######################################

def _bls_runner(times,
                mags,
                nfreq,
                freqmin,
                stepsize,
                nbins,
                minduration,
                maxduration):
    '''This runs the pyeebls.eebls function using the given inputs.

    Parameters
    ----------

    times,mags : np.array
        The input magnitude time-series to search for transits.

    nfreq : int
        The number of frequencies to use when searching for transits.

    freqmin : float
        The minimum frequency of the period-search -> max period that will be
        used for the search.

    stepsize : float
        The step-size in frequency to use to generate a frequency-grid.

    nbins : int
        The number of phase bins to use.

    minduration : float
        The minimum fractional transit duration that will be considered.

    maxduration : float
        The maximum fractional transit duration that will be considered.

    Returns
    -------

    dict
        Returns a dict of the form::

            {
                'power':           the periodogram power array,
                'bestperiod':      the best period found,
                'bestpower':       the highest peak of the periodogram power,
                'transdepth':      transit depth found by eebls.f,
                'transduration':   transit duration found by eebls.f,
                'transingressbin': transit ingress bin found by eebls.f,
                'transegressbin':  transit egress bin found by eebls.f,
            }

    '''

    workarr_u = npones(times.size)
    workarr_v = npones(times.size)

    blsresult = eebls(times, mags,
                      workarr_u, workarr_v,
                      nfreq, freqmin, stepsize,
                      nbins, minduration, maxduration)

    return {'power':blsresult[0],
            'bestperiod':blsresult[1],
            'bestpower':blsresult[2],
            'transdepth':blsresult[3],
            'transduration':blsresult[4],
            'transingressbin':blsresult[5],
            'transegressbin':blsresult[6]}


def _parallel_bls_worker(task):
    '''
    This wraps the BLS function for the parallel driver below.

    Parameters
    ----------

    tasks : tuple
        This is of the form::

            task[0] = times
            task[1] = mags
            task[2] = nfreq
            task[3] = freqmin
            task[4] = stepsize
            task[5] = nbins
            task[6] = minduration
            task[7] = maxduration

    Returns
    -------

    dict
        Returns a dict of the form::

            {
                'power':           the periodogram power array,
                'bestperiod':      the best period found,
                'bestpower':       the highest peak of the periodogram power,
                'transdepth':      transit depth found by eebls.f,
                'transduration':   transit duration found by eebls.f,
                'transingressbin': transit ingress bin found by eebls.f,
                'transegressbin':  transit egress bin found by eebls.f,
            }

    '''

    try:

        return _bls_runner(*task)

    except Exception:

        LOGEXCEPTION('BLS failed for task %s' % repr(task[2:]))

        return {
            'power':nparray([npnan for x in range(task[2])]),
            'bestperiod':npnan,
            'bestpower':npnan,
            'transdepth':npnan,
            'transduration':npnan,
            'transingressbin':npnan,
            'transegressbin':npnan
        }


[docs]def bls_serial_pfind( times, mags, errs, magsarefluxes=False, startp=0.1, # search from 0.1 d to... endp=100.0, # ... 100.0 d -- don't search full timebase stepsize=5.0e-4, mintransitduration=0.01, # minimum transit length in phase maxtransitduration=0.4, # maximum transit length in phase nphasebins=200, autofreq=True, # figure out f0, nf, and df automatically periodepsilon=0.1, nbestpeaks=5, sigclip=10.0, endp_timebase_check=True, verbose=True, get_stats=True, ): '''Runs the Box Least Squares Fitting Search for transit-shaped signals. Based on eebls.f from Kovacs et al. 2002 and python-bls from Foreman-Mackey et al. 2015. This is the serial version (which is good enough in most cases because BLS in Fortran is fairly fast). If nfreq > 5e5, this will take a while. Parameters ---------- times,mags,errs : np.array The magnitude/flux time-series to search for transits. magsarefluxes : bool If the input measurement values in `mags` and `errs` are in fluxes, set this to True. startp,endp : float The minimum and maximum periods to consider for the transit search. stepsize : float The step-size in frequency to use when constructing a frequency grid for the period search. mintransitduration,maxtransitduration : float The minimum and maximum transitdurations (in units of phase) to consider for the transit search. nphasebins : int The number of phase bins to use in the period search. autofreq : bool If this is True, the values of `stepsize` and `nphasebins` will be ignored, and these, along with a frequency-grid, will be determined based on the following relations:: nphasebins = int(ceil(2.0/mintransitduration)) if nphasebins > 3000: nphasebins = 3000 stepsize = 0.25*mintransitduration/(times.max()-times.min()) minfreq = 1.0/endp maxfreq = 1.0/startp nfreq = int(ceil((maxfreq - minfreq)/stepsize)) periodepsilon : float The fractional difference between successive values of 'best' periods when sorting by periodogram power to consider them as separate periods (as opposed to part of the same periodogram peak). This is used to avoid broad peaks in the periodogram and make sure the 'best' periods returned are all actually independent. nbestpeaks : int The number of 'best' peaks to return from the periodogram results, starting from the global maximum of the periodogram peak values. 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. endp_timebase_check : bool If True, will check if the ``endp`` value is larger than the time-base of the observations. If it is, will change the ``endp`` value such that it is half of the time-base. If False, will allow an ``endp`` larger than the time-base of the observations. verbose : bool If this is True, will indicate progress and details about the frequency grid used for the period search. get_stats : bool If True, runs :py:func:`.bls_stats_singleperiod` for each of the best periods in the output and injects the output into the output dict so you only have to run this function to get the periods and their stats. The output dict from this function will then contain a 'stats' key containing a list of dicts with statistics for each period in ``resultdict['nbestperiods']``. These dicts will contain fit values of transit parameters after a trapezoid transit model is fit to the phased light curve at each period in ``resultdict['nbestperiods']``, i.e. fit values for period, epoch, transit depth, duration, ingress duration, and the SNR of the transit. NOTE: make sure to check the 'fit_status' key for each ``resultdict['stats']`` item to confirm that the trapezoid transit model fit succeeded and that the stats calculated are valid. Returns ------- dict This function returns a dict, referred to as an `lspinfo` dict in other astrobase functions that operate on periodogram results. This is a standardized format across all astrobase period-finders, and is of the form below:: {'bestperiod': the best period value in the periodogram, 'bestlspval': the periodogram peak associated with the best period, 'nbestpeaks': the input value of nbestpeaks, 'nbestlspvals': nbestpeaks-size list of best period peak values, 'nbestperiods': nbestpeaks-size list of best periods, 'stats': BLS stats for each best period, 'lspvals': the full array of periodogram powers, 'frequencies': the full array of frequencies considered, 'periods': the full array of periods considered, 'blsresult': the result dict from the eebls.f wrapper function, 'stepsize': the actual stepsize used, 'nfreq': the actual nfreq used, 'nphasebins': the actual nphasebins used, 'mintransitduration': the input mintransitduration, 'maxtransitduration': the input maxtransitdurations, 'method':'bls' -> the name of the period-finder method, 'kwargs':{ dict of all of the input kwargs for record-keeping}} ''' # get rid of nans first and sigclip stimes, smags, serrs = sigclip_magseries(times, mags, errs, magsarefluxes=magsarefluxes, sigclip=sigclip) # resort by time stimes, smags, errs = resort_by_time(stimes, smags, serrs) # make sure there are enough points to calculate a spectrum if len(stimes) > 9 and len(smags) > 9 and len(serrs) > 9: # if we're setting up everything automatically if autofreq: # figure out the best number of phasebins to use nphasebins = int(npceil(2.0/mintransitduration)) if nphasebins > 3000: nphasebins = 3000 # use heuristic to figure out best timestep stepsize = 0.25*mintransitduration/(stimes.max()-stimes.min()) # now figure out the frequencies to use minfreq = 1.0/endp maxfreq = 1.0/startp nfreq = int(npceil((maxfreq - minfreq)/stepsize)) # say what we're using if verbose: LOGINFO('min P: %s, max P: %s, nfreq: %s, ' 'minfreq: %s, maxfreq: %s' % (startp, endp, nfreq, minfreq, maxfreq)) LOGINFO('autofreq = True: using AUTOMATIC values for ' 'freq stepsize: %s, nphasebins: %s, ' 'min transit duration: %s, max transit duration: %s' % (stepsize, nphasebins, mintransitduration, maxtransitduration)) else: minfreq = 1.0/endp maxfreq = 1.0/startp nfreq = int(npceil((maxfreq - minfreq)/stepsize)) # say what we're using if verbose: LOGINFO('min P: %s, max P: %s, nfreq: %s, ' 'minfreq: %s, maxfreq: %s' % (startp, endp, nfreq, minfreq, maxfreq)) LOGINFO('autofreq = False: using PROVIDED values for ' 'freq stepsize: %s, nphasebins: %s, ' 'min transit duration: %s, max transit duration: %s' % (stepsize, nphasebins, mintransitduration, maxtransitduration)) if nfreq > 5.0e5: if verbose: LOGWARNING('more than 5.0e5 frequencies to go through; ' 'this will take a while. ' 'you might want to use the ' 'periodbase.bls_parallel_pfind function instead') if ((minfreq < (1.0/(stimes.max() - stimes.min()))) and endp_timebase_check): LOGWARNING('the requested max P = %.3f is larger than ' 'the time base of the observations = %.3f, ' ' will make minfreq = 2 x 1/timebase' % (endp, stimes.max() - stimes.min())) minfreq = 2.0/(stimes.max() - stimes.min()) LOGWARNING('new minfreq: %s, maxfreq: %s' % (minfreq, maxfreq)) # # run BLS # try: blsresult = _bls_runner(stimes, smags, nfreq, minfreq, stepsize, nphasebins, mintransitduration, maxtransitduration) frequencies = minfreq + nparange(nfreq)*stepsize periods = 1.0/frequencies lsp = blsresult['power'] # find the nbestpeaks for the periodogram: 1. sort the lsp array # by highest value first 2. go down the values until we find # five values that are separated by at least periodepsilon in # period # make sure to get only the finite peaks in the periodogram # this is needed because BLS may produce infs for some peaks finitepeakind = npisfinite(lsp) finlsp = lsp[finitepeakind] finperiods = periods[finitepeakind] # make sure that finlsp has finite values before we work on it try: bestperiodind = npargmax(finlsp) except ValueError: LOGERROR('no finite periodogram values ' 'for this mag series, skipping...') return {'bestperiod':npnan, 'bestlspval':npnan, 'nbestpeaks':nbestpeaks, 'nbestlspvals':None, 'nbestperiods':None, 'lspvals':None, 'periods':None, 'method':'bls', 'kwargs':{'startp':startp, 'endp':endp, 'stepsize':stepsize, 'mintransitduration':mintransitduration, 'maxtransitduration':maxtransitduration, 'nphasebins':nphasebins, 'autofreq':autofreq, 'periodepsilon':periodepsilon, 'nbestpeaks':nbestpeaks, 'sigclip':sigclip, 'magsarefluxes':magsarefluxes}} sortedlspind = npargsort(finlsp)[::-1] sortedlspperiods = finperiods[sortedlspind] sortedlspvals = finlsp[sortedlspind] # now get the nbestpeaks nbestperiods, nbestlspvals, peakcount = ( [finperiods[bestperiodind]], [finlsp[bestperiodind]], 1 ) prevperiod = sortedlspperiods[0] # find the best nbestpeaks in the lsp and their periods for period, lspval in zip(sortedlspperiods, sortedlspvals): if peakcount == nbestpeaks: break perioddiff = abs(period - prevperiod) bestperiodsdiff = [abs(period - x) for x in nbestperiods] # this ensures that this period is different from the last # period and from all the other existing best periods by # periodepsilon to make sure we jump to an entire different # peak in the periodogram if (perioddiff > (periodepsilon*prevperiod) and all(x > (periodepsilon*period) for x in bestperiodsdiff)): nbestperiods.append(period) nbestlspvals.append(lspval) peakcount = peakcount + 1 prevperiod = period # generate the return dict resultdict = { 'bestperiod':finperiods[bestperiodind], 'bestlspval':finlsp[bestperiodind], 'nbestpeaks':nbestpeaks, 'nbestlspvals':nbestlspvals, 'nbestperiods':nbestperiods, 'lspvals':lsp, 'frequencies':frequencies, 'periods':periods, 'blsresult':blsresult, 'stepsize':stepsize, 'nfreq':nfreq, 'nphasebins':nphasebins, 'mintransitduration':mintransitduration, 'maxtransitduration':maxtransitduration, 'method':'bls', 'kwargs':{'startp':startp, 'endp':endp, 'stepsize':stepsize, 'mintransitduration':mintransitduration, 'maxtransitduration':maxtransitduration, 'nphasebins':nphasebins, 'autofreq':autofreq, 'periodepsilon':periodepsilon, 'nbestpeaks':nbestpeaks, 'sigclip':sigclip, 'magsarefluxes':magsarefluxes} } # get stats if requested if get_stats: resultdict['stats'] = [] for bp in nbestperiods: if verbose: LOGINFO("Getting stats for best period: %.6f" % bp) this_pstats = bls_stats_singleperiod( stimes, smags, serrs, bp, magsarefluxes=resultdict['kwargs']['magsarefluxes'], sigclip=resultdict['kwargs']['sigclip'], nphasebins=resultdict['nphasebins'], mintransitduration=resultdict['mintransitduration'], maxtransitduration=resultdict['maxtransitduration'], verbose=verbose, ) resultdict['stats'].append(this_pstats) return resultdict except Exception: LOGEXCEPTION('BLS failed!') return {'bestperiod':npnan, 'bestlspval':npnan, 'nbestpeaks':nbestpeaks, 'nbestlspvals':None, 'nbestperiods':None, 'lspvals':None, 'periods':None, 'blsresult':None, 'stepsize':stepsize, 'nfreq':nfreq, 'nphasebins':nphasebins, 'mintransitduration':mintransitduration, 'maxtransitduration':maxtransitduration, 'method':'bls', 'kwargs':{'startp':startp, 'endp':endp, 'stepsize':stepsize, 'mintransitduration':mintransitduration, 'maxtransitduration':maxtransitduration, 'nphasebins':nphasebins, 'autofreq':autofreq, 'periodepsilon':periodepsilon, 'nbestpeaks':nbestpeaks, 'sigclip':sigclip, 'magsarefluxes':magsarefluxes}} else: LOGERROR('no good detections for these times and mags, skipping...') return {'bestperiod':npnan, 'bestlspval':npnan, 'nbestpeaks':nbestpeaks, 'nbestlspvals':None, 'nbestperiods':None, 'lspvals':None, 'periods':None, 'blsresult':None, 'stepsize':stepsize, 'nfreq':None, 'nphasebins':None, 'mintransitduration':mintransitduration, 'maxtransitduration':maxtransitduration, 'method':'bls', 'kwargs':{'startp':startp, 'endp':endp, 'stepsize':stepsize, 'mintransitduration':mintransitduration, 'maxtransitduration':maxtransitduration, 'nphasebins':nphasebins, 'autofreq':autofreq, 'periodepsilon':periodepsilon, 'nbestpeaks':nbestpeaks, 'sigclip':sigclip, 'magsarefluxes':magsarefluxes}}
[docs]def bls_parallel_pfind( times, mags, errs, magsarefluxes=False, startp=0.1, # by default, search from 0.1 d to... endp=100.0, # ... 100.0 d -- don't search full timebase stepsize=1.0e-4, mintransitduration=0.01, # minimum transit length in phase maxtransitduration=0.4, # maximum transit length in phase nphasebins=200, autofreq=True, # figure out f0, nf, and df automatically nbestpeaks=5, periodepsilon=0.1, # 0.1 sigclip=10.0, endp_timebase_check=True, verbose=True, nworkers=None, get_stats=True, ): '''Runs the Box Least Squares Fitting Search for transit-shaped signals. Based on eebls.f from Kovacs et al. 2002 and python-bls from Foreman-Mackey et al. 2015. Breaks up the full frequency space into chunks and passes them to parallel BLS workers. NOTE: the combined BLS spectrum produced by this function is not identical to that produced by running BLS in one shot for the entire frequency space. There are differences on the order of 1.0e-3 or so in the respective peak values, but peaks appear at the same frequencies for both methods. This is likely due to different aliasing caused by smaller chunks of the frequency space used by the parallel workers in this function. When in doubt, confirm results for this parallel implementation by comparing to those from the serial implementation above. Parameters ---------- times,mags,errs : np.array The magnitude/flux time-series to search for transits. magsarefluxes : bool If the input measurement values in `mags` and `errs` are in fluxes, set this to True. startp,endp : float The minimum and maximum periods to consider for the transit search. stepsize : float The step-size in frequency to use when constructing a frequency grid for the period search. mintransitduration,maxtransitduration : float The minimum and maximum transitdurations (in units of phase) to consider for the transit search. nphasebins : int The number of phase bins to use in the period search. autofreq : bool If this is True, the values of `stepsize` and `nphasebins` will be ignored, and these, along with a frequency-grid, will be determined based on the following relations:: nphasebins = int(ceil(2.0/mintransitduration)) if nphasebins > 3000: nphasebins = 3000 stepsize = 0.25*mintransitduration/(times.max()-times.min()) minfreq = 1.0/endp maxfreq = 1.0/startp nfreq = int(ceil((maxfreq - minfreq)/stepsize)) periodepsilon : float The fractional difference between successive values of 'best' periods when sorting by periodogram power to consider them as separate periods (as opposed to part of the same periodogram peak). This is used to avoid broad peaks in the periodogram and make sure the 'best' periods returned are all actually independent. nbestpeaks : int The number of 'best' peaks to return from the periodogram results, starting from the global maximum of the periodogram peak values. 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. endp_timebase_check : bool If True, will check if the ``endp`` value is larger than the time-base of the observations. If it is, will change the ``endp`` value such that it is half of the time-base. If False, will allow an ``endp`` larger than the time-base of the observations. verbose : bool If this is True, will indicate progress and details about the frequency grid used for the period search. nworkers : int or None The number of parallel workers to launch for period-search. If None, nworkers = NCPUS. get_stats : bool If True, runs :py:func:`.bls_stats_singleperiod` for each of the best periods in the output and injects the output into the output dict so you only have to run this function to get the periods and their stats. The output dict from this function will then contain a 'stats' key containing a list of dicts with statistics for each period in ``resultdict['nbestperiods']``. These dicts will contain fit values of transit parameters after a trapezoid transit model is fit to the phased light curve at each period in ``resultdict['nbestperiods']``, i.e. fit values for period, epoch, transit depth, duration, ingress duration, and the SNR of the transit. NOTE: make sure to check the 'fit_status' key for each ``resultdict['stats']`` item to confirm that the trapezoid transit model fit succeeded and that the stats calculated are valid. Returns ------- dict This function returns a dict, referred to as an `lspinfo` dict in other astrobase functions that operate on periodogram results. This is a standardized format across all astrobase period-finders, and is of the form below:: {'bestperiod': the best period value in the periodogram, 'bestlspval': the periodogram peak associated with the best period, 'nbestpeaks': the input value of nbestpeaks, 'nbestlspvals': nbestpeaks-size list of best period peak values, 'nbestperiods': nbestpeaks-size list of best periods, 'stats': list of stats dicts returned for each best period, 'lspvals': the full array of periodogram powers, 'frequencies': the full array of frequencies considered, 'periods': the full array of periods considered, 'blsresult': list of result dicts from eebls.f wrapper functions, 'stepsize': the actual stepsize used, 'nfreq': the actual nfreq used, 'nphasebins': the actual nphasebins used, 'mintransitduration': the input mintransitduration, 'maxtransitduration': the input maxtransitdurations, 'method':'bls' -> the name of the period-finder method, 'kwargs':{ dict of all of the input kwargs for record-keeping}} ''' # get rid of nans first and sigclip stimes, smags, serrs = sigclip_magseries(times, mags, errs, magsarefluxes=magsarefluxes, sigclip=sigclip) # resort by time stimes, smags, errs = resort_by_time(stimes, smags, serrs) # make sure there are enough points to calculate a spectrum if len(stimes) > 9 and len(smags) > 9 and len(serrs) > 9: # if we're setting up everything automatically if autofreq: # figure out the best number of phasebins to use nphasebins = int(npceil(2.0/mintransitduration)) if nphasebins > 3000: nphasebins = 3000 # use heuristic to figure out best timestep stepsize = 0.25*mintransitduration/(stimes.max()-stimes.min()) # now figure out the frequencies to use minfreq = 1.0/endp maxfreq = 1.0/startp nfreq = int(npceil((maxfreq - minfreq)/stepsize)) # say what we're using if verbose: LOGINFO('min P: %s, max P: %s, nfreq: %s, ' 'minfreq: %s, maxfreq: %s' % (startp, endp, nfreq, minfreq, maxfreq)) LOGINFO('autofreq = True: using AUTOMATIC values for ' 'freq stepsize: %s, nphasebins: %s, ' 'min transit duration: %s, max transit duration: %s' % (stepsize, nphasebins, mintransitduration, maxtransitduration)) else: minfreq = 1.0/endp maxfreq = 1.0/startp nfreq = int(npceil((maxfreq - minfreq)/stepsize)) # say what we're using if verbose: LOGINFO('min P: %s, max P: %s, nfreq: %s, ' 'minfreq: %s, maxfreq: %s' % (startp, endp, nfreq, minfreq, maxfreq)) LOGINFO('autofreq = False: using PROVIDED values for ' 'freq stepsize: %s, nphasebins: %s, ' 'min transit duration: %s, max transit duration: %s' % (stepsize, nphasebins, mintransitduration, maxtransitduration)) # check the minimum frequency if ((minfreq < (1.0/(stimes.max() - stimes.min()))) and endp_timebase_check): LOGWARNING('the requested max P = %.3f is larger than ' 'the time base of the observations = %.3f, ' ' will make minfreq = 2 x 1/timebase' % (endp, stimes.max() - stimes.min())) minfreq = 2.0/(stimes.max() - stimes.min()) LOGWARNING('new minfreq: %s, maxfreq: %s' % (minfreq, maxfreq)) ############################# ## NOW RUN BLS IN PARALLEL ## ############################# # fix number of CPUs if needed if not nworkers or nworkers > NCPUS: nworkers = NCPUS if verbose: LOGINFO('using %s workers...' % nworkers) # the frequencies array to be searched frequencies = minfreq + nparange(nfreq)*stepsize # break up the tasks into chunks csrem = int(fmod(nfreq, nworkers)) csint = int(float(nfreq/nworkers)) chunk_minfreqs, chunk_nfreqs = [], [] for x in range(nworkers): this_minfreqs = frequencies[x*csint] # handle usual nfreqs if x < (nworkers - 1): this_nfreqs = frequencies[x*csint:x*csint+csint].size else: this_nfreqs = frequencies[x*csint:x*csint+csint+csrem].size chunk_minfreqs.append(this_minfreqs) chunk_nfreqs.append(this_nfreqs) # populate the tasks list tasks = [(stimes, smags, chunk_minf, chunk_nf, stepsize, nphasebins, mintransitduration, maxtransitduration) for (chunk_nf, chunk_minf) in zip(chunk_minfreqs, chunk_nfreqs)] if verbose: for ind, task in enumerate(tasks): LOGINFO('worker %s: minfreq = %.6f, nfreqs = %s' % (ind+1, task[3], task[2])) LOGINFO('running...') # return tasks # start the pool pool = Pool(nworkers) results = pool.map(_parallel_bls_worker, tasks) pool.close() pool.join() del pool # now concatenate the output lsp arrays lsp = npconcatenate([x['power'] for x in results]) periods = 1.0/frequencies # find the nbestpeaks for the periodogram: 1. sort the lsp array # by highest value first 2. go down the values until we find # five values that are separated by at least periodepsilon in # period # make sure to get only the finite peaks in the periodogram # this is needed because BLS may produce infs for some peaks finitepeakind = npisfinite(lsp) finlsp = lsp[finitepeakind] finperiods = periods[finitepeakind] # make sure that finlsp has finite values before we work on it try: bestperiodind = npargmax(finlsp) except ValueError: LOGERROR('no finite periodogram values ' 'for this mag series, skipping...') return {'bestperiod':npnan, 'bestlspval':npnan, 'nbestpeaks':nbestpeaks, 'nbestlspvals':None, 'nbestperiods':None, 'lspvals':None, 'periods':None, 'blsresult':None, 'method':'bls', 'kwargs':{'startp':startp, 'endp':endp, 'stepsize':stepsize, 'mintransitduration':mintransitduration, 'maxtransitduration':maxtransitduration, 'nphasebins':nphasebins, 'autofreq':autofreq, 'periodepsilon':periodepsilon, 'nbestpeaks':nbestpeaks, 'sigclip':sigclip, 'magsarefluxes':magsarefluxes}} sortedlspind = npargsort(finlsp)[::-1] sortedlspperiods = finperiods[sortedlspind] sortedlspvals = finlsp[sortedlspind] # now get the nbestpeaks nbestperiods, nbestlspvals, peakcount = ( [finperiods[bestperiodind]], [finlsp[bestperiodind]], 1 ) prevperiod = sortedlspperiods[0] # find the best nbestpeaks in the lsp and their periods for period, lspval in zip(sortedlspperiods, sortedlspvals): if peakcount == nbestpeaks: break perioddiff = abs(period - prevperiod) bestperiodsdiff = [abs(period - x) for x in nbestperiods] # this ensures that this period is different from the last # period and from all the other existing best periods by # periodepsilon to make sure we jump to an entire different # peak in the periodogram if (perioddiff > (periodepsilon*prevperiod) and all(x > (periodepsilon*period) for x in bestperiodsdiff)): nbestperiods.append(period) nbestlspvals.append(lspval) peakcount = peakcount + 1 prevperiod = period # generate the return dict resultdict = { 'bestperiod':finperiods[bestperiodind], 'bestlspval':finlsp[bestperiodind], 'nbestpeaks':nbestpeaks, 'nbestlspvals':nbestlspvals, 'nbestperiods':nbestperiods, 'lspvals':lsp, 'frequencies':frequencies, 'periods':periods, 'blsresult':results, 'stepsize':stepsize, 'nfreq':nfreq, 'nphasebins':nphasebins, 'mintransitduration':mintransitduration, 'maxtransitduration':maxtransitduration, 'method':'bls', 'kwargs':{'startp':startp, 'endp':endp, 'stepsize':stepsize, 'mintransitduration':mintransitduration, 'maxtransitduration':maxtransitduration, 'nphasebins':nphasebins, 'autofreq':autofreq, 'periodepsilon':periodepsilon, 'nbestpeaks':nbestpeaks, 'sigclip':sigclip, 'magsarefluxes':magsarefluxes} } # get stats if requested if get_stats: resultdict['stats'] = [] for bp in nbestperiods.copy(): if verbose: LOGINFO("Getting stats for best period: %.6f" % bp) this_pstats = bls_stats_singleperiod( stimes, smags, serrs, bp, magsarefluxes=resultdict['kwargs']['magsarefluxes'], sigclip=resultdict['kwargs']['sigclip'], nphasebins=resultdict['nphasebins'], mintransitduration=resultdict['mintransitduration'], maxtransitduration=resultdict['maxtransitduration'], verbose=verbose, ) resultdict['stats'].append(this_pstats) return resultdict else: LOGERROR('no good detections for these times and mags, skipping...') return {'bestperiod':npnan, 'bestlspval':npnan, 'nbestpeaks':nbestpeaks, 'nbestlspvals':None, 'nbestperiods':None, 'lspvals':None, 'periods':None, 'blsresult':None, 'stepsize':stepsize, 'nfreq':None, 'nphasebins':None, 'mintransitduration':mintransitduration, 'maxtransitduration':maxtransitduration, 'method':'bls', 'kwargs':{'startp':startp, 'endp':endp, 'stepsize':stepsize, 'mintransitduration':mintransitduration, 'maxtransitduration':maxtransitduration, 'nphasebins':nphasebins, 'autofreq':autofreq, 'periodepsilon':periodepsilon, 'nbestpeaks':nbestpeaks, 'sigclip':sigclip, 'magsarefluxes':magsarefluxes}}
def _get_bls_stats(stimes, smags, serrs, thistransdepth, thistransduration, ingressdurationfraction, nphasebins, thistransingressbin, thistransegressbin, thisbestperiod, thisnphasebins, magsarefluxes=False, verbose=False): ''' Actually calculates the stats. ''' try: # try getting the minimum light epoch using the phase bin method me_epochbin = int((thistransegressbin + thistransingressbin)/2.0) me_phases = ( (stimes - stimes.min())/thisbestperiod - npfloor((stimes - stimes.min())/thisbestperiod) ) me_phases_sortind = npargsort(me_phases) me_sorted_phases = me_phases[me_phases_sortind] me_sorted_times = stimes[me_phases_sortind] me_bins = nplinspace(0.0, 1.0, thisnphasebins) me_bininds = npdigitize(me_sorted_phases, me_bins) me_centertransit_ind = me_bininds == me_epochbin me_centertransit_phase = ( npmedian(me_sorted_phases[me_centertransit_ind]) ) me_centertransit_timeloc = npwhere( npabs(me_sorted_phases - me_centertransit_phase) == npmin(npabs(me_sorted_phases - me_centertransit_phase)) ) me_centertransit_time = me_sorted_times[ me_centertransit_timeloc ] if me_centertransit_time.size > 1: LOGWARNING('multiple possible times-of-center transits ' 'found for period %.7f, picking the first ' 'one from: %s' % (thisbestperiod, repr(me_centertransit_time))) thisminepoch = me_centertransit_time[0] except Exception: LOGEXCEPTION( 'could not determine the center time of transit for ' 'the phased LC, trying SavGol fit instead...' ) # fit a Savitsky-Golay instead and get its minimum savfit = savgol_fit_magseries(stimes, smags, serrs, thisbestperiod, magsarefluxes=magsarefluxes, verbose=verbose, sigclip=None) thisminepoch = savfit['fitinfo']['fitepoch'] if isinstance(thisminepoch, npndarray): if verbose: LOGWARNING('minimum epoch is actually an array:\n' '%s\n' 'instead of a float, ' 'are there duplicate time values ' 'in the original input? ' 'will use the first value in this array.' % repr(thisminepoch)) thisminepoch = thisminepoch[0] # make sure the INITIAL value of the ingress duration isn't more than half # of the total duration of the transit model_ingressduration = ingressdurationfraction*thistransduration if model_ingressduration > (0.5*thistransduration): model_ingressduration = 0.5*thistransduration # # set the depth bounds appropriately for the type of mag input # # require positive depth for fluxes if magsarefluxes: transit_depth_bounds = (0.0, npinf) # require negative depth for mags else: transit_depth_bounds = (-npinf, 0.0) # set up trapezoid transit model to fit for this LC transitparams = [ thisbestperiod, thisminepoch, thistransdepth, thistransduration, model_ingressduration ] # # run the model fit # modelfit = traptransit_fit_magseries( stimes, smags, serrs, transitparams, sigclip=None, magsarefluxes=magsarefluxes, param_bounds={ 'period':(0.0, npinf), 'epoch':(0.0, npinf), 'depth':transit_depth_bounds, 'duration':(0.0,1.0), # allow ingress duration during fitting to float up to 0.5 # FIXME: this should technically always be a fn of duration but # curve_fit doesn't support this kind of constraint # use scipy.optimize.minimize instead? 'ingressduration':(0.0,0.5), }, verbose=verbose ) # if the model fit succeeds, calculate SNR using the trapezoid model fit if modelfit and modelfit['fitinfo']['finalparams'] is not None: fitparams = modelfit['fitinfo']['finalparams'] fiterrs = modelfit['fitinfo']['finalparamerrs'] modelmags, actualmags, modelphase = ( modelfit['fitinfo']['fitmags'], modelfit['magseries']['mags'], modelfit['magseries']['phase'] ) subtractedmags = actualmags - modelmags subtractedrms = npstd(subtractedmags) fit_period, fit_epoch, fit_depth, fit_duration, fit_ingress_dur = ( fitparams ) npts_in_transit = modelfit['fitinfo']['ntransitpoints'] transit_snr = ( npsqrt(npts_in_transit) * npabs(fit_depth/subtractedrms) ) if verbose: LOGINFO('refit best period: %.6f, ' 'refit center of transit: %.5f' % (fit_period, fit_epoch)) LOGINFO('npoints in transit: %s' % npts_in_transit) LOGINFO('transit depth (delta): %.5f, ' 'frac transit length (q): %.3f, ' ' SNR: %.3f' % (fit_depth, fit_duration, transit_snr)) return {'period':fit_period, 'epoch':fit_epoch, 'snr':transit_snr, 'transitdepth':fit_depth, 'transitduration':fit_duration, 'ingressduration':fit_ingress_dur, 'npoints_in_transit':npts_in_transit, 'fitparams':fitparams, 'fiterrs':fiterrs, 'fit_status':'ok', 'nphasebins':nphasebins, 'transingressbin':thistransingressbin, 'transegressbin':thistransegressbin, 'blsmodel':modelmags, 'subtractedmags':subtractedmags, 'phasedmags':actualmags, 'phases':modelphase, 'fitinfo':modelfit} # if the model fit doesn't work, then do the SNR calculation the old way else: # phase using this epoch phased_magseries = phase_magseries_with_errs(stimes, smags, serrs, thisbestperiod, thisminepoch, wrap=False, sort=True) tphase = phased_magseries['phase'] tmags = phased_magseries['mags'] # use the transit depth and duration to subtract the BLS transit # model from the phased mag series. we're centered about 0.0 as the # phase of the transit minimum so we need to look at stuff from # [0.0, transitphase] and [1.0-transitphase, 1.0] transitphase = thistransduration/2.0 transitindices = ((tphase < transitphase) | (tphase > (1.0 - transitphase))) # this is the BLS model # constant = median(tmags) outside transit # constant = thistransitdepth inside transit blsmodel = npfull_like(tmags, npmedian(tmags)) if magsarefluxes: # eebls.f returns +ve transit depth for fluxes # so we need to subtract here to get fainter fluxes in transit blsmodel[transitindices] = ( blsmodel[transitindices] - thistransdepth ) else: # eebls.f returns -ve transit depth for magnitudes # so we need to subtract here to get fainter mags in transits blsmodel[transitindices] = ( blsmodel[transitindices] - thistransdepth ) # see __init__/get_snr_of_dip docstring for description of transit # SNR equation, which is what we use for `thissnr`. subtractedmags = tmags - blsmodel subtractedrms = npstd(subtractedmags) npts_in_transit = len(tmags[transitindices]) thissnr = ( npsqrt(npts_in_transit) * npabs(thistransdepth/subtractedrms) ) # tell user about stuff if verbose = True if verbose: LOGINFO('refit best period: %.6f, ' 'refit center of transit: %.5f' % (thisbestperiod, thisminepoch)) LOGINFO('transit ingress phase = %.3f to %.3f' % (1.0 - transitphase, 1.0)) LOGINFO('transit egress phase = %.3f to %.3f' % (0.0, transitphase)) LOGINFO('npoints in transit: %s' % tmags[transitindices].size) LOGINFO('transit depth (delta): %.5f, ' 'frac transit length (q): %.3f, ' ' SNR: %.3f' % (thistransdepth, thistransduration, thissnr)) return {'period':thisbestperiod, 'epoch':thisminepoch, 'snr':thissnr, 'transitdepth':thistransdepth, 'transitduration':thistransduration, 'npoints_in_transit':npts_in_transit, 'fit_status':'trapezoid model fit failed, using box model', 'nphasebins':nphasebins, 'transingressbin':thistransingressbin, 'transegressbin':thistransegressbin, 'blsmodel':blsmodel, 'subtractedmags':subtractedmags, 'phasedmags':tmags, 'phases':tphase}
[docs]def bls_stats_singleperiod(times, mags, errs, period, magsarefluxes=False, sigclip=10.0, perioddeltapercent=10, nphasebins=200, mintransitduration=0.01, maxtransitduration=0.4, ingressdurationfraction=0.1, verbose=True): '''This calculates the SNR, depth, duration, a refit period, and time of center-transit for a single period. The equation used for SNR is:: SNR = (transit model depth / RMS of LC with transit model subtracted) * sqrt(number of points in transit) NOTE: you should set the kwargs `sigclip`, `nphasebins`, `mintransitduration`, `maxtransitduration` to what you used for an initial BLS run to detect transits in the input light curve to match those input conditions. Parameters ---------- times,mags,errs : np.array These contain the magnitude/flux time-series and any associated errors. period : float The period to search around and refit the transits. This will be used to calculate the start and end periods of a rerun of BLS to calculate the stats. magsarefluxes : bool Set to True if the input measurements in `mags` are actually fluxes and not magnitudes. 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. perioddeltapercent : float The fraction of the period provided to use to search around this value. This is a percentage. The period range searched will then be:: [period - (perioddeltapercent/100.0)*period, period + (perioddeltapercent/100.0)*period] nphasebins : int The number of phase bins to use in the BLS run. mintransitduration : float The minimum transit duration in phase to consider. maxtransitduration : float The maximum transit duration to consider. ingressdurationfraction : float The fraction of the transit duration to use to generate an initial value of the transit ingress duration for the BLS model refit. This will be fit by this function. verbose : bool If True, will indicate progress and any problems encountered. Returns ------- dict A dict of the following form is returned:: {'period': the refit best period, 'epoch': the refit epoch (i.e. mid-transit time), 'snr':the SNR of the transit, 'transitdepth':the depth of the transit, 'transitduration':the duration of the transit, 'ingressduration':if trapezoid fit OK, is the ingress duration, 'npoints_in_transit':the number of LC points in transit, 'fit_status': 'ok' or 'trapezoid model fit failed,...', 'nphasebins':the input value of nphasebins, 'transingressbin':the phase bin containing transit ingress, 'transegressbin':the phase bin containing transit egress, 'blsmodel':the full BLS model used along with its parameters, 'subtractedmags':BLS model - phased light curve, 'phasedmags':the phase light curve, 'phases': the phase values} You should check the 'fit_status' key in this returned dict for a value of 'ok'. If it is 'trapezoid model fit failed, using box model', you may not want to trust the transit period and epoch found. ''' # get rid of nans first and sigclip stimes, smags, serrs = sigclip_magseries(times, mags, errs, magsarefluxes=magsarefluxes, sigclip=sigclip) # make sure there are enough points to calculate a spectrum if len(stimes) > 9 and len(smags) > 9 and len(serrs) > 9: # get the period interval startp = period - perioddeltapercent*period/100.0 if startp < 0: startp = period endp = period + perioddeltapercent*period/100.0 # rerun BLS in serial mode around the specified period to get the # transit depth, duration, ingress and egress bins blsres = bls_serial_pfind(stimes, smags, serrs, verbose=verbose, startp=startp, endp=endp, nphasebins=nphasebins, mintransitduration=mintransitduration, maxtransitduration=maxtransitduration, magsarefluxes=magsarefluxes, get_stats=False, sigclip=None) if (not blsres or 'blsresult' not in blsres or blsres['blsresult'] is None): LOGERROR("BLS failed during a period-search " "performed around the input best period: %.6f. " "Can't continue. " % period) return None thistransdepth = blsres['blsresult']['transdepth'] thistransduration = blsres['blsresult']['transduration'] thisbestperiod = blsres['bestperiod'] thistransingressbin = blsres['blsresult']['transingressbin'] thistransegressbin = blsres['blsresult']['transegressbin'] thisnphasebins = nphasebins stats = _get_bls_stats(stimes, smags, serrs, thistransdepth, thistransduration, ingressdurationfraction, nphasebins, thistransingressbin, thistransegressbin, thisbestperiod, thisnphasebins, magsarefluxes=magsarefluxes, verbose=verbose) return stats # if there aren't enough points in the mag series, bail out else: LOGERROR('no good detections for these times and mags, skipping...') return None
[docs]def bls_snr(blsdict, times, mags, errs, assumeserialbls=False, magsarefluxes=False, sigclip=10.0, npeaks=None, perioddeltapercent=10, ingressdurationfraction=0.1, verbose=True): '''Calculates the signal to noise ratio for each best peak in the BLS periodogram, along with transit depth, duration, and refit period and epoch. The following equation is used for SNR:: SNR = (transit model depth / RMS of LC with transit model subtracted) * sqrt(number of points in transit) Parameters ---------- blsdict : dict This is an lspinfo dict produced by either `bls_parallel_pfind` or `bls_serial_pfind` in this module, or by your own BLS function. If you provide results in a dict from an external BLS function, make sure this matches the form below:: {'bestperiod': the best period value in the periodogram, 'bestlspval': the periodogram peak associated with the best period, 'nbestpeaks': the input value of nbestpeaks, 'nbestlspvals': nbestpeaks-size list of best period peak values, 'nbestperiods': nbestpeaks-size list of best periods, 'lspvals': the full array of periodogram powers, 'frequencies': the full array of frequencies considered, 'periods': the full array of periods considered, 'blsresult': list of result dicts from eebls.f wrapper functions, 'stepsize': the actual stepsize used, 'nfreq': the actual nfreq used, 'nphasebins': the actual nphasebins used, 'mintransitduration': the input mintransitduration, 'maxtransitduration': the input maxtransitdurations, 'method':'bls' -> the name of the period-finder method, 'kwargs':{ dict of all of the input kwargs for record-keeping}} times,mags,errs : np.array These contain the magnitude/flux time-series and any associated errors. assumeserialbls : bool If this is True, this function will not rerun BLS around each best peak in the input lspinfo dict to refit the periods and epochs. This is usally required for `bls_parallel_pfind` so set this to False if you use results from that function. The parallel method breaks up the frequency space into chunks for speed, and the results may not exactly match those from a regular BLS run. magsarefluxes : bool Set to True if the input measurements in `mags` are actually fluxes and not magnitudes. npeaks : int or None This controls how many of the periods in `blsdict['nbestperiods']` to find the SNR for. If it's None, then this will calculate the SNR for all of them. If it's an integer between 1 and `len(blsdict['nbestperiods'])`, will calculate for only the specified number of peak periods, starting from the best period. perioddeltapercent : float The fraction of the period provided to use to search around this value. This is a percentage. The period range searched will then be:: [period - (perioddeltapercent/100.0)*period, period + (perioddeltapercent/100.0)*period] ingressdurationfraction : float The fraction of the transit duration to use to generate an initial value of the transit ingress duration for the BLS model refit. This will be fit by this function. verbose : bool If True, will indicate progress and any problems encountered. Returns ------- dict A dict of the following form is returned:: {'npeaks: the number of periodogram peaks requested to get SNR for, 'period': list of refit best periods for each requested peak, 'epoch': list of refit epochs (i.e. mid-transit times), 'snr':list of SNRs of the transit for each requested peak, 'transitdepth':list of depths of the transits, 'transitduration':list of durations of the transits, 'nphasebins':the input value of nphasebins, 'transingressbin':the phase bin containing transit ingress, 'transegressbin':the phase bin containing transit egress, 'allblsmodels':the full BLS models used along with its parameters, 'allsubtractedmags':BLS models - phased light curves, 'allphasedmags':the phase light curves, 'allphases': the phase values} ''' # figure out how many periods to work on if (npeaks and (0 < npeaks < len(blsdict['nbestperiods']))): nperiods = npeaks else: if verbose: LOGWARNING('npeaks not specified or invalid, ' 'getting SNR for all %s BLS peaks' % len(blsdict['nbestperiods'])) nperiods = len(blsdict['nbestperiods']) nbestperiods = blsdict['nbestperiods'][:nperiods] # get rid of nans first and sigclip stimes, smags, serrs = sigclip_magseries(times, mags, errs, magsarefluxes=magsarefluxes, sigclip=sigclip) # make sure there are enough points to calculate a spectrum if len(stimes) > 9 and len(smags) > 9 and len(serrs) > 9: nbestsnrs = [] transitdepth, transitduration = [], [] ingressduration, points_in_transit, fit_status = [], [], [] nphasebins, transingressbin, transegressbin = [], [], [] # keep these around for diagnostics allsubtractedmags = [] allphasedmags = [] allphases = [] allblsmodels = [] # these are refit periods and epochs refitperiods = [] refitepochs = [] for period in nbestperiods: # get the period interval startp = period - perioddeltapercent*period/100.0 if startp < 0: startp = period endp = period + perioddeltapercent*period/100.0 # see if we need to rerun bls_serial_pfind if not assumeserialbls: # run bls_serial_pfind with the kwargs copied over from the # initial run. replace only the startp, endp, verbose, sigclip # kwarg values prevkwargs = blsdict['kwargs'].copy() prevkwargs['verbose'] = verbose prevkwargs['startp'] = startp prevkwargs['endp'] = endp prevkwargs['sigclip'] = None blsres = bls_serial_pfind(stimes, smags, serrs, **prevkwargs) else: blsres = blsdict thistransdepth = blsres['blsresult']['transdepth'] thistransduration = blsres['blsresult']['transduration'] thisbestperiod = blsres['bestperiod'] thistransingressbin = blsres['blsresult']['transingressbin'] thistransegressbin = blsres['blsresult']['transegressbin'] thisnphasebins = blsdict['kwargs']['nphasebins'] stats = _get_bls_stats(stimes, smags, serrs, thistransdepth, thistransduration, ingressdurationfraction, nphasebins, thistransingressbin, thistransegressbin, thisbestperiod, thisnphasebins, magsarefluxes=magsarefluxes, verbose=verbose) # update the lists with results from this peak nbestsnrs.append(stats['snr']) transitdepth.append(stats['transitdepth']) transitduration.append(stats['transitduration']) ingressduration.append(stats['ingressduration']) points_in_transit.append(stats['npoints_in_transit']) fit_status.append(stats['fit_status']) transingressbin.append(stats['transingressbin']) transegressbin.append(stats['transegressbin']) nphasebins.append(stats['nphasebins']) # update the refit periods and epochs refitperiods.append(stats['period']) refitepochs.append(stats['epoch']) # update the diagnostics allsubtractedmags.append(stats['subtractedmags']) allphasedmags.append(stats['phasedmags']) allphases.append(stats['phases']) allblsmodels.append(stats['blsmodel']) # # done with working on each peak # # if there aren't enough points in the mag series, bail out else: LOGERROR('no good detections for these times and mags, skipping...') nbestsnrs = None transitdepth, transitduration = None, None ingressduration, points_in_transit, fit_status = None, None, None nphasebins, transingressbin, transegressbin = None, None, None allsubtractedmags, allphases, allphasedmags = None, None, None return {'npeaks':npeaks, 'period':refitperiods, 'epoch':refitepochs, 'snr':nbestsnrs, 'transitdepth':transitdepth, 'transitduration':transitduration, 'ingressduration':ingressduration, 'npoints_in_transit':points_in_transit, 'fit_status':fit_status, 'nphasebins':nphasebins, 'transingressbin':transingressbin, 'transegressbin':transegressbin, 'allblsmodels':allblsmodels, 'allsubtractedmags':allsubtractedmags, 'allphasedmags':allphasedmags, 'allphases':allphases}