signature:
def absolute_gaia_magnitude ( gaia_mag , gaia_parallax_mas , gaia_mag_err = None , gaia_parallax_err_mas = None ) :
    '''Calculates the GAIA absolute magnitude for object ( or array of objects ) .
 Given a G mag and the parallax measured by GAIA , gets the absolute mag using
 the usual equation : :
 G - M _ G = 5 x log10 ( d _ pc ) - 5
 M _ G = 5 - 5log10 ( d _ pc ) + G
 Parameters
 gaia _ mag : float or array - like
 The measured GAIA G magnitude .
 gaia _ parallax _ max : float or array - like
 The measured parallax of the object in mas .
 gaia _ mag _ err : float or array - like or None
 The measurement error in GAIA G magnitude .
 gaia _ parallax _ err _ mas : float or array - like or None
 The measurement error in GAIA parallax in mas .
 Returns
 float or array - like
 The absolute magnitude M _ G of the object ( s ) .
 If both ` _ err ` input kwargs are provided , will return a tuple of the form : :
 ( M _ G float or array - like , M _ G _ err float or array - like )'''

implementation:
# get the distance
# we ' re using the naive calculation of d . this is inaccurate as stated in
# Bailer - Jones 2015 ( http : / / arxiv . org / abs / 1507.02105 ) if the error in
# parallax is a significant fraction of parallax
d_pc = np . abs ( 1000.0 / gaia_parallax_mas )
# get the distance error
if gaia_parallax_err_mas is not None :
    d_pc_err = ( ( 1000.0 / ( gaia_parallax_mas * gaia_parallax_mas ) ) * gaia_parallax_err_mas )
else :
    d_pc_err = None
# calculate the absolute mag from the relation
# FIXME : this is NOT corrected for extinction in G mag . see Jordi + 2010
# ( http : / / adsabs . harvard . edu / abs / 2010A % 26A . . . 523A . . 48J ) to figure out
# A _ G / A _ V as a function of ( V - I ) _ 0 , then apply it here
M_G = 5 - 5.0 * np . log10 ( d_pc ) + gaia_mag
# calculate the err in M _ G
if d_pc_err is not None and gaia_mag_err is not None :
    M_G_err = np . sqrt ( ( ( 5.0 / ( d_pc * np . log ( 10.0 ) ) ) ** 2 * ( d_pc_err ) ** 2 ) + gaia_mag_err * gaia_mag_err )
else :
    M_G_err = None
if M_G_err is not None :
    return M_G , M_G_err
else :
    return M_G


signature:
def aov_theta ( times , mags , errs , frequency , binsize = 0.05 , minbin = 9 ) :
    '''Calculates the Schwarzenberg - Czerny AoV statistic at a test frequency .
 Parameters
 times , mags , errs : np . array
 The input time - series and associated errors .
 frequency : float
 The test frequency to calculate the theta statistic at .
 binsize : float
 The phase bin size to use .
 minbin : int
 The minimum number of items in a phase bin to consider in the
 calculation of the statistic .
 Returns
 theta _ aov : float
 The value of the AoV statistic at the specified ` frequency ` .'''

implementation:
period = 1.0 / frequency
fold_time = times [ 0 ]
phased = phase_magseries ( times , mags , period , fold_time , wrap = False , sort = True )
phases = phased [ 'phase' ]
pmags = phased [ 'mags' ]
bins = nparange ( 0.0 , 1.0 , binsize )
ndets = phases . size
binnedphaseinds = npdigitize ( phases , bins )
bin_s1_tops = [ ]
bin_s2_tops = [ ]
binndets = [ ]
goodbins = 0
all_xbar = npmedian ( pmags )
for x in npunique ( binnedphaseinds ) :
    thisbin_inds = binnedphaseinds == x
    thisbin_mags = pmags [ thisbin_inds ]
    if thisbin_mags . size > minbin :
        thisbin_ndet = thisbin_mags . size
        thisbin_xbar = npmedian ( thisbin_mags )
        # get s1
        thisbin_s1_top = ( thisbin_ndet * ( thisbin_xbar - all_xbar ) * ( thisbin_xbar - all_xbar ) )
        # get s2
        thisbin_s2_top = npsum ( ( thisbin_mags - all_xbar ) * ( thisbin_mags - all_xbar ) )
        bin_s1_tops . append ( thisbin_s1_top )
        bin_s2_tops . append ( thisbin_s2_top )
        binndets . append ( thisbin_ndet )
        goodbins = goodbins + 1
# turn the quantities into arrays
bin_s1_tops = nparray ( bin_s1_tops )
bin_s2_tops = nparray ( bin_s2_tops )
binndets = nparray ( binndets )
# calculate s1 first
s1 = npsum ( bin_s1_tops ) / ( goodbins - 1.0 )
# then calculate s2
s2 = npsum ( bin_s2_tops ) / ( ndets - goodbins )
theta_aov = s1 / s2
return theta_aov


signature:
def _aov_worker ( task ) :
    '''This is a parallel worker for the function below .
 Parameters
 task : tuple
 This is of the form below : :
 task [ 0 ] = times
 task [ 1 ] = mags
 task [ 2 ] = errs
 task [ 3 ] = frequency
 task [ 4 ] = binsize
 task [ 5 ] = minbin
 Returns
 theta _ aov : float
 The theta value at the specified frequency . nan if the calculation
 fails .'''

implementation:
times , mags , errs , frequency , binsize , minbin = task
try :
    theta = aov_theta ( times , mags , errs , frequency , binsize = binsize , minbin = minbin )
    return theta
except Exception as e :
    return npnan


signature:
def get_frequency_grid ( times , samplesperpeak = 5 , nyquistfactor = 5 , minfreq = None , maxfreq = None , returnf0dfnf = False ) :
    '''This calculates a frequency grid for the period finding functions in this
 module .
 Based on the autofrequency function in astropy . stats . lombscargle .
 http : / / docs . astropy . org / en / stable / _ modules / astropy / stats / lombscargle / core . html # LombScargle . autofrequency
 Parameters
 times : np . array
 The times to use to generate the frequency grid over .
 samplesperpeak : int
 The minimum sample coverage each frequency point in the grid will get .
 nyquistfactor : int
 The multiplier over the Nyquist rate to use .
 minfreq , maxfreq : float or None
 If not None , these will be the limits of the frequency grid generated .
 returnf0dfnf : bool
 If this is True , will return the values of ` f0 ` , ` df ` , and ` Nf `
 generated for this grid .
 Returns
 np . array
 A grid of frequencies .'''

implementation:
baseline = times . max ( ) - times . min ( )
nsamples = times . size
df = 1. / baseline / samplesperpeak
if minfreq is not None :
    f0 = minfreq
else :
    f0 = 0.5 * df
if maxfreq is not None :
    Nf = int ( np . ceil ( ( maxfreq - f0 ) / df ) )
else :
    Nf = int ( 0.5 * samplesperpeak * nyquistfactor * nsamples )
if returnf0dfnf :
    return f0 , df , Nf , f0 + df * np . arange ( Nf )
else :
    return f0 + df * np . arange ( Nf )


signature:
def independent_freq_count ( frequencies , times , conservative = True ) :
    '''This estimates M : the number of independent frequencies in the periodogram .
 This follows the terminology on page 3 of Zechmeister & Kurster ( 2009 ) : :
 M = DELTA _ f / delta _ f
 where : :
 DELTA _ f = freq . max ( ) - freq . min ( )
 delta _ f = 1.0 / ( times . max ( ) - times . min ( ) )
 Parameters
 frequencies : np . array
 The frequencies array used for the calculation of the GLS periodogram .
 times : np . array
 The array of input times used for the calculation of the GLS
 periodogram .
 conservative : bool
 If True , will follow the prescription given in Schwarzenberg - Czerny
 (2003 ) :
 http : / / adsabs . harvard . edu / abs / 2003ASPC . . 292 . . 383S
 and estimate the number of independent frequences as : :
 min ( N _ obs , N _ freq , DELTA _ f / delta _ f )
 Returns
 M : int
 The number of independent frequencies .'''

implementation:
M = frequencies . ptp ( ) * times . ptp ( )
if conservative :
    M_eff = min ( [ times . size , frequencies . size , M ] )
else :
    M_eff = M
return M_eff


signature:
def bootstrap_falsealarmprob ( lspinfo , times , mags , errs , nbootstrap = 250 , magsarefluxes = False , sigclip = 10.0 , npeaks = None ) :
    '''Calculates the false alarm probabilities of periodogram peaks using
 bootstrap resampling of the magnitude time series .
 The false alarm probability here is defined as : :
 (1.0 + sum ( trialbestpeaks [ i ] > peak [ j ] ) ) / ( ntrialbestpeaks + 1)
 for each best periodogram peak j . The index i is for each bootstrap
 trial . This effectively gives us a significance for the peak . Smaller FAP
 means a better chance that the peak is real .
 The basic idea is to get the number of trial best peaks that are larger than
 the current best peak and divide this by the total number of trials . The
 distribution of these trial best peaks is obtained after scrambling the mag
 values and rerunning the specified periodogram method for a bunch of trials .
 ` lspinfo ` is the output dict from a periodbase periodogram function and MUST
 contain a ' method ' key that corresponds to one of the keys in the LSPMETHODS
 dict above . This will let this function know which periodogram function to
 run to generate the bootstrap samples . The lspinfo SHOULD also have a
 ' kwargs ' key that corresponds to the input keyword arguments for the
 periodogram function as it was run originally , to keep everything the same
 during the bootstrap runs . If this is missing , default values will be used .
 FIXME : this may not be strictly correct ; must look more into bootstrap
 significance testing . Also look into if we ' re doing resampling correctly for
 time series because the samples are not iid . Look into moving block
 bootstrap .
 Parameters
 lspinfo : dict
 A dict of period - finder results from one of the period - finders in
 periodbase , or your own functions , provided it ' s of the form and
 contains at least the keys listed below : :
 { ' periods ' : np . array of all periods searched by the period - finder ,
 ' lspvals ' : np . array of periodogram power value for each period ,
 ' bestperiod ' : a float value that is the period with the highest
 peak in the periodogram , i . e . the most - likely actual
 period ,
 ' method ' : a three - letter code naming the period - finder used ; must
 be one of the keys in the
 ` astrobase . periodbase . METHODLABELS ` dict ,
 ' nbestperiods ' : a list of the periods corresponding to periodogram
 peaks ( ` nbestlspvals ` below ) to annotate on the
 periodogram plot so they can be called out
 visually ,
 ' nbestlspvals ' : a list of the power values associated with
 periodogram peaks to annotate on the periodogram
 plot so they can be called out visually ; should be
 the same length as ` nbestperiods ` above ,
 ' kwargs ' : dict of kwargs passed to your own period - finder function }
 If you provide your own function ' s period - finder results , you should add
 a corresponding key for it to the LSPMETHODS dict above so the bootstrap
 function can use it correctly . Your period - finder function should take
 ` times ` , ` mags ` , errs and any extra parameters as kwargs and return a
 dict of the form described above . A small worked example : :
 from your _ module import your _ periodfinder _ func
 from astrobase import periodbase
 periodbase . LSPMETHODS [ ' your - finder ' ] = your _ periodfinder _ func
 # run a period - finder session
 your _ pfresults = your _ periodfinder _ func ( times , mags , errs ,
 * * extra _ kwargs )
 # run bootstrap to find FAP
 falsealarm _ info = periodbase . bootstrap _ falsealarmprob (
 your _ pfresults ,
 times , mags , errs ,
 nbootstrap = 250,
 magsarefluxes = False ,
 times , mags , errs : np . arrays
 The magnitude / flux time - series to process along with their associated
 measurement errors .
 nbootstrap : int
 The total number of bootstrap trials to run . This is set to 250 by
 default , but should probably be around 1000 for realistic results .
 magsarefluxes : bool
 If True , indicates the input time - series is fluxes and not mags .
 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 .
 npeaks : int or None
 The number of peaks from the list of ' nbestlspvals ' in the period - finder
 result dict to run the bootstrap for . If None , all of the peaks in this
 list will have their FAP calculated .
 Returns
 dict
 Returns a dict of the form : :
 { ' peaks ' : allpeaks ,
 ' periods ' : allperiods ,
 ' probabilities ' : allfaps ,
 ' alltrialbestpeaks ' : alltrialbestpeaks }'''

