astrobase.lcfit.transits module

Fitting routines for planetary transits:

astrobase.lcfit.transits.traptransit_fit_magseries(times, mags, errs, transitparams, param_bounds=None, scale_errs_redchisq_unity=True, sigclip=10.0, plotfit=False, magsarefluxes=False, verbose=True, curve_fit_kwargs=None)[source]

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 = [transit_period (time),
                     transit_epoch (time),
                     transit_depth (flux or mags),
                     transit_duration (phase),
                     ingress_duration (phase)]
    
    • for magnitudes -> transit_depth should be < 0
    • for fluxes -> transit_depth 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.

  • param_bounds (dict or None) –

    This is a dict of the upper and lower bounds on each fit parameter. Should be of the form:

    {'period':          (lower_bound_period, upper_bound_period),
     'epoch':           (lower_bound_epoch, upper_bound_epoch),
     'depth':           (lower_bound_depth, upper_bound_depth),
     'duration':        (lower_bound_duration, upper_bound_duration),
     'ingressduration': (lower_bound_ingressduration,
                         upper_bound_ingressduration)}
    
    • To indicate that a parameter is fixed, use ‘fixed’ instead of a tuple providing its lower and upper bounds as tuple.
    • To indicate that a parameter has no bounds, don’t include it in the param_bounds dict.

    If this is None, the default value of this kwarg will be:

    {'period':(0.0,np.inf),       # period is between 0 and inf
     'epoch':(0.0, np.inf),       # epoch is between 0 and inf
     'depth':(-np.inf,np.inf),    # depth is between -np.inf and np.inf
     'duration':(0.0,1.0),        # duration is between 0.0 and 1.0
     'ingressduration':(0.0,0.5)} # ingress duration between 0.0 and 0.5
    
  • scale_errs_redchisq_unity (bool) – If True, the standard errors on the fit parameters will be scaled to make the reduced chi-sq = 1.0. This sets the absolute_sigma kwarg for the scipy.optimize.curve_fit function to False.
  • 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.
  • curve_fit_kwargs (dict or None) – If not None, this should be a dict containing extra kwargs to pass to the scipy.optimize.curve_fit function.
Returns:

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,
        '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
    }
}

Return type:

dict

astrobase.lcfit.transits.log_posterior_transit(theta, params, model, t, flux, err_flux, priorbounds)[source]

Evaluate posterior probability given proposed model parameters and the observed flux timeseries.

astrobase.lcfit.transits.log_posterior_transit_plus_line(theta, params, model, t, flux, err_flux, priorbounds)[source]

Evaluate posterior probability given proposed model parameters and the observed flux timeseries.

astrobase.lcfit.transits.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, exp_time_minutes=2, eps=0.0001, skipsampling=False, overwriteexistingsamples=False, mcmcprogressbar=False, plotfit=False, magsarefluxes=False, sigclip=10.0, verbose=True, nworkers=4)[source]

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.
  • exp_time_minutes (int) – Exposure time, in minutes, passed to transit model to smear observations.
  • 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:

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,
        'acceptancefraction': fraction of MCMC ensemble. low=bad.
        'autocorrtime': if autocorrtime ~= n_mcmc_steps, not good.
    },
    '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
    }
}

Return type:

dict

astrobase.lcfit.transits.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, exp_time_minutes=2, eps=0.0001, skipsampling=False, overwriteexistingsamples=False, mcmcprogressbar=False, plotfit=False, scatterxdata=None, scatteryaxes=None, magsarefluxes=True, sigclip=10.0, verbose=True, nworkers=4)[source]

The model fit by this function is: a Mandel & Agol (2002) transit, PLUS a line. You can fit and fix whatever parameters you want.

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.
  • exp_time_minutes (int) – Exposure time, in minutes, passed to transit model to smear observations.
  • 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:

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,
        'acceptancefraction': fraction of MCMC ensemble. low=bad.
        'autocorrtime': if autocorrtime ~= n_mcmc_steps, not good.
    },
    '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
    }
}

Return type:

dict

astrobase.lcfit.transits.fivetransitparam_fit_magseries(times, mags, errs, teff, rstar, logg, identifier, fit_savdir, chain_savdir, n_mcmc_steps=1, overwriteexistingsamples=False, burninpercent=0.3, n_transit_durations=5, make_tlsfit_plot=True, exp_time_minutes=30, bandpass='tess', magsarefluxes=True, nworkers=32)[source]

Wrapper to mandelagol_fit_magseries that fits out a line around each transit window in the given light curve, and then fits the entire light curve for (t0, period, a/Rstar, Rp/Rstar, inclination). Fixes e to 0, and uses theoretical quadratic limb darkening coefficients in the bandpass given by the user, as found with the stellar parameters. Figures out the priors for you.

Typical use case: you have a light curve with >=2 transits in it. You want to fit the entire light curve for the parameters noted above, but you don’t want to need to manually determine all the priors.

Parameters:
  • times,mags,errs (np.array) – The input flux time-series to fit.
  • teff,rstar,logg (float) – Stellar parameters [K, Rsun, cgs] used to get limb darkening coefficients.
  • identifier (str) – String that goes into file names to identify the object being fit. E.g., fit CSV file will be at {fit_savdir}/{identifier}_fivetransitparam_fitresults.csv
  • fit_savdir (str) – Path to directory where CSV results of fits, fit status files, and diagnostic plots are saved. If it doesn’t exist, this function tries to make it.
  • chain_savdir (str) – Path to directory where MCMC chains are saved.
  • n_mcmc_steps (int) – Number of steps to run MCMC. (Note: convergence not guaranteed).
  • overwriteexistingsamples (bool) – If False, and finds pickle file with saved parameters (in fit_savdir), no additoinal MCMC sampling is done.
  • exp_time_minutes (int) – Exposure time in minutes. Used for the model fitting.
  • n_transit_durations (int) – The points used in the fit are only those within +/- N transit durations of each transit mid-point. This is to prevent excessive out-of-transit data being used in the fit (these points do not inform the model’s parameters).
Returns:

(mafr, tlsr, is_converged)mafr is the Mandel-Agol fit result dictionary, which contains the same information as from mandelagol_and_line_fit_magseries. Fit parameters are accessed like maf_empc_errs['fitinfo']['finalparams']['sma'],

tlsr is the TLS result dictionary, containing keys documented in periodbase/htls.tls_parallel_pfind.

is_converged : boolean for whether the fitting converged, according to the chain autocorrelation time.

Return type:

tuple