def mandelagol_and_line_fit_magseries(<EOL>times, mags, errs,<EOL>fitparams,<EOL>priorbounds,<EOL>fixedparams,<EOL>trueparams=None,<EOL>burninpercent=<NUM_LIT>,<EOL>plotcorner=False,<EOL>timeoffset=<NUM_LIT:0>,<EOL>samplesavpath=False,<EOL>n_walkers=<NUM_LIT:50>,<EOL>n_mcmc_steps=<NUM_LIT>,<EOL>eps=<NUM_LIT>,<EOL>skipsampling=False,<EOL>overwriteexistingsamples=False,<EOL>mcmcprogressbar=False,<EOL>plotfit=False,<EOL>scatterxdata=None,<EOL>scatteryaxes=None,<EOL>magsarefluxes=True,<EOL>sigclip=<NUM_LIT>,<EOL>verbose=True,<EOL>nworkers=<NUM_LIT:4><EOL>):
from multiprocessing import Pool<EOL>fittype = '<STR_LIT>'<EOL>if not magsarefluxes:<EOL><INDENT>raise NotImplementedError('<STR_LIT>')<EOL><DEDENT>if not samplesavpath:<EOL><INDENT>raise ValueError(<EOL>'<STR_LIT>'<EOL>)<EOL><DEDENT>if not mandel_agol_dependencies:<EOL><INDENT>raise ImportError(<EOL>'<STR_LIT>'<EOL>)<EOL><DEDENT>stimes, smags, serrs = sigclip_magseries(times, mags, errs,<EOL>sigclip=sigclip,<EOL>magsarefluxes=magsarefluxes)<EOL>nzind = np.nonzero(serrs)<EOL>stimes, smags, serrs = stimes[nzind], smags[nzind], serrs[nzind]<EOL>init_period = _get_value('<STR_LIT>', fitparams, fixedparams)<EOL>init_epoch = _get_value('<STR_LIT>', fitparams, fixedparams)<EOL>init_rp = _get_value('<STR_LIT>', fitparams, fixedparams)<EOL>init_sma = _get_value('<STR_LIT>', fitparams, fixedparams)<EOL>init_incl = _get_value('<STR_LIT>', fitparams, fixedparams)<EOL>init_ecc = _get_value('<STR_LIT>', fitparams, fixedparams)<EOL>init_omega = _get_value('<STR_LIT>', fitparams, fixedparams)<EOL>limb_dark = _get_value('<STR_LIT>', fitparams, fixedparams)<EOL>init_u = _get_value('<STR_LIT:u>', fitparams, fixedparams)<EOL>init_poly_order0 = _get_value('<STR_LIT>', fitparams, fixedparams)<EOL>init_poly_order1 = _get_value('<STR_LIT>', fitparams, fixedparams)<EOL>if not limb_dark == '<STR_LIT>':<EOL><INDENT>raise ValueError(<EOL>'<STR_LIT>'<EOL>)<EOL><DEDENT>init_params, init_m = _transit_model(<EOL>stimes, init_epoch, init_period, init_rp, init_sma, init_incl,<EOL>init_ecc, init_omega, init_u, limb_dark)<EOL>init_flux = (<EOL>init_m.light_curve(init_params) +<EOL>init_poly_order0 + init_poly_order1*stimes<EOL>)<EOL>theta, fitparamnames = [], []<EOL>for k in np.sort(list(fitparams.keys())):<EOL><INDENT>if isinstance(fitparams[k], float) or isinstance(fitparams[k], int):<EOL><INDENT>theta.append(fitparams[k])<EOL>fitparamnames.append(fitparams[k])<EOL><DEDENT>elif isinstance(fitparams[k], list):<EOL><INDENT>if not len(fitparams[k]) == <NUM_LIT:2>:<EOL><INDENT>raise ValueError('<STR_LIT>')<EOL><DEDENT>theta.append(fitparams[k][<NUM_LIT:0>])<EOL>theta.append(fitparams[k][<NUM_LIT:1>])<EOL>fitparamnames.append(fitparams[k][<NUM_LIT:0>])<EOL>fitparamnames.append(fitparams[k][<NUM_LIT:1>])<EOL><DEDENT><DEDENT>n_dim = len(theta)<EOL>if not skipsampling:<EOL><INDENT>backend = emcee.backends.HDFBackend(samplesavpath)<EOL>if overwriteexistingsamples:<EOL><INDENT>LOGWARNING(<EOL>'<STR_LIT>'.format(samplesavpath)<EOL>)<EOL>backend.reset(n_walkers, n_dim)<EOL><DEDENT>def nll(*args):<EOL><INDENT>return -_log_likelihood_transit_plus_line(*args)<EOL><DEDENT>soln = spminimize(<EOL>nll, theta, method='<STR_LIT>',<EOL>args=(init_params, init_m, stimes, smags, serrs, priorbounds)<EOL>)<EOL>theta_ml = soln.x<EOL>ml_poly_order0 = theta_ml[<NUM_LIT:0>]<EOL>ml_poly_order1 = theta_ml[<NUM_LIT:1>]<EOL>ml_rp = theta_ml[<NUM_LIT:2>]<EOL>ml_t0 = theta_ml[<NUM_LIT:3>]<EOL>ml_params, ml_m = _transit_model(stimes, ml_t0, init_period,<EOL>ml_rp, init_sma, init_incl,<EOL>init_ecc, init_omega, init_u,<EOL>limb_dark)<EOL>ml_mags = (<EOL>ml_m.light_curve(ml_params) +<EOL>ml_poly_order0 + ml_poly_order1*stimes<EOL>)<EOL>initial_position_vec = [theta_ml + eps*np.random.randn(n_dim)<EOL>for i in range(n_walkers)]<EOL>starting_positions = initial_position_vec<EOL>isfirstrun = True<EOL>if os.path.exists(backend.filename):<EOL><INDENT>if backend.iteration > <NUM_LIT:1>:<EOL><INDENT>starting_positions = None<EOL>isfirstrun = False<EOL><DEDENT><DEDENT>if verbose and isfirstrun:<EOL><INDENT>LOGINFO(<EOL>'<STR_LIT>'.<EOL>format(fittype, n_dim, n_mcmc_steps, n_walkers) +<EOL>'<STR_LIT>'.format(nworkers)<EOL>)<EOL><DEDENT>elif verbose and not isfirstrun:<EOL><INDENT>LOGINFO(<EOL>'<STR_LIT>'.<EOL>format(fittype, n_dim, n_mcmc_steps, n_walkers) +<EOL>'<STR_LIT>'.format(nworkers)<EOL>)<EOL><DEDENT>with Pool(nworkers) as pool:<EOL><INDENT>sampler = emcee.EnsembleSampler(<EOL>n_walkers, n_dim, log_posterior_transit_plus_line,<EOL>args=(init_params, init_m, stimes, smags, serrs, priorbounds),<EOL>pool=pool,<EOL>backend=backend<EOL>)<EOL>sampler.run_mcmc(starting_positions, n_mcmc_steps,<EOL>progress=mcmcprogressbar)<EOL><DEDENT>if verbose:<EOL><INDENT>LOGINFO(<EOL>'<STR_LIT>'.format(<EOL>fittype, n_mcmc_steps, n_walkers<EOL>) + '<STR_LIT>'.format(nworkers)<EOL>)<EOL><DEDENT><DEDENT>reader = emcee.backends.HDFBackend(samplesavpath)<EOL>n_to_discard = int(burninpercent*n_mcmc_steps)<EOL>samples = reader.get_chain(discard=n_to_discard, flat=True)<EOL>log_prob_samples = reader.get_log_prob(discard=n_to_discard, flat=True)<EOL>log_prior_samples = reader.get_blobs(discard=n_to_discard, flat=True)<EOL>fit_statistics = list(<EOL>map(lambda v: (v[<NUM_LIT:1>], v[<NUM_LIT:2>]-v[<NUM_LIT:1>], v[<NUM_LIT:1>]-v[<NUM_LIT:0>]),<EOL>list(zip( *np.percentile(samples, [<NUM_LIT>, <NUM_LIT:50>, <NUM_LIT>], axis=<NUM_LIT:0>))))<EOL>)<EOL>medianparams, std_perrs, std_merrs = {}, {}, {}<EOL>for ix, k in enumerate(np.sort(list(priorbounds.keys()))):<EOL><INDENT>medianparams[k] = fit_statistics[ix][<NUM_LIT:0>]<EOL>std_perrs[k] = fit_statistics[ix][<NUM_LIT:1>]<EOL>std_merrs[k] = fit_statistics[ix][<NUM_LIT:2>]<EOL><DEDENT>stderrs = {'<STR_LIT>':std_perrs, '<STR_LIT>':std_merrs}<EOL>per = _get_value('<STR_LIT>', medianparams, fixedparams)<EOL>t0 = _get_value('<STR_LIT>', medianparams, fixedparams)<EOL>rp = _get_value('<STR_LIT>', medianparams, fixedparams)<EOL>sma = _get_value('<STR_LIT>', medianparams, fixedparams)<EOL>incl = _get_value('<STR_LIT>', medianparams, fixedparams)<EOL>ecc = _get_value('<STR_LIT>', medianparams, fixedparams)<EOL>omega = _get_value('<STR_LIT>', medianparams, fixedparams)<EOL>limb_dark = _get_value('<STR_LIT>', medianparams, fixedparams)<EOL>try:<EOL><INDENT>u = fixedparams['<STR_LIT:u>']<EOL><DEDENT>except Exception as e:<EOL><INDENT>u = [medianparams['<STR_LIT>'], medianparams['<STR_LIT>']]<EOL><DEDENT>poly_order0 = _get_value('<STR_LIT>', medianparams, fixedparams)<EOL>poly_order1 = _get_value('<STR_LIT>', medianparams, fixedparams)<EOL>fit_params, fit_m = _transit_model(stimes, t0, per, rp, sma, incl, ecc,<EOL>omega, u, limb_dark)<EOL>fitmags = (<EOL>fit_m.light_curve(fit_params) +<EOL>poly_order0 + poly_order1*stimes<EOL>)<EOL>fepoch = t0<EOL>medianparams['<STR_LIT>'] += timeoffset<EOL>returndict = {<EOL>'<STR_LIT>':fittype,<EOL>'<STR_LIT>':{<EOL>'<STR_LIT>':fitparams,<EOL>'<STR_LIT>':init_flux,<EOL>'<STR_LIT>':fixedparams,<EOL>'<STR_LIT>':medianparams,<EOL>'<STR_LIT>':stderrs,<EOL>'<STR_LIT>':fitmags,<EOL>'<STR_LIT>':fepoch+timeoffset,<EOL>},<EOL>'<STR_LIT>':None,<EOL>'<STR_LIT>':{<EOL>'<STR_LIT>':stimes+timeoffset,<EOL>'<STR_LIT>':smags,<EOL>'<STR_LIT>':serrs,<EOL>'<STR_LIT>':magsarefluxes,<EOL>},<EOL>}<EOL>if plotcorner:<EOL><INDENT>fig = corner.corner(<EOL>samples,<EOL>labels=['<STR_LIT>', '<STR_LIT>',<EOL>'<STR_LIT>','<STR_LIT>'.format(timeoffset)],<EOL>truths=[ml_poly_order0, ml_poly_order1, ml_rp, ml_t0],<EOL>quantiles=[<NUM_LIT>, <NUM_LIT:0.5>, <NUM_LIT>], show_titles=True<EOL>)<EOL>plt.savefig(plotcorner, dpi=<NUM_LIT>)<EOL>if verbose:<EOL><INDENT>LOGINFO('<STR_LIT>'.format(plotcorner))<EOL><DEDENT><DEDENT>if plotfit and isinstance(plotfit, str):<EOL><INDENT>plt.close('<STR_LIT:all>')<EOL>f, (a0, a1) = plt.subplots(nrows=<NUM_LIT:2>, ncols=<NUM_LIT:1>, sharex=True,<EOL>figsize=(<NUM_LIT:8>,<NUM_LIT:5>),<EOL>gridspec_kw={'<STR_LIT>':[<NUM_LIT:3>, <NUM_LIT:1>]})<EOL>a0.scatter(stimes, smags, c='<STR_LIT:k>', alpha=<NUM_LIT>, label='<STR_LIT:data>', zorder=<NUM_LIT:1>,<EOL>s=<NUM_LIT:10>, rasterized=True, linewidths=<NUM_LIT:0>)<EOL>DEBUGGING = False<EOL>if DEBUGGING:<EOL><INDENT>a0.scatter(stimes, init_flux, c='<STR_LIT:r>', alpha=<NUM_LIT:1>, s=<NUM_LIT>, zorder=<NUM_LIT:2>,<EOL>rasterized=True, linewidths=<NUM_LIT:0>,<EOL>label='<STR_LIT>')<EOL>a0.scatter(stimes, ml_mags, c='<STR_LIT:g>', alpha=<NUM_LIT:1>, s=<NUM_LIT>, zorder=<NUM_LIT:2>,<EOL>rasterized=True, linewidths=<NUM_LIT:0>, label='<STR_LIT>')<EOL><DEDENT>a0.plot(<EOL>stimes, fitmags, c='<STR_LIT:b>',<EOL>zorder=<NUM_LIT:0>, rasterized=True, lw=<NUM_LIT:2>, alpha=<NUM_LIT>,<EOL>label='<STR_LIT>'.format(fittype, len(fitparamnames))<EOL>)<EOL>a1.scatter(<EOL>stimes, smags-fitmags, c='<STR_LIT:k>', alpha=<NUM_LIT>,<EOL>rasterized=True, s=<NUM_LIT:10>, linewidths=<NUM_LIT:0><EOL>)<EOL>if scatterxdata and scatteryaxes:<EOL><INDENT>import matplotlib.transforms as transforms<EOL>for a in [a0, a1]:<EOL><INDENT>transform = transforms.blended_transform_factory(<EOL>a.transData, a.transAxes<EOL>)<EOL>a.scatter(scatterxdata, scatteryaxes, c='<STR_LIT:r>', alpha=<NUM_LIT>,<EOL>zorder=<NUM_LIT:2>, s=<NUM_LIT:10>, rasterized=True, linewidths=<NUM_LIT:0>,<EOL>marker="<STR_LIT>", transform=transform)<EOL><DEDENT><DEDENT>a1.set_xlabel('<STR_LIT>')<EOL>a0.set_ylabel('<STR_LIT>')<EOL>a1.set_ylabel('<STR_LIT>')<EOL>a0.legend(loc='<STR_LIT>', fontsize='<STR_LIT>')<EOL>for a in [a0, a1]:<EOL><INDENT>a.get_yaxis().set_tick_params(which='<STR_LIT>', direction='<STR_LIT>')<EOL>a.get_xaxis().set_tick_params(which='<STR_LIT>', direction='<STR_LIT>')<EOL><DEDENT>f.tight_layout(h_pad=<NUM_LIT:0>, w_pad=<NUM_LIT:0>)<EOL>f.savefig(plotfit, dpi=<NUM_LIT>, bbox_inches='<STR_LIT>')<EOL>if verbose:<EOL><INDENT>LOGINFO('<STR_LIT>'.format(plotfit))<EOL><DEDENT>returndict['<STR_LIT>'] = plotfit<EOL><DEDENT>return returndict<EOL>
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
                }
            }