implementation:
# figure out how many periods to work on
if ( npeaks and ( 0 < npeaks < len ( lspinfo [ 'nbestperiods' ] ) ) ) :
    nperiods = npeaks
else :
    LOGWARNING ( 'npeaks not specified or invalid, ' 'getting FAP for all %s periodogram peaks' % len ( lspinfo [ 'nbestperiods' ] ) )
    nperiods = len ( lspinfo [ 'nbestperiods' ] )
nbestperiods = lspinfo [ 'nbestperiods' ] [ : nperiods ]
nbestpeaks = lspinfo [ 'nbestlspvals' ] [ : nperiods ]
# get rid of nans first and sigclip
stimes , smags , serrs = sigclip_magseries ( times , mags , errs , magsarefluxes = magsarefluxes , sigclip = sigclip )
allpeaks = [ ]
allperiods = [ ]
allfaps = [ ]
alltrialbestpeaks = [ ]
# make sure there are enough points to calculate a spectrum
if len ( stimes ) > 9 and len ( smags ) > 9 and len ( serrs ) > 9 :
    for ind , period , peak in zip ( range ( len ( nbestperiods ) ) , nbestperiods , nbestpeaks ) :
        LOGINFO ( 'peak %s: running %s trials...' % ( ind + 1 , nbootstrap ) )
        trialbestpeaks = [ ]
        for _trial in range ( nbootstrap ) : # get a scrambled index
            tindex = np . random . randint ( 0 , high = mags . size , size = mags . size )
            # get the kwargs dict out of the lspinfo
            if 'kwargs' in lspinfo :
                kwargs = lspinfo [ 'kwargs' ]
                # update the kwargs with some local stuff
                kwargs . update ( { 'magsarefluxes' : magsarefluxes , 'sigclip' : sigclip , 'verbose' : False } )
            else :
                kwargs = { 'magsarefluxes' : magsarefluxes , 'sigclip' : sigclip , 'verbose' : False }
            # run the periodogram with scrambled mags and errs
            # and the appropriate keyword arguments
            lspres = LSPMETHODS [ lspinfo [ 'method' ] ] ( times , mags [ tindex ] , errs [ tindex ] , ** kwargs )
            trialbestpeaks . append ( lspres [ 'bestlspval' ] )
        trialbestpeaks = np . array ( trialbestpeaks )
        alltrialbestpeaks . append ( trialbestpeaks )
        # calculate the FAP for a trial peak j = FAP [ j ] =
        # (1.0 + sum ( trialbestpeaks [ i ] > peak [ j ] ) ) / ( ntrialbestpeaks + 1)
        if lspinfo [ 'method' ] != 'pdm' :
            falsealarmprob = ( ( 1.0 + trialbestpeaks [ trialbestpeaks > peak ] . size ) / ( trialbestpeaks . size + 1.0 ) )
        # for PDM , we ' re looking for a peak smaller than the best peak
        # because values closer to 0.0 are more significant
        else :
            falsealarmprob = ( ( 1.0 + trialbestpeaks [ trialbestpeaks < peak ] . size ) / ( trialbestpeaks . size + 1.0 ) )
        LOGINFO ( 'FAP for peak %s, period: %.6f = %.3g' % ( ind + 1 , period , falsealarmprob ) )
        allpeaks . append ( peak )
        allperiods . append ( period )
        allfaps . append ( falsealarmprob )
    return { 'peaks' : allpeaks , 'periods' : allperiods , 'probabilities' : allfaps , 'alltrialbestpeaks' : alltrialbestpeaks }
else :
    LOGERROR ( 'not enough mag series points to calculate periodogram' )
    return None


signature:
def make_combined_periodogram ( pflist , outfile , addmethods = False ) :
    '''This just puts all of the period - finders on a single periodogram .
 This will renormalize all of the periodograms so their values lie between 0
 and 1 , with values lying closer to 1 being more significant . Periodograms
 that give the same best periods will have their peaks line up together .
 Parameters
 pflist : list of dict
 This is a list of result dicts from any of the period - finders in
 periodbase . To use your own period - finders ' results here , make sure the
 result dict is of the form and has at least the keys below : :
 { ' periods ' : np . array of all periods searched by the period - finder ,
 ' lspvals ' : np . array of periodogram power value for each period ,
 ' bestperiod ' : a float value that is the period with the highest
 peak in the periodogram , i . e . the most - likely actual
 period ,
 ' method ' : a three - letter code naming the period - finder used ; must
 be one of the keys in the
 ` astrobase . periodbase . METHODLABELS ` dict ,
 ' nbestperiods ' : a list of the periods corresponding to periodogram
 peaks ( ` nbestlspvals ` below ) to annotate on the
 periodogram plot so they can be called out
 visually ,
 ' nbestlspvals ' : a list of the power values associated with
 periodogram peaks to annotate on the periodogram
 plot so they can be called out visually ; should be
 the same length as ` nbestperiods ` above ,
 ' kwargs ' : dict of kwargs passed to your own period - finder function }
 outfile : str
 This is the output file to write the output to . NOTE : EPS / PS won ' t work
 because we use alpha transparency to better distinguish between the
 various periodograms .
 addmethods : bool
 If this is True , will add all of the normalized periodograms together ,
 then renormalize them to between 0 and 1 . In this way , if all of the
 period - finders agree on something , it ' ll stand out easily . FIXME :
 implement this kwarg .
 Returns
 str
 The name of the generated plot file .'''

implementation:
import matplotlib . pyplot as plt
for pf in pflist :
    if pf [ 'method' ] == 'pdm' :
        plt . plot ( pf [ 'periods' ] , np . max ( pf [ 'lspvals' ] ) / pf [ 'lspvals' ] - 1.0 , label = '%s P=%.5f' % ( pf [ 'method' ] , pf [ 'bestperiod' ] ) , alpha = 0.5 )
    else :
        plt . plot ( pf [ 'periods' ] , pf [ 'lspvals' ] / np . max ( pf [ 'lspvals' ] ) , label = '%s P=%.5f' % ( pf [ 'method' ] , pf [ 'bestperiod' ] ) , alpha = 0.5 )
plt . xlabel ( 'period [days]' )
plt . ylabel ( 'normalized periodogram power' )
plt . xscale ( 'log' )
plt . legend ( )
plt . tight_layout ( )
plt . savefig ( outfile )
plt . close ( 'all' )
return outfile


signature:
def read_original_textlc ( lcpath ) :
    '''Read . epdlc , and . tfalc light curves and return a corresponding labelled
 dict ( if LC from < 2012 ) or astropy table ( if > = 2012 ) . Each has different
 keys that can be accessed via . keys ( )
 Input :
 lcpath : path ( string ) to light curve data , which is a textfile with HAT
 LC data .
 Example :
 dat = read _ original _ textlc ( ' HAT - 115-0003266 . epdlc ' )'''

implementation:
LOGINFO ( 'reading original HAT text LC: {:s}' . format ( lcpath ) )
N_lines_to_parse_comments = 50
with open ( lcpath , 'rb' ) as file :
    head = [ next ( file ) for ind in range ( N_lines_to_parse_comments ) ]
N_comment_lines = len ( [ l for l in head if l . decode ( 'UTF-8' ) [ 0 ] == '#' ] )
# if there are too many comment lines , fail out
if N_comment_lines < N_lines_to_parse_comments :
    LOGERROR ( 'LC file {fpath} has too many comment lines' . format ( fpath = lcpath ) )
    return None
first_data_line = list ( filter ( None , head [ N_comment_lines ] . decode ( 'UTF-8' ) . split ( ) ) )
N_cols = len ( first_data_line )
# There are different column formats depending on when HAT pipeline was run
# also different formats for different types of LCs :
# pre - 2012 : . epdlc - > 17 columns
# pre - 2012 : . tfalc - > 20 columns
# post - 2012 : . epdlc or . tfalc - > 32 columns
if N_cols == 17 :
    colformat = 'pre2012-epdlc'
elif N_cols == 20 :
    colformat = 'pre2012-tfalc'
elif N_cols == 32 :
    colformat = 'post2012-hatlc'
else :
    LOGERROR ( "can't handle this column format yet, " "file: {fpath}, ncols: {ncols}" . format ( fpath = lcpath , ncols = N_cols ) )
    return None
# deal with pre - 2012 column format
if colformat == 'pre2012-epdlc' :
    col_names = [ 'framekey' , 'rjd' , 'aim_000' , 'aie_000' , 'aiq_000' , 'aim_001' , 'aie_001' , 'aiq_001' , 'aim_002' , 'aie_002' , 'aiq_002' , 'arm_000' , 'arm_001' , 'arm_002' , 'aep_000' , 'aep_001' , 'aep_002' ]
    col_dtypes = [ 'U8' , float , float , float , 'U1' , float , float , 'U1' , float , float , 'U1' , float , float , float , float , float , float ]
    dtype_pairs = [ el for el in zip ( col_names , col_dtypes ) ]
    data = np . genfromtxt ( lcpath , names = col_names , dtype = col_dtypes , skip_header = N_comment_lines , delimiter = None )
    out = { }
    for ix in range ( len ( data . dtype . names ) ) :
        out [ data . dtype . names [ ix ] ] = data [ data . dtype . names [ ix ] ]
elif colformat == 'pre2012-tfalc' :
    col_names = [ 'framekey' , 'rjd' , 'aim_000' , 'aie_000' , 'aiq_000' , 'aim_001' , 'aie_001' , 'aiq_001' , 'aim_002' , 'aie_002' , 'aiq_002' , 'arm_000' , 'arm_001' , 'arm_002' , 'aep_000' , 'aep_001' , 'aep_002' , 'atf_000' , 'atf_001' , 'atf_002' ]
    col_dtypes = [ 'U8' , float , float , float , 'U1' , float , float , 'U1' , float , float , 'U1' , float , float , float , float , float , float , float , float , float ]
    dtype_pairs = [ el for el in zip ( col_names , col_dtypes ) ]
    data = np . genfromtxt ( lcpath , names = col_names , dtype = col_dtypes , skip_header = N_comment_lines , delimiter = None )
    out = { }
    for ix in range ( len ( data . dtype . names ) ) :
        out [ data . dtype . names [ ix ] ] = data [ data . dtype . names [ ix ] ]
