<jupyter_start><jupyter_text>KIC 9651065<jupyter_code>%run setup.py
t, y = np.loadtxt('../lc/9651065_lc.txt', usecols=(0,1)).T
ms = Maelstrom(t, y, max_peaks=5, fmin=5, fmax=48)
ms.first_look()
period_guess = 300
a_guess = 200

time, flux = ms.time, ms.flux
freq = ms.freq
weights = ms.get_weights(norm=False)
pg = ms.period_search()
periods = np.linspace(100, 300, 300)
results = pg.fit(periods)
ys = np.array([[r[0] for r in row] for row in results])
sm = np.sum(ys, axis=0)
period_ind = np.argmax(sm)
plt.plot(periods[:-2], sm[:-2]);
from maelstrom.utils import unique_colors

hh = unique_colors(len(ms.freq), cmap='Blues')
plt.figure(figsize=mnras_size(240.))

ys = np.array([[np.exp(r[1]["logasini"]) for r in row] for row in results])
for i, c in zip(ys, hh):
    plt.plot(periods, i, alpha=1, linewidth=0.8, c=c);

plt.xlabel('Period (d)')
plt.ylabel(r'a$\sin{i}$ (s)')
plt.ylim(0, None)
plt.xlim(100,300)
plt.axhline(184., c='r', linestyle='dashed', linewidth=0.7)
plt.axvline(272., c='r', linestyle='dashed', linewidth=0.7)

plt.savefig(overleaf_path + '9651065_period_search.pdf', dpi=300, bbox_inches='tight', pad_inches=0)<jupyter_output><empty_output><jupyter_text>Maelstrom model<jupyter_code>ms.setup_orbit_model(period=period_guess)
# opt = ms.optimize()
pb1 = ms.pin_orbit_model()
opt = pb1.optimize()
opt
# with pb1:
#     trace = pm.load_trace('traces/9651065_FINAL_VERSION2/')
with pb1:
    trace = pm.sample(
        tune=1000,
        draws=2000,
        start=opt,
        chains=2,
        step=xo.get_dense_nuts_step(target_accept=0.9),
    )
pm.save_trace(trace, 'trace/NEW/9651065')
with pb1:
    trace = pm.load_trace('trace/NEW/9651065')
pm.summary(trace)
from tqdm import tqdm

taus = []
with pb1:
    for samp in tqdm(xo.utils.get_samples_from_trace(trace, size=1000), total=1000):
            taus.append(xo.eval_in_model(pb1.orbit.get_time_delay(time), samp) * 86400)
med_td = np.median(taus, axis=0)
sd_td = np.std(taus, axis=0)
mean = np.mean(taus)
mean
np.random.seed(23)
fig, ax = plt.subplots(figsize=mnras_size(240), constrained_layout=True)

ax.set_rasterized(True)
#ax.set_rasterization_zorder(1)

with pb1:
    for samp in xo.utils.get_samples_from_trace(trace, size=10):
        taumod = xo.eval_in_model(pb1.orbit.get_time_delay(time), samp) * 86400
        #ttime = (ms.time_mid + time - samp['tref']) % samp['period'] / samp['period']
        ttime = (time) % samp['PB1_period'] / samp['PB1_period']
        #ttime = ((ms.time_mid + time) + (samp['phi'] * samp['period'] / (2*np.pi))) % samp['period'] / samp['period']
        sort = np.argsort(ttime)
        ax.plot(ttime[sort], (taumod - np.mean(taumod))[sort], color=blue, linewidth=0.4, alpha=1, 
#                 rasterized=True,
                zorder=1)


ax.set_xlabel('Orbital phase')
ax.set_ylabel('Time delay (s)', c=blue)


times = time# + xo.eval_in_model(phi * period / (2*np.pi), samp)
fold = times % np.median(trace['PB1_period']) / np.median(trace['PB1_period'])
sort = np.argsort(fold)
plt.fill_between(fold[sort], (med_td - sd_td - 108.12878089572754)[:,0][sort], (med_td+sd_td - 108.12878089572754)[:,0][sort], alpha=0.2, color=blue)


ax.set_xlim(0, 1)

