<jupyter_start><jupyter_text>Basic FULMAR usage. Using K2-109 as an exampleThis system contains two known exoplanets. For more information, have a look at [Barros et al. 2017](https://www.aanda.org/articles/aa/pdf/2017/12/aa31276-17.pdf)Let's start by importing common modules<jupyter_code>import fulmar
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from transitleastsquares import cleaned_array, transitleastsquares, transit_mask<jupyter_output>WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.<jupyter_text>Functions to plot tls results. To be included in the module in a later version<jupyter_code>def perioplot(tls_results, target, folder, num, maxper=None, savefig=False):
    if maxper is None:
        maxper = np.max(tls_results.periods)
    plt.figure()
    ax = plt.gca()
    ax.axvline(tls_results.period, alpha=0.4, lw=3, color='xkcd:green')
    plt.xlim(np.min(tls_results.periods), maxper)
    for n in range(2, 20):
        ax.axvline(n * tls_results.period, alpha=0.4, lw=1,
                   linestyle="dashed", color='xkcd:green')
    for n in range(2, int(tls_results.period / min(tls_results.periods))):
        ax.axvline(tls_results.period / n, alpha=0.4, lw=1,
                   linestyle="dashed", color='xkcd:green')
    plt.ylabel(r'SDE')
    plt.title(target)
    plt.xlabel('Period (days)')
    plt.annotate("best period = {0:.5f} days".format(tls_results.period), (.02, .02),
            xycoords="axes fraction",
            fontsize=12)
    ax.plot(tls_results.periods, tls_results.power, color='black', lw=0.5)
    plt.xlim(0, maxper)
    plt.tight_layout()
    if savefig is True:
        plt.savefig(
            folder + 'tls_periodogram_{n}_{Pmax}_days'.format(n=num,Pmax=int(maxper)), facecolor='white', dpi=240)
    plt.show()

def modelplot(tls_results, xlim=(0.48, 0.52)):
    plt.figure()
    plt.plot(tls_results.model_folded_phase,
             tls_results.model_folded_model, color='xkcd:green')
    plt.scatter(tls_results.folded_phase, tls_results.folded_y,
                color='xkcd:charcoal', s=10, alpha=0.5, zorder=2)
    plt.xlim(xlim[0], xlim[1])
    plt.xlabel('Phase')
    plt.ylabel('Relative flux')
    plt.show()<jupyter_output><empty_output><jupyter_text>Define target<jupyter_code>lc_targ = fulmar.target('K2-109') # aka. HD 106315 
lc_targ.R_star
lc_targ.mission<jupyter_output><empty_output><jupyter_text>Update stellar parameters if necessary<jupyter_code># K2-109 parameters from Barros et al. 2017 (https://www.aanda.org/articles/aa/pdf/2017/12/aa31276-17.pdf)
lc_targ.R_star = 1.296
lc_targ.R_star_max = 0.058
lc_targ.R_star_min = 0.058

lc_targ.M_star = 1.091
lc_targ.M_star_max = 0.036
lc_targ.M_star_min = 0.036<jupyter_output><empty_output><jupyter_text>Build the lightcurve<jupyter_code>lc_targ.build_lightcurve(author='EVEREST')<jupyter_output>Warning: 25% (851/3384) of the cadences will be ignored due to the quality mask (quality_bitmask=1130799).
Warning: 25% (851/3384) of the cadences will be ignored due to the quality mask (quality_bitmask=1130799).<jupyter_text>Alternatively, build the lightcurve from files, for example data reduced with POLAR<jupyter_code>filelist = [your_files_here]
lc_targ.build_lightcurve(filelist=filelist, author='AUTHOR')<jupyter_output><empty_output><jupyter_text>Plot the light curve for a quick visualisation<jupyter_code>plt.plot(lc_targ.ts_stitch.time.value, lc_targ.ts_stitch['flux'],'k.', markersize=1.8, alpha = 0.25)

plt.xlabel('Time [{}]'.format(str(lc_targ.ts_stitch.time.format).upper()))
plt.ylabel('Flux')
plt.title(lc_targ.K2 + ' Lightcurve')
# plt.xlim(2350,2356)
# plt.ylim(0.9985,1.0015)
plt.show()<jupyter_output><empty_output><jupyter_text>Mask Outliers<jupyter_code>m1 = lc_targ.mask_outliers(sigma=4)<jupyter_output><empty_output><jupyter_text>Plot the light curve with outliers highlighted<jupyter_code>plt.plot(lc_targ.ts_stitch.time.value, lc_targ.ts_stitch['flux'],'k.', markersize=1.8, alpha = 0.25)
plt.plot(lc_targ.ts_stitch.time.value[~m1], lc_targ.ts_stitch['flux'][~m1],'o', color='xkcd:orange', markersize=1.8, alpha = 0.5)

plt.xlabel('Time [{}]'.format(str(lc_targ.ts_stitch.time.format).upper()))
plt.ylabel('Flux')
# plt.xlim(2350,2356)
# plt.ylim(0.985,1.015)
plt.show()<jupyter_output><empty_output><jupyter_text>Remove the outliers if necessary.<jupyter_code>lc_targ.ts_stitch = lc_targ.ts_stitch[m1]<jupyter_output><empty_output><jupyter_text>Clean the light curve / correct activity  Using a Savistzky-Golay filter:<jupyter_code>lc_targ.clean_subt_activity_flatten(sigma=3,
        wl=37,
        time_window=18*u.h,
        polyorder=2,
        return_trend=False,
        remove_outliers=True,
        break_tolerance=5,
        niters=3,
        mask=None)

# You might face issues with Astropy >= 5.0 where flatten gives TypeError: "cannot write to unmasked output"
# The issue comes from Lightkurve and a pull request has been made: https://github.com/lightkurve/lightkurve/pull/1162/commits/05034f4bbd8e6a40f8efd9cfa1f108a36ddc0872<jupyter_output><empty_output><jupyter_text>If you encounter "TypeError: "cannot write to unmasked output"" The temporary fix is to correct this way:<jupyter_code>m_fix=np.full_like(lc_targ.ts_stitch.time.value, False, dtype=bool)
m_fix[0]=True # make the first point invalid

lc_targ.clean_subt_activity_flatten(sigma=3,
        wl=701,
        polyorder=2,
        return_trend=False,
        remove_outliers=True,
        break_tolerance=5,
        niters=3,
        mask=m_fix)

# You might face issues with Astropy >= 5.0 where flatten gives TypeError: "cannot write to unmasked output"
# The issue comes from Lightkurve and a pull request has been made: https://github.com/lightkurve/lightkurve/pull/1162/commits/05034f4bbd8e6a40f8efd9cfa1f108a36ddc0872<jupyter_output><empty_output><jupyter_text>Alternatively, using GP:<jupyter_code>lc_targ.clean_subt_activity_GP(
            timeseries=None,
            bin_duration=60 * u.min,
            period_min=0.2,
            period_max=100.0,
            tune=2500,
            draws=2500,
            chains=2,
            target_accept=0.95,
            ncores=None,
            return_trend=False,
            remove_outliers=True,
            sigma_out=3,
            mask=None)<jupyter_output>/home/jrodrigues/fulmar_venv/lib64/python3.9/site-packages/astropy/timeseries/downsample.py:28: RuntimeWarning: Mean of empty slice /home/jrodrigues/fulmar_venv/lib64/python3.9/site-packages/astropy/timeseries/downsample.py:26: RuntimeWarning: Mean of empty slice /home/jrodrigues/fulmar_venv/lib64/python3.9/site-packages/astropy/timeseries/downsample.py:29: RuntimeWarning: Mean of empty slice<jupyter_text>Plot the corrected light curve<jupyter_code>plt.plot(lc_targ.ts_clean.time.value, lc_targ.ts_clean['flux'],'k.', markersize=1.8, alpha = 0.25)