elif colformat == 'post2012-hatlc' :
    col_names = [ 'hatid' , 'framekey' , 'fld' , 'bjd' , 'aim_000' , 'aie_000' , 'aiq_000' , 'aim_001' , 'aie_001' , 'aiq_001' , 'aim_002' , 'aie_002' , 'aiq_002' , 'arm_000' , 'arm_001' , 'arm_002' , 'aep_000' , 'aep_001' , 'aep_002' , 'atf_000' , 'atf_001' , 'atf_002' , 'xcc' , 'ycc' , 'bgv' , 'bge' , 'fsv' , 'fdv' , 'fkv' , 'iha' , 'izd' , 'rjd' ]
    out = astascii . read ( lcpath , names = col_names , comment = '#' )
return out


signature:
def traptransit_fit_magseries ( times , mags , errs , transitparams , sigclip = 10.0 , plotfit = False , magsarefluxes = False , verbose = True ) :
    '''This fits a trapezoid transit model to a magnitude time series .
 Parameters
 times , mags , errs : np . array
 The input mag / flux time - series to fit a trapezoid planet - transit model
 to .
 period : float
 The period to use for the model fit .
 transitparams : list of floats
 These are initial parameters for the transit model fit . A list of the
 following form is required : :
 transitparams = [ transitperiod ( time ) ,
 transitepoch ( time ) ,
 transitdepth ( flux or mags ) ,
 transitduration ( phase ) ,
 ingressduration ( phase ) ]
 - for magnitudes - > ` transitdepth ` should be < 0
 - for fluxes - > ` transitdepth ` should be > 0
 If ` transitepoch ` is None , this function will do an initial spline fit
 to find an approximate minimum of the phased light curve using the given
 period .
 The ` transitdepth ` provided is checked against the value of
 ` magsarefluxes ` . if ` magsarefluxes = True ` , the ` transitdepth ` is forced
 to be > 0 ; if ` magsarefluxes ` = False , the ` transitdepth ` is forced to
 be < 0.
 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 .
 magsarefluxes : bool
 If True , will treat the input values of ` mags ` as fluxes for purposes of
 plotting the fit and sig - clipping .
 plotfit : str or False
 If this is a string , this function will make a plot for the fit to the
 mag / flux time - series and writes the plot to the path specified here .
 ignoreinitfail : bool
 If this is True , ignores the initial failure to find a set of optimized
 Fourier parameters using the global optimization function and proceeds
 to do a least - squares fit anyway .
 verbose : bool
 If True , will indicate progress and warn of any problems .
 Returns
 dict
 This function returns a dict containing the model fit parameters , the
 minimized chi - sq value and the reduced chi - sq value . The form of this
 dict is mostly standardized across all functions in this module : :
 ' fittype ' : ' traptransit ' ,
 ' fitinfo ' : {
 ' initialparams ' : the initial transit params provided ,
 ' finalparams ' : the final model fit transit params ,
 ' finalparamerrs ' : formal errors in the params ,
 ' leastsqfit ' : the full tuple returned by scipy . leastsq ,
 ' fitmags ' : the model fit mags ,
 ' fitepoch ' : the epoch of minimum light for the fit ,
 ' ntransitpoints ' : the number of LC points in transit phase
 ' fitchisq ' : the minimized value of the fit ' s chi - sq ,
 ' fitredchisq ' : the reduced chi - sq value ,
 ' fitplotfile ' : the output fit plot if fitplot is not None ,
 ' magseries ' : {
 ' times ' : input times in phase order of the model ,
 ' phase ' : the phases of the model mags ,
 ' mags ' : input mags / fluxes in the phase order of the model ,
 ' errs ' : errs in the phase order of the model ,
 ' magsarefluxes ' : input value of magsarefluxes kwarg'''

implementation:
stimes , smags , serrs = sigclip_magseries ( times , mags , errs , sigclip = sigclip , magsarefluxes = magsarefluxes )
# get rid of zero errs
nzind = np . nonzero ( serrs )
stimes , smags , serrs = stimes [ nzind ] , smags [ nzind ] , serrs [ nzind ]
# check the transitparams
transitperiod , transitepoch , transitdepth = transitparams [ 0 : 3 ]
# check if we have a transitepoch to use
if transitepoch is None :
    if verbose :
        LOGWARNING ( 'no transitepoch given in transitparams, ' 'trying to figure it out automatically...' )
    # do a spline fit to figure out the approximate min of the LC
    try :
        spfit = spline_fit_magseries ( times , mags , errs , transitperiod , sigclip = sigclip , magsarefluxes = magsarefluxes , verbose = verbose )
        transitepoch = spfit [ 'fitinfo' ] [ 'fitepoch' ]
    # if the spline - fit fails , try a savgol fit instead
    except Exception as e :
        sgfit = savgol_fit_magseries ( times , mags , errs , transitperiod , sigclip = sigclip , magsarefluxes = magsarefluxes , verbose = verbose )
        transitepoch = sgfit [ 'fitinfo' ] [ 'fitepoch' ]
    # if everything failed , then bail out and ask for the transitepoch
    finally :
        if transitepoch is None :
            LOGERROR ( "couldn't automatically figure out the transit epoch, " "can't continue. please provide it in transitparams." )
            # assemble the returndict
            returndict = { 'fittype' : 'traptransit' , 'fitinfo' : { 'initialparams' : transitparams , 'finalparams' : None , 'leastsqfit' : None , 'fitmags' : None , 'fitepoch' : None , } , 'fitchisq' : np . nan , 'fitredchisq' : np . nan , 'fitplotfile' : None , 'magseries' : { 'phase' : None , 'times' : None , 'mags' : None , 'errs' : None , 'magsarefluxes' : magsarefluxes , } , }
            return returndict
        else : # check the case when there are more than one transitepochs
        # returned
            if transitepoch . size > 1 :
                if verbose :
                    LOGWARNING ( "could not auto-find a single minimum in LC for " "transitepoch, using the first one returned" )
                transitparams [ 1 ] = transitepoch [ 0 ]
            else :
                if verbose :
                    LOGWARNING ( 'using automatically determined transitepoch = %.5f' % transitepoch )
                transitparams [ 1 ] = transitepoch . item ( )
# next , check the transitdepth and fix it to the form required
if magsarefluxes :
    if transitdepth < 0.0 :
        transitparams [ 2 ] = - transitdepth
else :
    if transitdepth > 0.0 :
        transitparams [ 2 ] = - transitdepth
# finally , do the fit
try :
    leastsqfit = spleastsq ( transits . trapezoid_transit_residual , transitparams , args = ( stimes , smags , serrs ) , full_output = True )
except Exception as e :
    leastsqfit = None
# if the fit succeeded , then we can return the final parameters
if leastsqfit and leastsqfit [ - 1 ] in ( 1 , 2 , 3 , 4 ) :
    finalparams = leastsqfit [ 0 ]
    covxmatrix = leastsqfit [ 1 ]
    # calculate the chisq and reduced chisq
    fitmags , phase , ptimes , pmags , perrs , n_transitpoints = ( transits . trapezoid_transit_func ( finalparams , stimes , smags , serrs , get_ntransitpoints = True ) )
    fitchisq = np . sum ( ( ( fitmags - pmags ) * ( fitmags - pmags ) ) / ( perrs * perrs ) )
    fitredchisq = fitchisq / ( len ( pmags ) - len ( finalparams ) - 1 )
    # get the residual variance and calculate the formal 1 - sigma errs on the
    # final parameters
    residuals = leastsqfit [ 2 ] [ 'fvec' ]
    residualvariance = ( np . sum ( residuals * residuals ) / ( pmags . size - finalparams . size ) )
    if covxmatrix is not None :
        covmatrix = residualvariance * covxmatrix
        stderrs = np . sqrt ( np . diag ( covmatrix ) )
    else :
        LOGERROR ( 'covxmatrix not available, fit probably failed!' )
        stderrs = None
    if verbose :
        LOGINFO ( 'final fit done. chisq = %.5f, reduced chisq = %.5f' % ( fitchisq , fitredchisq ) )
    # get the fit epoch
    fperiod , fepoch = finalparams [ : 2 ]
    # assemble the returndict
    returndict = { 'fittype' : 'traptransit' , 'fitinfo' : { 'initialparams' : transitparams , 'finalparams' : finalparams , 'finalparamerrs' : stderrs , 'leastsqfit' : leastsqfit , 'fitmags' : fitmags , 'fitepoch' : fepoch , 'ntransitpoints' : n_transitpoints } , 'fitchisq' : fitchisq , 'fitredchisq' : fitredchisq , 'fitplotfile' : None , 'magseries' : { 'phase' : phase , 'times' : ptimes , 'mags' : pmags , 'errs' : perrs , 'magsarefluxes' : magsarefluxes , } , }
    # make the fit plot if required
    if plotfit and isinstance ( plotfit , str ) :
        make_fit_plot ( phase , pmags , perrs , fitmags , fperiod , ptimes . min ( ) , fepoch , plotfit , magsarefluxes = magsarefluxes )
        returndict [ 'fitplotfile' ] = plotfit
    return returndict
# if the leastsq fit failed , return nothing
else :
    LOGERROR ( 'trapezoid-fit: least-squared fit to the light curve failed!' )
    # assemble the returndict
    returndict = { 'fittype' : 'traptransit' , 'fitinfo' : { 'initialparams' : transitparams , 'finalparams' : None , 'finalparamerrs' : None , 'leastsqfit' : leastsqfit , 'fitmags' : None , 'fitepoch' : None , 'ntransitpoints' : 0 } , 'fitchisq' : np . nan , 'fitredchisq' : np . nan , 'fitplotfile' : None , 'magseries' : { 'phase' : None , 'times' : None , 'mags' : None , 'errs' : None , 'magsarefluxes' : magsarefluxes , } , }
    return returndict


signature:
def _get_value ( quantitystr , fitparams , fixedparams ) :
    """This decides if a value is to be fit for or is fixed in a model fit .
 When you want to get the value of some parameter , but you ' re not sure if
 it ' s being fit or if it is fixed . then , e . g . for ` period ` : :
 period _ value = _ get _ value ( ' period ' , fitparams , fixedparams )"""

implementation:
# for Mandel - Agol fitting , sometimes we want to fix some parameters ,
# and fit others . this function allows that flexibility .
fitparamskeys , fixedparamskeys = fitparams . keys ( ) , fixedparams . keys ( )
if quantitystr in fitparamskeys :
    quantity = fitparams [ quantitystr ]
elif quantitystr in fixedparamskeys :
    quantity = fixedparams [ quantitystr ]
return quantity


