<jupyter_start><jupyter_text>Let's measure the time delay curve using a linear spline. First, load the data:<jupyter_code>kicid = 11754974
data = np.loadtxt("data/kic{0}_lc.txt".format(kicid))
fulltimes = data[:, 0] # days
tmid = 0.5*(fulltimes[0] + fulltimes[-1])
times = fulltimes - tmid
dmmags = data[:, 1] * 1000. # mmags

metadata = np.loadtxt("data/kic{0}_metadata_more.csv".format(kicid), delimiter=",", skiprows=1)
nu_arr = metadata[::6]

plt.plot(times, dmmags, "k")
plt.xlabel(r"$t$")
plt.ylabel(r"$L(t)$");
nnu = 8
nu_arr<jupyter_output><empty_output><jupyter_text>Now, let's set up the model. We'll use 1000 control points in the spline and we'll start with them randomly distributed with a standard deviation of 1e-5.<jupyter_code>N = 1000
model = InterpMaelstrom(times, dmmags, nu_arr[:nnu], log_sigma2=0.0,
                        interp_x=np.linspace(times.min(), times.max(), N),
                        interp_y=1e-5*np.random.randn(N))<jupyter_output><empty_output><jupyter_text>Here's the initial time delay curve:<jupyter_code>plt.plot(times, model.run(model.tau)[:, 0], ".k", ms=1)
plt.ylabel(r"$\tau(t)$")
plt.xlabel("$t$");<jupyter_output><empty_output><jupyter_text>Now let's optimize the frequency vector and the noise variance:<jupyter_code>model.optimize([model.log_sigma2, model.nu])<jupyter_output>INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 259485.844005
  Number of iterations: 109
  Number of functions evaluations: 142<jupyter_text>Then we'll fit for the time delay curve at a set of 1000 control points as defined above. We're using a gradient optimizer that tends to be better if you have a lot of parameters (it was designed for neural networks).<jupyter_code>opt = tf.train.AdamOptimizer().minimize(
    model.chi2, var_list=[model.interp_y, model.log_sigma2])
model.run(tf.global_variables_initializer())

bar = tqdm.trange(1000)
for i in bar:
    chi2, _ = model.run([model.chi2, opt])
    bar.set_postfix(chi2="{0:.0f}".format(chi2))<jupyter_output>100%|██████████| 1000/1000 [00:39<00:00, 25.32it/s, chi2=621761]<jupyter_text>Here's the new time delay curve:<jupyter_code>plt.plot(times, model.run(model.tau)[:, 0], ".k", ms=1)
plt.ylabel(r"$\tau(t)$")
plt.xlabel("$t$");<jupyter_output><empty_output><jupyter_text>That looks interesting! Let's do one more round of optimization to make sure that we're converged to something.<jupyter_code>bar = tqdm.trange(1000)
for i in bar:
    chi2, _ = model.run([model.chi2, opt])
    bar.set_postfix(chi2="{0:.0f}".format(chi2))<jupyter_output>100%|██████████| 1000/1000 [00:40<00:00, 24.70it/s, chi2=392232]<jupyter_text>And here's the new time delay curve:<jupyter_code>plt.plot(times, model.run(model.tau)[:, 0], ".k", ms=1)
plt.ylabel(r"$\tau(t)$")
plt.xlabel("$t$");<jupyter_output><empty_output><jupyter_text>Cool. What if we zoom in?<jupyter_code>plt.plot(times, model.run(model.tau)[:, 0], ".k", ms=1)
plt.xlim(-600, -500)
plt.ylabel(r"$\tau(t)$")
plt.xlabel("$t$");<jupyter_output><empty_output><jupyter_text>That looks pretty legit. What about the PSD of the time delay curve?<jupyter_code>x, y = model.run([model.time, model.psi])
df = 0.5 / (times.max() - times.min())
freq = np.arange(df, 1.0, 0.1*df)
power = LombScargle(x, y).power(freq)

peaks = estimate_frequencies(x, y, max_peaks=2)
print(1.0 / peaks)

plt.loglog(freq, power, "k")
plt.ylabel("power")
plt.xlabel("frequency [1 / day]");<jupyter_output>INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
  Objective function value: 0.004375
  Number of iterations: 5
  Number of functions evaluations: 13