plt.xlabel('Time [{}]'.format(str(lc_targ.ts_stitch.time.format).upper()))
plt.ylabel('Flux')
# plt.xlim(2350,2356)
# plt.ylim(0.985,1.015)
plt.show()<jupyter_output><empty_output><jupyter_text>Looks like there's something wrong at the begginning of our data. For now, lets mask it out.<jupyter_code>debut_mask = np.array(lc_targ.ts_clean.time.value < 2752.5)
# np.where(debut_mask)<jupyter_output><empty_output><jupyter_text>Look for transits using TLS  Look for the first exoplanet<jupyter_code>perio_results1 = lc_targ.tls_periodogram(cleaned=True, n_transits_min=2, mask=debut_mask)<jupyter_output>Transit Least Squares TLS 1.0.31 (22 Nov 2021)
Creating model cache for 44 durations
Searching 2342 data points, 6881 periods from 0.602 to 33.099 days
Using all 8 CPU threads<jupyter_text>Plot the TLS periodogram<jupyter_code>perioplot(perio_results1, lc_targ.K2, folder="K2-109/", num=1, maxper=None, savefig=False)
# perioplot(results, lc_targ.K2, folder="K2-109/", num=1, maxper=50, savefig=False)<jupyter_output><empty_output><jupyter_text>Plot the resulting model to have a quick visual check<jupyter_code>modelplot(lc_targ.tls_results)<jupyter_output><empty_output><jupyter_text>Optimize the transit parameters. **_Be careful it assumes circular orbit_**<jupyter_code>from fulmar import time_flux_err
import multiprocessing
import exoplanet as xo
def params_optimizer(timeseries, period_guess, t0_guess, depth_guess, ab, r_star, target_id, tran_window=0.25, ncores=None, mask=None):
    if ncores is None:
        ncores = multiprocessing.cpu_count()
    print('running on {} cores'.format(ncores))