signature:
def _transit_model ( times , t0 , per , rp , a , inc , ecc , w , u , limb_dark , exp_time_minutes = 2 , supersample_factor = 7 ) :
    '''This returns a BATMAN planetary transit model .
 Parameters
 times : np . array
 The times at which the model will be evaluated .
 t0 : float
 The time of periastron for the transit .
 per : float
 The orbital period of the planet .
 rp : float
 The stellar radius of the planet ' s star ( in Rsun ) .
 a : float
 The semi - major axis of the planet ' s orbit ( in Rsun ) .
 inc : float
 The orbital inclination ( in degrees ) .
 ecc : float
 The eccentricity of the orbit .
 w : float
 The longitude of periastron ( in degrees ) .
 u : list of floats
 The limb darkening coefficients specific to the limb darkening model
 used .
 limb _ dark : { " uniform " , " linear " , " quadratic " , " square - root " , " logarithmic " , " exponential " , " power2 " , " custom " }
 The type of limb darkening model to use . See the full list here :
 https : / / www . cfa . harvard . edu / ~ lkreidberg / batman / tutorial . html # limb - darkening - options
 exp _ time _ minutes : float
 The amount of time to ' smear ' the transit LC points over to simulate a
 long exposure time .
 supersample _ factor : int
 The number of supersampled time data points to average the lightcurve
 model over .
 Returns
 ( params , batman _ model ) : tuple
 The returned tuple contains the params list and the generated
 ` batman . TransitModel ` object .'''

implementation:
params = batman . TransitParams ( )
# object to store transit parameters
params . t0 = t0
# time of periastron
params . per = per
# orbital period
params . rp = rp
# planet radius ( in stellar radii )
params . a = a
# semi - major axis ( in stellar radii )
params . inc = inc
# orbital inclination ( in degrees )
params . ecc = ecc
# the eccentricity of the orbit
params . w = w
# longitude of periastron ( in degrees )
params . u = u
# limb darkening coefficient list
params . limb_dark = limb_dark
# limb darkening model to use
t = times
m = batman . TransitModel ( params , t , exp_time = exp_time_minutes / 60. / 24. , supersample_factor = supersample_factor )
return params , m


signature:
def _log_prior_transit ( theta , priorbounds ) :
    '''Assume priors on all parameters have uniform probability .'''

implementation:
# priorbounds contains the input priors , and because of how we previously
# sorted theta , its sorted keys tell us which parts of theta correspond to
# which physical quantities .
allowed = True
for ix , key in enumerate ( np . sort ( list ( priorbounds . keys ( ) ) ) ) :
    if priorbounds [ key ] [ 0 ] < theta [ ix ] < priorbounds [ key ] [ 1 ] :
        allowed = True and allowed
    else :
        allowed = False
if allowed :
    return 0.
return - np . inf


signature:
def _log_likelihood_transit ( theta , params , model , t , flux , err_flux , priorbounds ) :
    '''Given a batman TransitModel and its proposed parameters ( theta ) , update the
 batman params object with the proposed parameters and evaluate the gaussian
 likelihood .
 Note : the priorbounds are only needed to parse theta .'''

implementation:
u = [ ]
for ix , key in enumerate ( sorted ( priorbounds . keys ( ) ) ) :
    if key == 'rp' :
        params . rp = theta [ ix ]
    elif key == 't0' :
        params . t0 = theta [ ix ]
    elif key == 'sma' :
        params . a = theta [ ix ]
    elif key == 'incl' :
        params . inc = theta [ ix ]
    elif key == 'period' :
        params . per = theta [ ix ]
    elif key == 'ecc' :
        params . per = theta [ ix ]
    elif key == 'omega' :
        params . w = theta [ ix ]
    elif key == 'u_linear' :
        u . append ( theta [ ix ] )
    elif key == 'u_quadratic' :
        u . append ( theta [ ix ] )
        params . u = u
lc = model . light_curve ( params )
residuals = flux - lc
log_likelihood = - 0.5 * ( np . sum ( ( residuals / err_flux ) ** 2 + np . log ( 2 * np . pi * ( err_flux ) ** 2 ) ) )
return log_likelihood


signature:
def _log_likelihood_transit_plus_line ( theta , params , model , t , data_flux , err_flux , priorbounds ) :
    '''Given a batman TransitModel and its proposed parameters ( theta ) , update the
 batman params object with the proposed parameters and evaluate the gaussian
 likelihood .
 Note : the priorbounds are only needed to parse theta .'''

implementation:
u = [ ]
for ix , key in enumerate ( sorted ( priorbounds . keys ( ) ) ) :
    if key == 'rp' :
        params . rp = theta [ ix ]
    elif key == 't0' :
        params . t0 = theta [ ix ]
    elif key == 'sma' :
        params . a = theta [ ix ]
    elif key == 'incl' :
        params . inc = theta [ ix ]
    elif key == 'period' :
        params . per = theta [ ix ]
    elif key == 'ecc' :
        params . per = theta [ ix ]
    elif key == 'omega' :
        params . w = theta [ ix ]
    elif key == 'u_linear' :
        u . append ( theta [ ix ] )
    elif key == 'u_quadratic' :
        u . append ( theta [ ix ] )
        params . u = u
    elif key == 'poly_order0' :
        poly_order0 = theta [ ix ]
    elif key == 'poly_order1' :
        poly_order1 = theta [ ix ]
try :
    poly_order0
except Exception as e :
    poly_order0 = 0
else :
    pass
transit = model . light_curve ( params )
line = poly_order0 + t * poly_order1
model = transit + line
residuals = data_flux - model
log_likelihood = - 0.5 * ( np . sum ( ( residuals / err_flux ) ** 2 + np . log ( 2 * np . pi * ( err_flux ) ** 2 ) ) )
return log_likelihood


signature:
def log_posterior_transit ( theta , params , model , t , flux , err_flux , priorbounds ) :
    '''Evaluate posterior probability given proposed model parameters and
 the observed flux timeseries .'''

implementation:
lp = _log_prior_transit ( theta , priorbounds )
if not np . isfinite ( lp ) :
    return - np . inf
else :
    return lp + _log_likelihood_transit ( theta , params , model , t , flux , err_flux , priorbounds )


signature:
def log_posterior_transit_plus_line ( theta , params , model , t , flux , err_flux , priorbounds ) :
    '''Evaluate posterior probability given proposed model parameters and
 the observed flux timeseries .'''

implementation:
lp = _log_prior_transit_plus_line ( theta , priorbounds )
if not np . isfinite ( lp ) :
    return - np . inf
else :
    return ( lp + _log_likelihood_transit_plus_line ( theta , params , model , t , flux , err_flux , priorbounds ) )


signature:
def mandelagol_fit_magseries ( times , mags , errs , fitparams , priorbounds , fixedparams , trueparams = None , burninpercent = 0.3 , plotcorner = False , samplesavpath = False , n_walkers = 50 , n_mcmc_steps = 400 , eps = 1e-4 , skipsampling = False , overwriteexistingsamples = False , mcmcprogressbar = False , plotfit = False , magsarefluxes = False , sigclip = 10.0 , verbose = True , nworkers = 4 ) :
    '''This fits a Mandel & Agol ( 2002 ) planetary transit model to a flux time
 series . You can fit and fix whatever parameters you want .
 It relies on Kreidberg ( 2015 ) ' s BATMAN implementation for the transit model ,
 emcee to sample the posterior ( Foreman - Mackey et al 2013 ) , ` corner ` to plot
 it , and ` h5py ` to save the samples . See e . g . , Claret ' s work for good guesses
 of star - appropriate limb - darkening parameters .
 NOTE : this only works for flux time - series at the moment .
 NOTE : Between the ` fitparams ` , ` priorbounds ` , and ` fixedparams ` dicts , you
 must specify all of the planetary transit parameters required by BATMAN :
 ` [ ' t0 ' , ' rp ' , ' sma ' , ' incl ' , ' u ' , ' rp ' , ' ecc ' , ' omega ' , ' period ' ] ` , or the
 BATMAN model will fail to initialize .
 Parameters
 times , mags , errs : np . array
 The input flux time - series to fit a Fourier cosine series to .
 fitparams : dict
 This is the initial parameter guesses for MCMC , found e . g . , by
 BLS . The key string format must not be changed , but any parameter can be
 either " fit " or " fixed " . If it is " fit " , it must have a corresponding
 prior . For example : :
 fitparams = { ' t0 ' : 1325.9 , ' rp ' : np . sqrt ( fitd [ ' transitdepth ' ] ) ,
 ' sma ' : 6.17 , ' incl ' : 85 , ' u ' : [ 0.3 , 0.2 ] }
 where ' u ' is a list of the limb darkening parameters , Linear first , then
 quadratic . Quadratic limb darkening is the only form implemented .
 priorbounds : dict
 This sets the lower & upper bounds on uniform prior , e . g . : :
 priorbounds = { ' rp ' : ( 0.135 , 0.145 ) , ' u _ linear ' : ( 0.3-1 , 0.3 + 1 ) ,
 ' u _ quad ' : ( 0.2-1 , 0.2 + 1 ) , ' t0 ' : ( np . min ( time ) ,
 np . max ( time ) ) , ' sma ' : ( 6,6.4 ) , ' incl ' : ( 80,90 ) }
 fixedparams : dict
 This sets which parameters are fixed , and their values . For example : :
 fixedparams = { ' ecc ' : 0 . ,
 ' omega ' : 90 . ,
 ' limb _ dark ' : ' quadratic ' ,
 ' period ' : fitd [ ' period ' ] }
 ` limb _ dark ` must be " quadratic " . It ' s " fixed " , because once you
 choose your limb - darkening model , it ' s fixed .
 trueparams : list of floats
 The true parameter values you ' re fitting for , if they ' re known ( e . g . , a
 known planet , or fake data ) . Only for plotting purposes .
 burninpercent : float
 The percent of MCMC samples to discard as burn - in .
 plotcorner : str or False
 If this is a str , points to the path of output corner plot that will be
 generated for this MCMC run .
 samplesavpath : str
 This must be provided so ` emcee ` can save its MCMC samples to disk as
 HDF5 files . This will set the path of the output HDF5file written .
 n _ walkers : int
 The number of MCMC walkers to use .
 n _ mcmc _ steps : int
 The number of MCMC steps to take .
 eps : float
 The radius of the ` n _ walkers - dimensional ` Gaussian ball used to
 initialize the MCMC .
 skipsampling : bool
 If you ' ve already collected MCMC samples , and you do not want any more
 sampling ( e . g . , just make the plots ) , set this to be True .
 overwriteexistingsamples : bool
 If you ' ve collected samples , but you want to overwrite them , set this to
 True . Usually , it should be False , which appends samples to
 ` samplesavpath ` HDF5 file .
 mcmcprogressbar : bool
 If True , will show a progress bar for the MCMC process .
 plotfit : str or bool
 If a str , indicates the path of the output fit plot file . If False , no
 fit plot will be made .
 magsarefluxes : bool
 This indicates if the input measurements in ` mags ` are actually fluxes .
 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 .
 verbose : bool
 If True , will indicate MCMC progress .
 nworkers : int
 The number of parallel workers to launch for MCMC .
 Returns
 dict
 This function returns a dict containing the model fit parameters and
 other fit information . The form of this dict is mostly standardized
 across all functions in this module : :
 ' fittype ' : ' mandelagol ' ,
 ' fitinfo ' : {
 ' initialparams ' : the initial transit params provided ,
 ' fixedparams ' : the fixed transit params provided ,
 ' finalparams ' : the final model fit transit params ,
 ' finalparamerrs ' : formal errors in the params ,
 ' fitmags ' : the model fit mags ,
 ' fitepoch ' : the epoch of minimum light for the fit ,
 ' fitplotfile ' : the output fit plot if fitplot is not None ,
 ' magseries ' : {
 ' times ' : input times in phase order of the model ,
 ' phase ' : the phases of the model mags ,
 ' mags ' : input mags / fluxes in the phase order of the model ,
 ' errs ' : errs in the phase order of the model ,
 ' magsarefluxes ' : input value of magsarefluxes kwarg'''