[342.28448333  10.49866768]<jupyter_text>This has signal at 346 days as previously reported, and a new signal at 10.5 days.Now let's fit for the primary Keplerian orbit as we were doing earlier this week:<jupyter_code>kep_model_1 = Maelstrom(times, dmmags, nu_arr[:nnu])
kep_model_1.init_from_orbit(period=1.0 / peaks[0], lighttime=np.std(y), tref=0.0)
kep_model_1.pin_lighttime_values()
for i in range(3):
    kep_model_1.optimize([kep_model_1.log_sigma2])
    kep_model_1.optimize([kep_model_1.lighttime, kep_model_1.tref])
    kep_model_1.optimize([kep_model_1.period, kep_model_1.lighttime, kep_model_1.tref])
    kep_model_1.optimize([kep_model_1.period, kep_model_1.lighttime, kep_model_1.tref,
                          kep_model_1.eccen_param, kep_model_1.varpi])
    kep_model_1.optimize([kep_model_1.nu])<jupyter_output>INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
  Objective function value: 281433.380446
  Number of iterations: 12
  Number of functions evaluations: 13
INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 251031.851312
  Number of iterations: 4
  Number of functions evaluations: 9
INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 247216.491114
  Number of iterations: 8
  Number of functions evaluations: 14
INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 247216.491114
  Number of iterations: 1
  Number of functions evaluations: 4
INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 246456.3[...]<jupyter_text>Then overplot the time delay curve on the one from before. Looks about right modulo some "slope" in the spline version.<jupyter_code>plt.plot(times, model.run(model.tau)[:, 0], ".k", ms=1)
plt.plot(times, kep_model_1.run(kep_model_1.tau)[:, 0], ".g", ms=1)
plt.ylabel(r"$\tau(t)$")
plt.xlabel("$t$");<jupyter_output><empty_output><jupyter_text>Now, let's subtract out the primary time delay and fit for the second signal.<jupyter_code>tau1 = kep_model_1.run(kep_model_1.tau)[:, 0]
kep_model_2 = Maelstrom(times - tau1, dmmags, nu_arr[:nnu])
kep_model_2.init_from_orbit(period=1.0 / peaks[1], lighttime=0.5*np.std(y), tref=0.0)
kep_model_2.pin_lighttime_values()
for i in range(3):
    kep_model_2.optimize([kep_model_2.log_sigma2])
    kep_model_2.optimize([kep_model_2.lighttime, kep_model_2.tref])
    kep_model_2.optimize([kep_model_2.period, kep_model_2.lighttime, kep_model_2.tref])
    kep_model_2.optimize([kep_model_2.period, kep_model_2.lighttime, kep_model_2.tref,
                          kep_model_2.eccen_param, kep_model_2.varpi, kep_model_2.tref])
    kep_model_2.optimize([kep_model_2.nu])<jupyter_output>INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 239097.610264
  Number of iterations: 10
  Number of functions evaluations: 11
INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 221753.894709
  Number of iterations: 21
  Number of functions evaluations: 28
INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 221753.149060
  Number of iterations: 2
  Number of functions evaluations: 6
INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 221753.149040
  Number of iterations: 1
  Number of functions evaluations: 4
INFO:tensorflow:Optimization terminated with:
  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
  Objective function value: 219959.5[...]<jupyter_text>This time, let's plot the time delay curve folded on itself. (Note: I've subtracted the first Keplerian fit out of the spline results).<jupyter_code>period = kep_model_2.run(kep_model_2.period)
plt.plot(times % period, model.run(model.tau)[:, 0] - tau1, ".k", ms=1)
plt.plot(times % period, kep_model_2.run(kep_model_2.tau)[:, 0], ".g", ms=1)
plt.ylabel(r"$\tau(t) - \tau_1(t)$")
plt.xlabel("orbital phase");<jupyter_output><empty_output><jupyter_text>Dashy cray. It looks like this orbit has the following parameters:<jupyter_code>print("period: {0:.5f} day".format(kep_model_2.run(kep_model_2.period)))
print("a1 sin(i) / c: {0:.5f} sec".format(86400.0 * np.abs(kep_model_2.run(kep_model_2.lighttime[0]))))
resid = dmmags - model.run(model.model)
x = (times - model.run(model.tau)[:, 0]) % period

bins = np.linspace(0, period, 200)
num, _ = np.histogram(x, bins, weights=resid)
denom, _ = np.histogram(x, bins)

plt.scatter(x, resid, c=times, s=2)
plt.plot(0.5*(bins[1:]+bins[:-1]), num / denom)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13778083447959652, 'fraction_numerical': 0.05055020632737277, 'mean_word_length': 6.053354890864996, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 55, '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': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '437009', 'license': 'MIT', 'n_tokens_mistral': 3759, 'n_tokens_neox': 3404, 'n_words': 685, 'path': 'notebooks/trinary.ipynb', 'repo_name': 'dfm/maelstrom', 'size': '8724'}