#     x = ts_stitch.time.value
#     y = ts_stitch[flux_kw + '_clean'].value
#     yerr = ts_stitch[flux_err_kw+'_clean'].value
        
    if mask is not None:
        x, y, yerr = time_flux_err(timeseries[mask])
        
    x, y, yerr = time_flux_err(timeseries)

    # x = time.copy()
    # y = flux.copy()
    # yerr = flux_err.copy()

    transitMask = (np.abs(
        (x - t0_guess + 0.5 * period_guess) % period_guess - 0.5 * period_guess) < tran_window)
    x = np.ascontiguousarray(x[transitMask])
    y = np.ascontiguousarray(y[transitMask]) - 1
    yerr = np.ascontiguousarray(yerr[transitMask])


#     plt.figure(figsize=(8, 4))
#     x_fold = (
#         x - t0_guess + 0.5 * period_guess
#     ) % period_guess - 0.5 * period_guess
#     plt.scatter(x_fold, y, c=x, s=3)
#     plt.xlabel("time since transit [days]")
#     plt.ylabel("relative flux [ppt]")
#     plt.colorbar(label="time [days]")
#     _ = plt.xlim(-tran_window, tran_window)

    import pymc3 as pm
    import aesara_theano_fallback.tensor as tt

    import pymc3_ext as pmx
#     from celerite2.theano import terms, GaussianProcess

    with pm.Model() as model:

        # Stellar parameters
        mean = pm.Normal("mean", mu=0.0, sigma=10.0)