implementation:
from multiprocessing import Pool
fittype = 'mandelagol'
if not magsarefluxes :
    raise NotImplementedError ( 'magsarefluxes is not implemented yet.' )
if not samplesavpath :
    raise ValueError ( 'This function requires that you save the samples somewhere' )
if not mandel_agol_dependencies :
    raise ImportError ( 'This function depends on BATMAN, emcee>3.0, corner, and h5py.' )
# sigma clip and get rid of zero errs
stimes , smags , serrs = sigclip_magseries ( times , mags , errs , sigclip = sigclip , magsarefluxes = magsarefluxes )
nzind = np . nonzero ( serrs )
stimes , smags , serrs = stimes [ nzind ] , smags [ nzind ] , serrs [ nzind ]
init_period = _get_value ( 'period' , fitparams , fixedparams )
init_epoch = _get_value ( 't0' , fitparams , fixedparams )
init_rp = _get_value ( 'rp' , fitparams , fixedparams )
init_sma = _get_value ( 'sma' , fitparams , fixedparams )
init_incl = _get_value ( 'incl' , fitparams , fixedparams )
init_ecc = _get_value ( 'ecc' , fitparams , fixedparams )
init_omega = _get_value ( 'omega' , fitparams , fixedparams )
limb_dark = _get_value ( 'limb_dark' , fitparams , fixedparams )
init_u = _get_value ( 'u' , fitparams , fixedparams )
if not limb_dark == 'quadratic' :
    raise ValueError ( 'only quadratic limb-darkening is supported at the moment' )
# initialize the model and calculate the initial model light - curve
init_params , init_m = _transit_model ( stimes , init_epoch , init_period , init_rp , init_sma , init_incl , init_ecc , init_omega , init_u , limb_dark )
init_flux = init_m . light_curve ( init_params )
# guessed initial params . give nice guesses , or else emcee struggles .
theta , fitparamnames = [ ] , [ ]
for k in np . sort ( list ( fitparams . keys ( ) ) ) :
    if isinstance ( fitparams [ k ] , float ) or isinstance ( fitparams [ k ] , int ) :
        theta . append ( fitparams [ k ] )
        fitparamnames . append ( fitparams [ k ] )
    elif isinstance ( fitparams [ k ] , list ) :
        if not len ( fitparams [ k ] ) == 2 :
            raise ValueError ( 'should only be quadratic LD coeffs' )
        theta . append ( fitparams [ k ] [ 0 ] )
        theta . append ( fitparams [ k ] [ 1 ] )
        fitparamnames . append ( fitparams [ k ] [ 0 ] )
        fitparamnames . append ( fitparams [ k ] [ 1 ] )
# initialize sampler
n_dim = len ( theta )
initial_position_vec = [ theta + eps * np . random . randn ( n_dim ) for i in range ( n_walkers ) ]
# run the MCMC , unless you just want to load the available samples
if not skipsampling :
    backend = emcee . backends . HDFBackend ( samplesavpath )
    if overwriteexistingsamples :
        LOGWARNING ( 'erased samples previously at {:s}' . format ( samplesavpath ) )
        backend . reset ( n_walkers , n_dim )
    # if this is the first run , then start from a gaussian ball .
    # otherwise , resume from the previous samples .
    starting_positions = initial_position_vec
    isfirstrun = True
    if os . path . exists ( backend . filename ) :
        if backend . iteration > 1 :
            starting_positions = None
            isfirstrun = False
    if verbose and isfirstrun :
        LOGINFO ( 'start {:s} MCMC with {:d} dims, {:d} steps, {:d} walkers,' . format ( fittype , n_dim , n_mcmc_steps , n_walkers ) + ' {:d} threads' . format ( nworkers ) )
    elif verbose and not isfirstrun :
        LOGINFO ( 'continue {:s} with {:d} dims, {:d} steps, {:d} walkers, ' . format ( fittype , n_dim , n_mcmc_steps , n_walkers ) + '{:d} threads' . format ( nworkers ) )
    import sys
    if sys . version_info >= ( 3 , 3 ) :
        with Pool ( nworkers ) as pool :
            sampler = emcee . EnsembleSampler ( n_walkers , n_dim , log_posterior_transit , args = ( init_params , init_m , stimes , smags , serrs , priorbounds ) , pool = pool , backend = backend )
            sampler . run_mcmc ( starting_positions , n_mcmc_steps , progress = mcmcprogressbar )
    elif sys . version_info < ( 3 , 3 ) :
        sampler = emcee . EnsembleSampler ( n_walkers , n_dim , log_posterior_transit , args = ( init_params , init_m , stimes , smags , serrs , priorbounds ) , threads = nworkers , backend = backend )
        sampler . run_mcmc ( starting_positions , n_mcmc_steps , progress = mcmcprogressbar )
    if verbose :
        LOGINFO ( 'ended {:s} MCMC run with {:d} steps, {:d} walkers, ' . format ( fittype , n_mcmc_steps , n_walkers ) + '{:d} threads' . format ( nworkers ) )
reader = emcee . backends . HDFBackend ( samplesavpath )
n_to_discard = int ( burninpercent * n_mcmc_steps )
samples = reader . get_chain ( discard = n_to_discard , flat = True )
log_prob_samples = reader . get_log_prob ( discard = n_to_discard , flat = True )
log_prior_samples = reader . get_blobs ( discard = n_to_discard , flat = True )
# Get best - fit parameters and their 1 - sigma error bars
fit_statistics = list ( map ( lambda v : ( v [ 1 ] , v [ 2 ] - v [ 1 ] , v [ 1 ] - v [ 0 ] ) , list ( zip ( * np . percentile ( samples , [ 15.85 , 50 , 84.15 ] , axis = 0 ) ) ) ) )
medianparams , std_perrs , std_merrs = { } , { } , { }
for ix , k in enumerate ( np . sort ( list ( priorbounds . keys ( ) ) ) ) :
    medianparams [ k ] = fit_statistics [ ix ] [ 0 ]
    std_perrs [ k ] = fit_statistics [ ix ] [ 1 ]
    std_merrs [ k ] = fit_statistics [ ix ] [ 2 ]
stderrs = { 'std_perrs' : std_perrs , 'std_merrs' : std_merrs }
per = _get_value ( 'period' , medianparams , fixedparams )
t0 = _get_value ( 't0' , medianparams , fixedparams )
rp = _get_value ( 'rp' , medianparams , fixedparams )
sma = _get_value ( 'sma' , medianparams , fixedparams )
incl = _get_value ( 'incl' , medianparams , fixedparams )
ecc = _get_value ( 'ecc' , medianparams , fixedparams )
omega = _get_value ( 'omega' , medianparams , fixedparams )
limb_dark = _get_value ( 'limb_dark' , medianparams , fixedparams )
try :
    u = fixedparams [ 'u' ]
except Exception as e :
    u = [ medianparams [ 'u_linear' ] , medianparams [ 'u_quad' ] ]
fit_params , fit_m = _transit_model ( stimes , t0 , per , rp , sma , incl , ecc , omega , u , limb_dark )
fitmags = fit_m . light_curve ( fit_params )
fepoch = t0
# assemble the return dictionary
returndict = { 'fittype' : fittype , 'fitinfo' : { 'initialparams' : fitparams , 'initialmags' : init_flux , 'fixedparams' : fixedparams , 'finalparams' : medianparams , 'finalparamerrs' : stderrs , 'fitmags' : fitmags , 'fitepoch' : fepoch , } , 'fitplotfile' : None , 'magseries' : { 'times' : stimes , 'mags' : smags , 'errs' : serrs , 'magsarefluxes' : magsarefluxes , } , }
# make the output corner plot , and lightcurve plot if desired
if plotcorner :
    if isinstance ( trueparams , dict ) :
        trueparamkeys = np . sort ( list ( trueparams . keys ( ) ) )
        truelist = [ trueparams [ k ] for k in trueparamkeys ]
        fig = corner . corner ( samples , labels = trueparamkeys , truths = truelist , quantiles = [ 0.1585 , 0.5 , .8415 ] , show_titles = True )
    else :
        fig = corner . corner ( samples , labels = fitparamnames , quantiles = [ 0.1585 , 0.5 , .8415 ] , show_titles = True )
    plt . savefig ( plotcorner , dpi = 300 )
    if verbose :
        LOGINFO ( 'saved {:s}' . format ( plotcorner ) )
if plotfit and isinstance ( plotfit , str ) :
    f , ax = plt . subplots ( figsize = ( 8 , 6 ) )
    ax . scatter ( stimes , smags , c = 'k' , alpha = 0.5 , label = 'observed' , zorder = 1 , s = 1.5 , rasterized = True , linewidths = 0 )
    ax . scatter ( stimes , init_flux , c = 'r' , alpha = 1 , s = 3.5 , zorder = 2 , rasterized = True , linewidths = 0 , label = 'initial guess' )
    ax . scatter ( stimes , fitmags , c = 'b' , alpha = 1 , s = 1.5 , zorder = 3 , rasterized = True , linewidths = 0 , label = 'fit {:d} dims' . format ( len ( fitparamnames ) ) )
    ax . legend ( loc = 'best' )
    ax . set ( xlabel = 'time [days]' , ylabel = 'relative flux' )
    f . savefig ( plotfit , dpi = 300 , bbox_inches = 'tight' )
    if verbose :
        LOGINFO ( 'saved {:s}' . format ( plotfit ) )
    returndict [ 'fitplotfile' ] = plotfit
return returndict