plt.savefig(overleaf_path + '9651065.pdf', dpi=300, bbox_inches='tight', pad_inches=0)
#plt.savefig('rast.pdf', dpi=300, bbox_inches='tight')
from maelstrom.utils import mass_function
import astropy.units as u
rounding = 3
samples = pm.trace_to_dataframe(trace, varnames=['PB1_period', 'PB1_asini'])
mfs = mass_function(samples['PB1_period'].values * u.day, samples['PB1_asini'].values*u.s)
#mfs = np.array(mfs)
upper, med, lower = np.percentile(mfs.value, [84.13, 50, 15.86])
print('mass_func', ': ', np.round(med,rounding), ' + ', np.round(upper - med,rounding), ' - ', np.round(med - lower,rounding))
varnames = ["period", "asini", "eccen", "omega", "phi"]
for var in varnames:
    percentiles = np.percentile(trace['PB1_' + var], q=[15.87, 50, 84.13])
    print(f'{var}: {percentiles[1]:.2f} + {percentiles[1] - percentiles[0]:.2f} - {percentiles[2] - percentiles[1]:.2f}')<jupyter_output>period: 272.27 + 0.42 - 0.42
asini: 185.22 + 3.48 - 3.70
eccen: 0.45 + 0.03 - 0.03
omega: 2.15 + 0.07 - 0.07
phi: 0.52 + 0.07 - 0.07<jupyter_text>Subdividing model<jupyter_code>td_time, td_td, td_err = np.loadtxt('../data/kic9651065_uncertainties-plus-time-delay_Q99_llc.txt', delimiter=',', usecols=(0,1,2)).T
td_time += 2400000
td_time -= 2454833
plt.scatter(td_time, td_td)
import theano.tensor as tt
from maelstrom.orbit import Orbit

with pm.Model() as subdivide_model:
    logP = pm.Normal("logP", mu=np.log(272), sd=1.0, testval=np.log(272))
    period = pm.Deterministic("period", pm.math.exp(logP))
    
    # The time of conjunction
    phi = xo.distributions.Angle("phi", testval=0.5691498)
    logs_lc = pm.Normal('logs_lc', mu=np.log(np.std(flux)), sd=10, testval=0.)
    logasini = pm.Normal('logasini', mu=np.log(184), sd=10, testval=np.log(184))
    asini = pm.Normal("asini", mu=184, sd=10, testval=184)
    drift = pm.Normal('drift', mu=0., sd=0.1, testval=0)
    # Periastron sampled from uniform angle
    omega = xo.distributions.Angle("omega", testval=-0.94)
    # Eccentricity
    eccen = pm.Uniform("eccen", lower=0, upper=0.9, testval=0.45)
    # The baseline flux
    mean = pm.Normal("mean", mu=0.0, sd=10.0, testval=0.003)
    
    # Here, we generate an Orbit instance and pass in our priors. 
    orbit = Orbit(period=period, 
                  lighttime=asini, 
                  omega=omega, 
                  eccen=eccen, 
                  phi=phi, 
                  freq=0)
    
    # psi is defined to be negative but the light curve model takes 2*pi*f * (time - tau), so 
    # we must flip tau here to phase it on the same values
    td = -1*tt.squeeze(orbit.get_time_delay(td_time) * 86400) # Convert to s
    td += drift * td_time
    taumodel = pm.Deterministic('taumodel', td - tt.mean(td))
    
    pm.Normal('obs', mu=taumodel, sd=tt.exp(logs_lc), observed=td_td)
    
    plt.plot(td_time, xo.eval_in_model(taumodel))
    plt.plot(td_time, td_td)
with subdivide_model:
    opt = xo.optimize()
opt
with subdivide_model:
    trace = pm.sample(draws=2000, tune=2000, chains=2, start=opt)
pm.summary(trace)
varnames=['period', 'phi', 'eccen', 'asini', 'omega', 'phi', 'drift']
rounding = 2
for varname in varnames:
    upper, med, lower = np.percentile(trace[varname], [84.13, 50, 15.86])
    print(varname, ': ', np.round(med,rounding), ' + ', np.round(upper - med,rounding), ' - ', np.round(med - lower,rounding))
from maelstrom.utils import mass_function
import astropy.units as u
rounding = 3
samples = pm.trace_to_dataframe(trace, varnames=['period', 'asini'])
mfs = mass_function(samples['period'].values * u.day, samples['asini'].values*u.s)
#mfs = np.array(mfs)
upper, med, lower = np.percentile(mfs.value, [84.13, 50, 15.86])
print('mass_func', ': ', np.round(med,rounding), ' + ', np.round(upper - med,rounding), ' - ', np.round(med - lower,rounding))
plt.scatter(td_time % 272 / 272, np.median(trace['taumodel'], axis=0))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1762126997595814, 'fraction_numerical': 0.05939753924480271, 'mean_word_length': 4.177159590043924, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '23081', 'license': 'MIT', 'n_tokens_mistral': 3157, 'n_tokens_neox': 2874, 'n_words': 534, 'path': 'paper/notebooks/(4.1) KIC 9651065.ipynb', 'repo_name': 'danielhey/Maelstrom', 'size': '7071'}