#         u = xo.distributions.QuadLimbDark("u", testval=np.array(ab))
#         star_params = [mean, u]
        u = ab
        star_params = [mean]

        # Planet parameters
        log_ror = pm.Normal(
            "log_ror", mu=0.5 * np.log(depth_guess), sigma=10.0
        )
        ror = pm.Deterministic("ror", tt.exp(log_ror))
        r_pl = pm.Deterministic("r_pl", ror * r_star)
        # Orbital parameters
        log_period = pm.Normal(
            "log_period", mu=np.log(period_guess), sigma=1.0)
        period = pm.Deterministic("period", tt.exp(log_period))
        t0 = pm.Normal("t0", mu=t0_guess, sigma=1.0)
        log_dur = pm.Normal("log_dur", mu=np.log(0.06), sigma=10.0)
        dur = pm.Deterministic("dur", tt.exp(log_dur))
        b = xo.distributions.ImpactParameter("b", ror=ror)

        # Set up the orbit
        orbit = xo.orbits.KeplerianOrbit(
            period=period, duration=dur, ror=ror, t0=t0, b=b)

        # We're going to track the implied density
        pm.Deterministic("rho_circ", orbit.rho_star)

        # Set up the mean transit model
        light_curves = xo.LimbDarkLightCurve(
            u).get_light_curve(orbit=orbit, r=ror, t=x)

        light_curve = pm.math.sum(light_curves, axis=-1) + mean

        # Here we track the value of the model light curve for plotting
        # purposes
        pm.Deterministic("light_curves", light_curves)

        # Finally the GP observation model
    #     gp = GaussianProcess(
    #         kernel, t=x, diag=yerr ** 2 + sigma ** 2, mean=lc_model
    #     )
    #     gp.marginal("obs", observed=y)
    #     pm.Deterministic("gp_pred", gp.predict(y))

        pm.Normal("obs", mu=light_curve, sd=np.median(yerr), observed=y)

        # Double check that everything looks good - we shouldn't see any NaNs!
        print(model.check_test_point())

        # Optimize the model
        map_soln = model.test_point

        map_soln = pmx.optimize(map_soln, [ror, b, dur])

        map_soln = pmx.optimize(map_soln, star_params)
        map_soln = pmx.optimize(map_soln)
        map_soln = pmx.optimize()


#         plt.figure(figsize=(9, 5))
#         x_fold = (x - map_soln["t0"] + 0.5 * map_soln["period"]) % map_soln[
#             "period"
#         ] - 0.5 * map_soln["period"]
#         inds = np.argsort(x_fold)
#         plt.scatter(x_fold, 1 + y - map_soln["mean"], c=x, s=3)
#         plt.plot(x_fold[inds], 1 + map_soln["light_curves"][inds] - map_soln["mean"], "k")
#         plt.xlabel("time since transit [days]")
#         plt.ylabel("relative flux [ppt]")
#         plt.colorbar(label="time [days]")
#         _ = plt.xlim(-tran_window, tran_window)
#         plt.show()

        np.random.seed()
        with model:
            trace = pmx.sample(
                tune=2500,
                draws=2000,
                start=map_soln,
                chains=2,
                cores=ncores,
                target_accept=0.96,
                return_inferencedata=True,
            )

        import arviz as az
        az.summary(trace,
                   var_names=[
                       "period",
                       "t0",
                       "ror",
                       'dur',
                       'b',
                       #                     "u",
                       "mean"
                   ],)

        flat_samps = trace.posterior.stack(sample=("chain", "draw"))
        p = np.median(flat_samps["period"])
        t0 = np.median(flat_samps["t0"])
        dur = np.median(flat_samps["dur"])
        depth = np.median(flat_samps['ror'])**2