signature:
def mandelagol_and_line_fit_magseries ( times , mags , errs , fitparams , priorbounds , fixedparams , trueparams = None , burninpercent = 0.3 , plotcorner = False , timeoffset = 0 , samplesavpath = False , n_walkers = 50 , n_mcmc_steps = 400 , eps = 1e-4 , skipsampling = False , overwriteexistingsamples = False , mcmcprogressbar = False , plotfit = False , scatterxdata = None , scatteryaxes = None , magsarefluxes = True , sigclip = 10.0 , verbose = True , nworkers = 4 ) :
    '''The model fit by this function is : a Mandel & Agol ( 2002 ) transit , PLUS a
 line . You can fit and fix whatever parameters you want .
 A typical use case : you want to measure transit times of individual SNR > ~
 50 transits . You fix all the transit parameters except for the mid - time ,
 and also fit for a line locally .
 NOTE : this only works for flux time - series at the moment .
 NOTE : Between the ` fitparams ` , ` priorbounds ` , and ` fixedparams ` dicts , you
 must specify all of the planetary transit parameters required by BATMAN and
 the parameters for the line fit : ` [ ' t0 ' , ' rp ' , ' sma ' , ' incl ' , ' u ' , ' rp ' ,
 ' ecc ' , ' omega ' , ' period ' , ' poly _ order0 ' , poly _ order1 ' ] ` , or the BATMAN model
 will fail to initialize .
 Parameters
 times , mags , errs : np . array
 The input flux time - series to fit a Fourier cosine series to .
 fitparams : dict
 This is the initial parameter guesses for MCMC , found e . g . , by
 BLS . The key string format must not be changed , but any parameter can be
 either " fit " or " fixed " . If it is " fit " , it must have a corresponding
 prior . For example : :
 fitparams = { ' t0 ' : 1325.9,
 ' poly _ order0 ' : 1,
 ' poly _ order1 ' : 0 . }
 where ` t0 ` is the time of transit - center for a reference transit .
 ` poly _ order0 ` corresponds to the intercept of the line , ` poly _ order1 ` is
 the slope .
 priorbounds : dict
 This sets the lower & upper bounds on uniform prior , e . g . : :
 priorbounds = { ' t0 ' : ( np . min ( time ) , np . max ( time ) ) ,
 ' poly _ order0 ' : ( 0.5,1.5 ) ,
 ' poly _ order1 ' : ( - 0.5,0.5 ) }
 fixedparams : dict
 This sets which parameters are fixed , and their values . For example : :
 fixedparams = { ' ecc ' : 0 . ,
 ' omega ' : 90 . ,
 ' limb _ dark ' : ' quadratic ' ,
 ' period ' : fitd [ ' period ' ] ,
 ' rp ' : np . sqrt ( fitd [ ' transitdepth ' ] ) ,
 ' sma ' : 6.17 , ' incl ' : 85 , ' u ' : [ 0.3 , 0.2 ] }
 ` limb _ dark ` must be " quadratic " . It ' s " fixed " , because once you
 choose your limb - darkening model , it ' s fixed .
 trueparams : list of floats
 The true parameter values you ' re fitting for , if they ' re known ( e . g . , a
 known planet , or fake data ) . Only for plotting purposes .
 burninpercent : float
 The percent of MCMC samples to discard as burn - in .
 plotcorner : str or False
 If this is a str , points to the path of output corner plot that will be
 generated for this MCMC run .
 timeoffset : float
 If input times are offset by some constant , and you want saved pickles
 to fix that .
 samplesavpath : str
 This must be provided so ` emcee ` can save its MCMC samples to disk as
 HDF5 files . This will set the path of the output HDF5file written .
 n _ walkers : int
 The number of MCMC walkers to use .
 n _ mcmc _ steps : int
 The number of MCMC steps to take .
 eps : float
 The radius of the ` n _ walkers - dimensional ` Gaussian ball used to
 initialize the MCMC .
 skipsampling : bool
 If you ' ve already collected MCMC samples , and you do not want any more
 sampling ( e . g . , just make the plots ) , set this to be True .
 overwriteexistingsamples : bool
 If you ' ve collected samples , but you want to overwrite them , set this to
 True . Usually , it should be False , which appends samples to
 ` samplesavpath ` HDF5 file .
 mcmcprogressbar : bool
 If True , will show a progress bar for the MCMC process .
 plotfit : str or bool
 If a str , indicates the path of the output fit plot file . If False , no
 fit plot will be made .
 scatterxdata : np . array or None
 Use this to overplot x , y scatter points on the output model / data
 lightcurve ( e . g . , to highlight bad data , or to indicate an ephemeris ) ,
 this can take a ` np . ndarray ` with the same units as ` times ` .
 scatteryaxes : np . array or None
 Use this to provide the y - values for scatterxdata , in units of fraction
 of an axis .
 magsarefluxes : bool
 This indicates if the input measurements in ` mags ` are actually fluxes .
 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 .
 verbose : bool
 If True , will indicate MCMC progress .
 nworkers : int
 The number of parallel workers to launch for MCMC .
 Returns
 dict
 This function returns a dict containing the model fit parameters and
 other fit information . The form of this dict is mostly standardized
 across all functions in this module : :
 ' fittype ' : ' mandelagol _ and _ line ' ,
 ' fitinfo ' : {
 ' initialparams ' : the initial transit params provided ,
 ' fixedparams ' : the fixed transit params provided ,
 ' finalparams ' : the final model fit transit params ,
 ' finalparamerrs ' : formal errors in the params ,
 ' fitmags ' : the model fit mags ,
 ' fitepoch ' : the epoch of minimum light for the fit ,
 ' fitplotfile ' : the output fit plot if fitplot is not None ,
 ' magseries ' : {
 ' times ' : input times in phase order of the model ,
 ' phase ' : the phases of the model mags ,
 ' mags ' : input mags / fluxes in the phase order of the model ,
 ' errs ' : errs in the phase order of the model ,
 ' magsarefluxes ' : input value of magsarefluxes kwarg'''

implementation:
from multiprocessing import Pool
fittype = 'mandelagol_and_line'
if not magsarefluxes :
    raise NotImplementedError ( 'magsarefluxes is not implemented yet.' )
if not samplesavpath :
    raise ValueError ( 'This function requires that you save the samples somewhere' )
if not mandel_agol_dependencies :
    raise ImportError ( 'This function depends on BATMAN, emcee>3.0, corner, and h5py.' )
# sigma clip and get rid of zero errs
stimes , smags , serrs = sigclip_magseries ( times , mags , errs , sigclip = sigclip , magsarefluxes = magsarefluxes )
nzind = np . nonzero ( serrs )
stimes , smags , serrs = stimes [ nzind ] , smags [ nzind ] , serrs [ nzind ]
init_period = _get_value ( 'period' , fitparams , fixedparams )
init_epoch = _get_value ( 't0' , fitparams , fixedparams )
init_rp = _get_value ( 'rp' , fitparams , fixedparams )
init_sma = _get_value ( 'sma' , fitparams , fixedparams )
init_incl = _get_value ( 'incl' , fitparams , fixedparams )
init_ecc = _get_value ( 'ecc' , fitparams , fixedparams )
init_omega = _get_value ( 'omega' , fitparams , fixedparams )
limb_dark = _get_value ( 'limb_dark' , fitparams , fixedparams )
init_u = _get_value ( 'u' , fitparams , fixedparams )
init_poly_order0 = _get_value ( 'poly_order0' , fitparams , fixedparams )
init_poly_order1 = _get_value ( 'poly_order1' , fitparams , fixedparams )
if not limb_dark == 'quadratic' :
    raise ValueError ( 'only quadratic limb-darkening is supported at the moment' )
# initialize the model and calculate the initial model light - curve
init_params , init_m = _transit_model ( stimes , init_epoch , init_period , init_rp , init_sma , init_incl , init_ecc , init_omega , init_u , limb_dark )
init_flux = ( init_m . light_curve ( init_params ) + init_poly_order0 + init_poly_order1 * stimes )
# guessed initial params . give nice guesses , or else emcee struggles .
theta , fitparamnames = [ ] , [ ]
for k in np . sort ( list ( fitparams . keys ( ) ) ) :
    if isinstance ( fitparams [ k ] , float ) or isinstance ( fitparams [ k ] , int ) :
        theta . append ( fitparams [ k ] )
        fitparamnames . append ( fitparams [ k ] )
    elif isinstance ( fitparams [ k ] , list ) :
        if not len ( fitparams [ k ] ) == 2 :
            raise ValueError ( 'should only be quadratic LD coeffs' )
        theta . append ( fitparams [ k ] [ 0 ] )
        theta . append ( fitparams [ k ] [ 1 ] )
        fitparamnames . append ( fitparams [ k ] [ 0 ] )
        fitparamnames . append ( fitparams [ k ] [ 1 ] )
# initialize sampler
n_dim = len ( theta )
# run the MCMC , unless you just want to load the available samples
if not skipsampling :
    backend = emcee . backends . HDFBackend ( samplesavpath )
    if overwriteexistingsamples :
        LOGWARNING ( 'erased samples previously at {:s}' . format ( samplesavpath ) )
        backend . reset ( n_walkers , n_dim )
    # if this is the first run , then start from a gaussian ball , centered
    # on the maximum likelihood solution . otherwise , resume from the
    # previous samples .
    def nll ( * args ) :
        return - _log_likelihood_transit_plus_line ( * args )
    soln = spminimize ( nll , theta , method = 'BFGS' , args = ( init_params , init_m , stimes , smags , serrs , priorbounds ) )
    theta_ml = soln . x
    ml_poly_order0 = theta_ml [ 0 ]
    ml_poly_order1 = theta_ml [ 1 ]
    ml_rp = theta_ml [ 2 ]
    ml_t0 = theta_ml [ 3 ]
    ml_params , ml_m = _transit_model ( stimes , ml_t0 , init_period , ml_rp , init_sma , init_incl , init_ecc , init_omega , init_u , limb_dark )
    ml_mags = ( ml_m . light_curve ( ml_params ) + ml_poly_order0 + ml_poly_order1 * stimes )
    initial_position_vec = [ theta_ml + eps * np . random . randn ( n_dim ) for i in range ( n_walkers ) ]
    starting_positions = initial_position_vec
    isfirstrun = True
    if os . path . exists ( backend . filename ) :
        if backend . iteration > 1 :
            starting_positions = None
            isfirstrun = False
    if verbose and isfirstrun :
        LOGINFO ( 'start {:s} MCMC with {:d} dims, {:d} steps, {:d} walkers,' . format ( fittype , n_dim , n_mcmc_steps , n_walkers ) + ' {:d} threads' . format ( nworkers ) )
    elif verbose and not isfirstrun :
        LOGINFO ( 'continue {:s} with {:d} dims, {:d} steps, {:d} walkers, ' . format ( fittype , n_dim , n_mcmc_steps , n_walkers ) + '{:d} threads' . format ( nworkers ) )
    with Pool ( nworkers ) as pool :
        sampler = emcee . EnsembleSampler ( n_walkers , n_dim , log_posterior_transit_plus_line , args = ( init_params , init_m , stimes , smags , serrs , priorbounds ) , pool = pool , backend = backend )
        sampler . run_mcmc ( starting_positions , n_mcmc_steps , progress = mcmcprogressbar )
    if verbose :
        LOGINFO ( 'ended {:s} MCMC run with {:d} steps, {:d} walkers, ' . format ( fittype , n_mcmc_steps , n_walkers ) + '{:d} threads' . format ( nworkers ) )
reader = emcee . backends . HDFBackend ( samplesavpath )
n_to_discard = int ( burninpercent * n_mcmc_steps )
samples = reader . get_chain ( discard = n_to_discard , flat = True )
log_prob_samples = reader . get_log_prob ( discard = n_to_discard , flat = True )
log_prior_samples = reader . get_blobs ( discard = n_to_discard , flat = True )
# Get best - fit parameters and their 1 - sigma error bars
fit_statistics = list ( map ( lambda v : ( v [ 1 ] , v [ 2 ] - v [ 1 ] , v [ 1 ] - v [ 0 ] ) , list ( zip ( * np . percentile ( samples , [ 15.85 , 50 , 84.15 ] , axis = 0 ) ) ) ) )
medianparams , std_perrs , std_merrs = { } , { } , { }
for ix , k in enumerate ( np . sort ( list ( priorbounds . keys ( ) ) ) ) :
    medianparams [ k ] = fit_statistics [ ix ] [ 0 ]
    std_perrs [ k ] = fit_statistics [ ix ] [ 1 ]
    std_merrs [ k ] = fit_statistics [ ix ] [ 2 ]
stderrs = { 'std_perrs' : std_perrs , 'std_merrs' : std_merrs }
per = _get_value ( 'period' , medianparams , fixedparams )
t0 = _get_value ( 't0' , medianparams , fixedparams )
rp = _get_value ( 'rp' , medianparams , fixedparams )
sma = _get_value ( 'sma' , medianparams , fixedparams )
incl = _get_value ( 'incl' , medianparams , fixedparams )
ecc = _get_value ( 'ecc' , medianparams , fixedparams )
omega = _get_value ( 'omega' , medianparams , fixedparams )
limb_dark = _get_value ( 'limb_dark' , medianparams , fixedparams )
try :
    u = fixedparams [ 'u' ]
except Exception as e :
    u = [ medianparams [ 'u_linear' ] , medianparams [ 'u_quad' ] ]
poly_order0 = _get_value ( 'poly_order0' , medianparams , fixedparams )
poly_order1 = _get_value ( 'poly_order1' , medianparams , fixedparams )
# initialize the model and calculate the initial model light - curve
fit_params , fit_m = _transit_model ( stimes , t0 , per , rp , sma , incl , ecc , omega , u , limb_dark )
fitmags = ( fit_m . light_curve ( fit_params ) + poly_order0 + poly_order1 * stimes )
fepoch = t0
# assemble the return dictionary
medianparams [ 't0' ] += timeoffset
returndict = { 'fittype' : fittype , 'fitinfo' : { 'initialparams' : fitparams , 'initialmags' : init_flux , 'fixedparams' : fixedparams , 'finalparams' : medianparams , 'finalparamerrs' : stderrs , 'fitmags' : fitmags , 'fitepoch' : fepoch + timeoffset , } , 'fitplotfile' : None , 'magseries' : { 'times' : stimes + timeoffset , 'mags' : smags , 'errs' : serrs , 'magsarefluxes' : magsarefluxes , } , }
# make the output corner plot , and lightcurve plot if desired
if plotcorner :
    fig = corner . corner ( samples , labels = [ 'line intercept-1' , 'line slope' , 'rp' , 't0-{:.4f}' . format ( timeoffset ) ] , truths = [ ml_poly_order0 , ml_poly_order1 , ml_rp , ml_t0 ] , quantiles = [ 0.1585 , 0.5 , .8415 ] , show_titles = True )
    plt . savefig ( plotcorner , dpi = 300 )
    if verbose :
        LOGINFO ( 'saved {:s}' . format ( plotcorner ) )
if plotfit and isinstance ( plotfit , str ) :
    plt . close ( 'all' )
    f , ( a0 , a1 ) = plt . subplots ( nrows = 2 , ncols = 1 , sharex = True , figsize = ( 8 , 5 ) , gridspec_kw = { 'height_ratios' : [ 3 , 1 ] } )
    a0 . scatter ( stimes , smags , c = 'k' , alpha = 0.9 , label = 'data' , zorder = 1 , s = 10 , rasterized = True , linewidths = 0 )
    DEBUGGING = False
    if DEBUGGING :
        a0 . scatter ( stimes , init_flux , c = 'r' , alpha = 1 , s = 3.5 , zorder = 2 , rasterized = True , linewidths = 0 , label = 'initial guess for ml' )
        a0 . scatter ( stimes , ml_mags , c = 'g' , alpha = 1 , s = 3.5 , zorder = 2 , rasterized = True , linewidths = 0 , label = 'max likelihood' )
    a0 . plot ( stimes , fitmags , c = 'b' , zorder = 0 , rasterized = True , lw = 2 , alpha = 0.4 , label = '{:s} fit, {:d} dims' . format ( fittype , len ( fitparamnames ) ) )
    a1 . scatter ( stimes , smags - fitmags , c = 'k' , alpha = 0.9 , rasterized = True , s = 10 , linewidths = 0 )
    if scatterxdata and scatteryaxes :
        import matplotlib . transforms as transforms
        for a in [ a0 , a1 ] :
            transform = transforms . blended_transform_factory ( a . transData , a . transAxes )
            a . scatter ( scatterxdata , scatteryaxes , c = 'r' , alpha = 0.9 , zorder = 2 , s = 10 , rasterized = True , linewidths = 0 , marker = "^" , transform = transform )
    a1 . set_xlabel ( 'time-t0 [days]' )
    a0 . set_ylabel ( 'relative flux' )
    a1 . set_ylabel ( 'residual' )
    a0 . legend ( loc = 'best' , fontsize = 'x-small' )
    for a in [ a0 , a1 ] :
        a . get_yaxis ( ) . set_tick_params ( which = 'both' , direction = 'in' )
        a . get_xaxis ( ) . set_tick_params ( which = 'both' , direction = 'in' )
    f . tight_layout ( h_pad = 0 , w_pad = 0 )
    f . savefig ( plotfit , dpi = 300 , bbox_inches = 'tight' )
    if verbose :
        LOGINFO ( 'saved {:s}' . format ( plotfit ) )
    returndict [ 'fitplotfile' ] = plotfit
return returndict


signature:
def list_trilegal_filtersystems ( ) :
    '''This just lists all the filter systems available for TRILEGAL .'''

implementation:
print ( '%-40s %s' % ( 'FILTER SYSTEM NAME' , 'DESCRIPTION' ) )
print ( '%-40s %s' % ( '------------------' , '-----------' ) )
for key in sorted ( TRILEGAL_FILTER_SYSTEMS . keys ( ) ) :
    print ( '%-40s %s' % ( key , TRILEGAL_FILTER_SYSTEMS [ key ] [ 'desc' ] ) )


signature:
def query_galcoords ( gal_lon , gal_lat , filtersystem = 'sloan_2mass' , field_deg2 = 1.0 , usebinaries = True , extinction_sigma = 0.1 , magnitude_limit = 26.0 , maglim_filtercol = 4 , trilegal_version = 1.6 , extraparams = None , forcefetch = False , cachedir = '~/.astrobase/trilegal-cache' , verbose = True , timeout = 60.0 , refresh = 150.0 , maxtimeout = 700.0 ) :
    '''This queries the TRILEGAL model form , downloads results , and parses them .
 Parameters
 gal _ lon , gal _ lat : float
 These are the center galactic longitude and latitude in degrees .
 filtersystem : str
 This is a key in the TRILEGAL _ FILTER _ SYSTEMS dict . Use the function
 : py : func : ` astrobase . services . trilegal . list _ trilegal _ filtersystems ` to
 see a nicely formatted table with the key and description for each of
 these .
 field _ deg2 : float
 The area of the simulated field in square degrees .
 usebinaries : bool
 If this is True , binaries will be present in the model results .
 extinction _ sigma : float
 This is the applied std dev around the ` Av _ extinction ` value for the
 galactic coordinates requested .
 magnitude _ limit : float
 This is the limiting magnitude of the simulation in the
 ` maglim _ filtercol ` band index of the filter system chosen .
 maglim _ filtercol : int
 The index in the filter system list of the magnitude limiting band .
 trilegal _ version : float
 This is the the version of the TRILEGAL form to use . This can usually be
 left as - is .
 extraparams : dict or None
 This is a dict that can be used to override parameters of the model
 other than the basic ones used for input to this function . All
 parameters are listed in ` TRILEGAL _ DEFAULT _ PARAMS ` above . See :
 http : / / stev . oapd . inaf . it / cgi - bin / trilegal
 for explanations of these parameters .
 forcefetch : bool
 If this is True , the query will be retried even if cached results for
 it exist .
 cachedir : str
 This points to the directory where results will be downloaded .
 verbose : bool
 If True , will indicate progress and warn of any issues .
 timeout : float
 This sets the amount of time in seconds to wait for the service to
 respond to our initial request .
 refresh : float
 This sets the amount of time in seconds to wait before checking if the
 result file is available . If the results file isn ' t available after
 ` refresh ` seconds have elapsed , the function will wait for ` refresh `
 seconds continuously , until ` maxtimeout ` is reached or the results file
 becomes available .
 maxtimeout : float
 The maximum amount of time in seconds to wait for a result to become
 available after submitting our query request .
 Returns
 dict
 This returns a dict of the form : :
 { ' params ' : the input param dict used ,
 ' extraparams ' : any extra params used ,
 ' provenance ' : ' cached ' or ' new download ' ,
 ' tablefile ' : the path on disk to the downloaded model text file }'''

implementation:
# these are the default parameters
inputparams = copy . deepcopy ( TRILEGAL_INPUT_PARAMS )
# update them with the input params
inputparams [ 'binary_kind' ] = '1' if usebinaries else '0'
inputparams [ 'extinction_sigma' ] = '%.2f' % extinction_sigma
inputparams [ 'field' ] = '%.2f' % field_deg2
inputparams [ 'icm_lim' ] = str ( maglim_filtercol )
inputparams [ 'mag_lim' ] = '%.2f' % magnitude_limit
inputparams [ 'trilegal_version' ] = str ( trilegal_version )
# get the coordinates
inputparams [ 'gc_l' ] = '%.3f' % gal_lon
inputparams [ 'gc_b' ] = '%.3f' % gal_lat
# check if the area is less than 10 deg ^ 2
if field_deg2 > 10.0 :
    LOGERROR ( "can't have an area > 10 square degrees" )
    return None
# get the extinction parameter . this is by default A [ inf ] in V . we ' ll use
# the value from SF11 generated by the 2MASS DUST service
extinction_info = dust . extinction_query ( gal_lon , gal_lat , coordtype = 'galactic' , forcefetch = forcefetch , verbose = verbose , timeout = timeout )
try :
    Av_infinity = extinction_info [ 'Amag' ] [ 'CTIO V' ] [ 'sf11' ]
    inputparams [ 'extinction_infty' ] = '%.5f' % Av_infinity
except Exception as e :
    LOGEXCEPTION ( 'could not get A_V_SF11 from 2MASS DUST ' 'for Galactic coords: (%.3f, %.3f), ' 'using default value of %s' % ( gal_lon , gal_lat , inputparams [ 'extinction_infty' ] ) )
# get the filter system table
if filtersystem in TRILEGAL_FILTER_SYSTEMS :
    inputparams [ 'photsys_file' ] = ( TRILEGAL_FILTER_SYSTEMS [ filtersystem ] [ 'table' ] )
else :
    LOGERROR ( 'filtersystem name: %s is not in the table of known ' 'filter systems.\n' 'Try the trilegal.list_trilegal_filtersystems() function ' 'to see all available filter systems.' % filtersystem )
    return None
# override the complete form param dict now with our params
trilegal_params = copy . deepcopy ( TRILEGAL_DEFAULT_PARAMS )
trilegal_params . update ( inputparams )
# override the final params with any extraparams
if extraparams and isinstance ( extraparams , dict ) :
    trilegal_params . update ( extraparams )
# see if the cachedir exists
if '~' in cachedir :
    cachedir = os . path . expanduser ( cachedir )