#         ab = tuple((np.median(flat_samps['u'], axis=-1)))

        # Plot the folded data
        x_fold = (x - t0 + 0.5 * p) % p - 0.5 * p
        plt.plot(x_fold, 1 + y, ".k", alpha=0.4, label="data", zorder=-1000)

        # Overplot the phase binned light curve
        bins = np.linspace(-0.41, 0.41, 50)
        denom, _ = np.histogram(x_fold, bins)
        num, _ = np.histogram(x_fold, bins, weights=y)
        denom[num == 0] = 1.0
        plt.plot(
            0.5 * (bins[1:] + bins[:-1]), 1 + num / denom, "o", color="C1", label="binned", alpha=0.7
        )

        # Plot the folded model
        inds = np.argsort(x_fold)
        inds = inds[np.abs(x_fold)[inds] < 0.3]
        pred = np.percentile(
            flat_samps["light_curves"][inds, 0], [16, 50, 84], axis=-1
        )
        plt.plot(x_fold[inds], 1 + pred[1], color="xkcd:green", label="model")
        art = plt.fill_between(
            x_fold[inds], 1 + pred[0], 1 + pred[2], color="xkcd:green", alpha=0.2, zorder=1000
        )
        art.set_edgecolor("none")

        # Annotate the plot with the planet's period
        txt = "period = {0:.5f} +/- {1:.5f} d".format(
            np.mean(flat_samps["period"].values), np.std(
                flat_samps["period"].values)
        )
        plt.annotate(
            txt,
            (0, 0),
            xycoords="axes fraction",
            xytext=(5, 5),
            textcoords="offset points",
            ha="left",
            va="bottom",
            fontsize=12,
        )

        plt.legend(fontsize=10, loc=4)
        plt.title(target_id)
        plt.xlim(-0.5 * p, 0.5 * p)
        plt.xlabel("time since transit [days]")
        plt.ylabel("de-trended flux")
        _ = plt.xlim(-tran_window, tran_window)
        plt.show()

        return p, t0, dur, depth, ab, flat_samps
p, t0, dur, depth, ab, flat_samps = params_optimizer(lc_targ.ts_clean, lc_targ.tls_results.period, lc_targ.tls_results.T0, 1-lc_targ.tls_results.depth, lc_targ.ab, lc_targ.R_star, lc_targ.K2, tran_window=0.25, ncores=None, mask=None)<jupyter_output>running on 8 cores
mean           -3.22
log_ror        -3.22
log_period     -0.92
t0             -0.92
log_dur        -3.22
b_impact__     -1.39
obs           215.39
Name: Log-probability of test_point, dtype: float64<jupyter_text>Mask intransit data<jupyter_code>intransit = transit_mask(lc_targ.ts_clean.time.value, p, 1.3*dur, t0) # True when datapoints are in transit<jupyter_output><empty_output><jupyter_text>Combine intransit mask with the initial mask<jupyter_code>p1_mask = np.logical_or(debut_mask, intransit)<jupyter_output><empty_output><jupyter_text>Plot the lightcurve with transits highlighted<jupyter_code>plt.plot(lc_targ.ts_clean.time.value, lc_targ.ts_clean['flux'],'k.', markersize=1.8, alpha = 0.25)
plt.plot(lc_targ.ts_clean.time.value[intransit], lc_targ.ts_clean['flux'][intransit],'o', color='xkcd:green', markersize=4.8, alpha = 0.5)

plt.xlabel('Time [{}]'.format(str(lc_targ.ts_stitch.time.format).upper()))
plt.ylabel('Flux')
#plt.xlim(2749,2752)
# plt.ylim(0.985,1.015)
plt.show()<jupyter_output><empty_output><jupyter_text>Look for the second exoplanet<jupyter_code>perio_results2 = lc_targ.tls_periodogram(cleaned=True, n_transits_min=2, period_min=0.9, mask=p1_mask)<jupyter_output>Transit Least Squares TLS 1.0.31 (22 Nov 2021)
Creating model cache for 44 durations
Searching 2318 data points, 5707 periods from 0.9 to 33.099 days
Using all 8 CPU threads<jupyter_text>Plot the TLS periodogram<jupyter_code>perioplot(perio_results2, lc_targ.K2, folder="K2-109/", num=2, maxper=None, savefig=False)
# perioplot(results, lc_targ.K2, folder="K2-109/", num=1, maxper=50, savefig=False)
modelplot(lc_targ.tls_results)<jupyter_output><empty_output><jupyter_text>Visualize where the transits lay in the lightcurve<jupyter_code>plt.plot(lc_targ.ts_clean.time.value[~debut_mask], lc_targ.ts_clean['flux'][~debut_mask],'ko', markersize=1.8, alpha = 0.2)