if not os . path . exists ( cachedir ) :
    os . makedirs ( cachedir )
# generate the cachefname and look for it
cachekey = repr ( inputparams )
cachekey = hashlib . sha256 ( cachekey . encode ( ) ) . hexdigest ( )
cachefname = os . path . join ( cachedir , '%s.txt.gz' % cachekey )
provenance = 'cache'
lockfile = os . path . join ( cachedir , 'LOCK-%s' % cachekey )
# run the query if results not found in the cache
if forcefetch or ( not os . path . exists ( cachefname ) ) : # first , check if a query like this is running already
    if os . path . exists ( lockfile ) :
        with open ( lockfile , 'r' ) as infd :
            lock_contents = infd . read ( )
        lock_contents = lock_contents . replace ( '\n' , '' )
        LOGERROR ( 'this query appears to be active since %s' 'in another instance, not running it again' % lock_contents )
        return None
    else :
        with open ( lockfile , 'w' ) as outfd :
            outfd . write ( datetime . utcnow ( ) . isoformat ( ) )
    provenance = 'new download'
    try :
        if verbose :
            LOGINFO ( 'submitting TRILEGAL request for input params: %s' % repr ( inputparams ) )
        posturl = TRILEGAL_POSTURL . format ( formversion = trilegal_version )
        req = requests . post ( posturl , data = trilegal_params , timeout = timeout )
        resp = req . text
        # get the URL of the result file
        resultfile = TRILEGAL_REGEX . search ( resp )
        if resultfile :
            resultfile = resultfile [ 0 ]
            waitdone = False
            timeelapsed = 0.0
            resultfileurl = '%s/%s' % ( TRILEGAL_BASEURL , resultfile . replace ( 'a href=..' , '' ) )
            if verbose :
                LOGINFO ( 'request submitted sucessfully, waiting for results...' )
            # wait for 2 minutes , then try to download the result file
            while not waitdone :
                if timeelapsed > maxtimeout :
                    LOGERROR ( 'TRILEGAL timed out after waiting for results,' ' request was: ' '%s' % repr ( inputparams ) )
                    # remove the lock file
                    if os . path . exists ( lockfile ) :
                        os . remove ( lockfile )
                    return None
                time . sleep ( refresh )
                timeelapsed = timeelapsed + refresh
                try :
                    resreq = requests . get ( resultfileurl )
                    resreq . raise_for_status ( )
                    if verbose :
                        LOGINFO ( 'TRILEGAL completed, retrieving results...' )
                    # stream the response to the output cache file
                    with gzip . open ( cachefname , 'wb' ) as outfd :
                        for chunk in resreq . iter_content ( chunk_size = 65536 ) :
                            outfd . write ( chunk )
                    tablefname = cachefname
                    waitdone = True
                    if verbose :
                        LOGINFO ( 'done.' )
                except Exception as e :
                    if verbose :
                        LOGINFO ( 'elapsed time: %.1f, result file: %s ' 'not ready yet...' % ( timeelapsed , resultfileurl ) )
                    continue
        else :
            LOGERROR ( 'no result file URL found in TRILEGAL output, ' 'this is probably an error with the input. ' 'HTML of error page follows:\n' )
            LOGINFO ( resp )
            # remove the lock file
            if os . path . exists ( lockfile ) :
                os . remove ( lockfile )
            return None
    except requests . exceptions . Timeout as e :
        LOGERROR ( 'TRILEGAL submission timed out, ' 'site is probably down. Request was: ' '%s' % repr ( inputparams ) )
        return None
    except Exception as e :
        LOGEXCEPTION ( 'TRILEGAL request failed for ' '%s' % repr ( inputparams ) )
        return None
    finally : # remove the lock file
        if os . path . exists ( lockfile ) :
            os . remove ( lockfile )
# otherwise , get the file from the cache
else :
    if verbose :
        LOGINFO ( 'getting cached TRILEGAL model result for ' 'request: %s' % ( repr ( inputparams ) ) )
    tablefname = cachefname
# return a dict pointing to the result file
# we ' ll parse this later
resdict = { 'params' : inputparams , 'extraparams' : extraparams , 'provenance' : provenance , 'tablefile' : tablefname }
return resdict


signature:
def query_radecl ( ra , decl , filtersystem = 'sloan_2mass' , field_deg2 = 1.0 , usebinaries = True , extinction_sigma = 0.1 , magnitude_limit = 26.0 , maglim_filtercol = 4 , trilegal_version = 1.6 , extraparams = None , forcefetch = False , cachedir = '~/.astrobase/trilegal-cache' , verbose = True , timeout = 60.0 , refresh = 150.0 , maxtimeout = 700.0 ) :
    '''This runs the TRILEGAL query for decimal equatorial coordinates .
 Parameters
 ra , decl : float
 These are the center equatorial coordinates in decimal degrees
 filtersystem : str
 This is a key in the TRILEGAL _ FILTER _ SYSTEMS dict . Use the function
 : py : func : ` astrobase . services . trilegal . list _ trilegal _ filtersystems ` to
 see a nicely formatted table with the key and description for each of
 these .
 field _ deg2 : float
 The area of the simulated field in square degrees . This is in the
 Galactic coordinate system .
 usebinaries : bool
 If this is True , binaries will be present in the model results .
 extinction _ sigma : float
 This is the applied std dev around the ` Av _ extinction ` value for the
 galactic coordinates requested .
 magnitude _ limit : float
 This is the limiting magnitude of the simulation in the
 ` maglim _ filtercol ` band index of the filter system chosen .
 maglim _ filtercol : int
 The index in the filter system list of the magnitude limiting band .
 trilegal _ version : float
 This is the the version of the TRILEGAL form to use . This can usually be
 left as - is .
 extraparams : dict or None
 This is a dict that can be used to override parameters of the model
 other than the basic ones used for input to this function . All
 parameters are listed in ` TRILEGAL _ DEFAULT _ PARAMS ` above . See :
 http : / / stev . oapd . inaf . it / cgi - bin / trilegal
 for explanations of these parameters .
 forcefetch : bool
 If this is True , the query will be retried even if cached results for
 it exist .
 cachedir : str
 This points to the directory where results will be downloaded .
 verbose : bool
 If True , will indicate progress and warn of any issues .
 timeout : float
 This sets the amount of time in seconds to wait for the service to
 respond to our initial request .
 refresh : float
 This sets the amount of time in seconds to wait before checking if the
 result file is available . If the results file isn ' t available after
 ` refresh ` seconds have elapsed , the function will wait for ` refresh `
 seconds continuously , until ` maxtimeout ` is reached or the results file
 becomes available .
 maxtimeout : float
 The maximum amount of time in seconds to wait for a result to become
 available after submitting our query request .
 Returns
 dict
 This returns a dict of the form : :
 { ' params ' : the input param dict used ,
 ' extraparams ' : any extra params used ,
 ' provenance ' : ' cached ' or ' new download ' ,
 ' tablefile ' : the path on disk to the downloaded model text file }'''

implementation:
# convert the ra / decl to gl , gb
radecl = SkyCoord ( ra = ra * u . degree , dec = decl * u . degree )
gl = radecl . galactic . l . degree
gb = radecl . galactic . b . degree
return query_galcoords ( gl , gb , filtersystem = filtersystem , field_deg2 = field_deg2 , usebinaries = usebinaries , extinction_sigma = extinction_sigma , magnitude_limit = magnitude_limit , maglim_filtercol = maglim_filtercol , trilegal_version = trilegal_version , extraparams = extraparams , forcefetch = forcefetch , cachedir = cachedir , verbose = verbose , timeout = timeout , refresh = refresh , maxtimeout = maxtimeout )


signature:
def read_model_table ( modelfile ) :
    '''This reads a downloaded TRILEGAL model file .
 Parameters
 modelfile : str
 Path to the downloaded model file to read .
 Returns
 np . recarray
 Returns the model table as a Numpy record array .'''

implementation:
infd = gzip . open ( modelfile )
model = np . genfromtxt ( infd , names = True )
infd . close ( )
return model


signature:
def _time_independent_equals ( a , b ) :
    '''This compares two values in constant time .
 Taken from tornado :
 https : / / github . com / tornadoweb / tornado / blob /
 d4eb8eb4eb5cc9a6677e9116ef84ded8efba8859 / tornado / web . py # L3060'''

implementation:
if len ( a ) != len ( b ) :
    return False
result = 0
if isinstance ( a [ 0 ] , int ) : # python3 byte strings
    for x , y in zip ( a , b ) :
        result |= x ^ y
else : # python2
    for x , y in zip ( a , b ) :
        result |= ord ( x ) ^ ord ( y )
return result == 0


signature:
def default ( self , obj ) :
    '''Overrides the default serializer for ` JSONEncoder ` .
 This can serialize the following objects in addition to what
 ` JSONEncoder ` can already do .
 - ` np . array `
 - ` bytes `
 - ` complex `
 - ` np . float64 ` and other ` np . dtype ` objects
 Parameters
 obj : object
 A Python object to serialize to JSON .
 Returns
 str
 A JSON encoded representation of the input object .'''

implementation:
if isinstance ( obj , np . ndarray ) :
    return obj . tolist ( )
elif isinstance ( obj , bytes ) :
    return obj . decode ( )
elif isinstance ( obj , complex ) :
    return ( obj . real , obj . imag )
elif ( isinstance ( obj , ( float , np . float64 , np . float_ ) ) and not np . isfinite ( obj ) ) :
    return None
elif isinstance ( obj , ( np . int8 , np . int16 , np . int32 , np . int64 ) ) :
    return int ( obj )
else :
    return json . JSONEncoder . default ( self , obj )


signature:
def initialize ( self , currentdir , assetpath , cplist , cplistfile , executor , readonly , baseurl ) :
    '''handles initial setup .'''

implementation:
self . currentdir = currentdir
self . assetpath = assetpath
self . currentproject = cplist
self . cplistfile = cplistfile
self . executor = executor
self . readonly = readonly
self . baseurl = baseurl


signature:
def get ( self ) :
    '''This handles GET requests to the index page .
 TODO : provide the correct baseurl from the checkplotserver options dict ,
 so the frontend JS can just read that off immediately .'''

implementation:
# generate the project ' s list of checkplots
project_checkplots = self . currentproject [ 'checkplots' ]
project_checkplotbasenames = [ os . path . basename ( x ) for x in project_checkplots ]
project_checkplotindices = range ( len ( project_checkplots ) )
# get the sortkey and order
project_cpsortkey = self . currentproject [ 'sortkey' ]
if self . currentproject [ 'sortorder' ] == 'asc' :
    project_cpsortorder = 'ascending'
elif self . currentproject [ 'sortorder' ] == 'desc' :
    project_cpsortorder = 'descending'
# get the filterkey and condition
project_cpfilterstatements = self . currentproject [ 'filterstatements' ]
self . render ( 'cpindex.html' , project_checkplots = project_checkplots , project_cpsortorder = project_cpsortorder , project_cpsortkey = project_cpsortkey , project_cpfilterstatements = project_cpfilterstatements , project_checkplotbasenames = project_checkplotbasenames , project_checkplotindices = project_checkplotindices , project_checkplotfile = self . cplistfile , readonly = self . readonly , baseurl = self . baseurl )