plt.xlabel('Time [{}]'.format(str(lc_targ.ts_stitch.time.format).upper()))
plt.ylabel('Flux')
#plt.xlim(2749,2752)
# plt.ylim(0.985,1.015)

for n in perio_results1.transit_times:
    plt.axvline(n, alpha=0.4, lw=2, ymax=0.15, color='xkcd:green')
for n in perio_results2.transit_times:
    plt.axvline(n, alpha=0.3, lw=2, ymax=0.45, color='xkcd:orange')

plt.title(lc_targ.K2)

plt.ylim(0.9985,1.00028)
plt.show()<jupyter_output><empty_output><jupyter_text>Optimize the transit parameters. **_Be careful it still assumes circular orbit_**<jupyter_code>p_1, t0_1, dur_1, depth_1, ab_1, flat_samps_1 = params_optimizer(lc_targ.ts_clean, lc_targ.tls_results.period, lc_targ.tls_results.T0, 1-lc_targ.tls_results.depth, lc_targ.ab, lc_targ.R_star, lc_targ.K2, tran_window=0.25, ncores=None, mask=~p1_mask)<jupyter_output>running on 8 cores
mean           -3.22
log_ror        -3.22
log_period     -0.92
t0             -0.92
log_dur        -3.22
b_impact__     -1.39
obs           958.91
Name: Log-probability of test_point, dtype: float64<jupyter_text>Mask intransit data<jupyter_code>intransit2 = transit_mask(lc_targ.ts_clean.time.value, p_1, 1.3*dur_1, t0_1) # True when datapoints are in transit<jupyter_output><empty_output><jupyter_text>Combine intransit mask with the initial mask<jupyter_code>p2_mask = np.logical_or(p1_mask, intransit2)<jupyter_output><empty_output><jupyter_text>Plot the lightcurve with transits highlighted<jupyter_code>plt.plot(lc_targ.ts_clean.time.value[~debut_mask], lc_targ.ts_clean['flux'][~debut_mask],'k.', markersize=1.8, alpha = 0.25)
plt.plot(lc_targ.ts_clean.time.value[intransit2], lc_targ.ts_clean['flux'][intransit2],'o', color='xkcd:orange', markersize=4.8, alpha = 0.25, label='K2-279b')
plt.plot(lc_targ.ts_clean.time.value[intransit], lc_targ.ts_clean['flux'][intransit],'o', color='xkcd:green', markersize=4.8, alpha = 0.25, label='K2-279c')


plt.xlabel('Time [{}]'.format(str(lc_targ.ts_stitch.time.format).upper()))
plt.ylabel('Flux')
#plt.xlim(2749,2752)
# plt.ylim(0.985,1.015)
plt.title(lc_targ.K2)
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Look for a possible third exoplanet ?<jupyter_code>perio_results3 = lc_targ.tls_periodogram(cleaned=True, n_transits_min=2, mask=p2_mask)<jupyter_output>Transit Least Squares TLS 1.0.31 (22 Nov 2021)
Creating model cache for 44 durations
Searching 2265 data points, 6881 periods from 0.602 to 33.099 days
Using all 8 CPU threads<jupyter_text>Plot the TLS periodogram<jupyter_code>perioplot(perio_results3, lc_targ.K2, folder="K2-109/", num=3, maxper=None, savefig=False)
# perioplot(results, lc_targ.K2, folder="K2-109/", num=1, maxper=50, savefig=False)
modelplot(lc_targ.tls_results)
perio_results3.SDE, perio_results3.FAP<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13363640905230553, 'fraction_numerical': 0.04161462756656842, 'mean_word_length': 3.5630271255983588, 'pattern_counts': {'":': 0, '<': 120, '<?xml version=': 0, '>': 119, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '31', 'dir': 'jupyter-structured-clean-dedup', 'id': '96571', 'license': 'MIT', 'n_tokens_mistral': 8221, 'n_tokens_neox': 7585, 'n_words': 1344, 'path': 'Examples/fulmar_example.ipynb', 'repo_name': 'astrojose9/Fulmar', 'size': '20017'}