Note
This tutorial was generated from an IPython notebook that can be downloaded here.
# ref: https://radvel.readthedocs.io/en/latest/tutorials/K2-24_Fitting+MCMC.html
# ref: https://arxiv.org/abs/1511.04497
%matplotlib inline
%config InlineBackend.figure_format = "retina"
from matplotlib import rcParams
rcParams["savefig.dpi"] = 100
rcParams["figure.dpi"] = 100
rcParams["font.size"] = 20
import corner
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
from astropy.io import fits
from astropy.stats import BoxLeastSquares, LombScargle
from astropy import constants, units
import pymc3 as pm
import pymc3.distributions.transforms as tr
import theano
import theano.tensor as tt
import exoplanet
from exoplanet import distributions
from exoplanet.gp import terms, GP
import lightkurve
from lightkurve import search_targetpixelfile
lc_url = "https://archive.stsci.edu/hlsps/everest/v2/c02/203700000/71098/hlsp_everest_k2_llc_203771098-c02_kepler_v2.0_lc.fits"
with fits.open(lc_url) as hdus:
lc = hdus[1].data
lc_hdr = hdus[1].header
texp = lc_hdr["FRAMETIM"] * lc_hdr["NUM_FRM"]
texp /= 60.0 * 60.0 * 24.0
m = (np.arange(len(lc)) > 100) & np.isfinite(lc["FLUX"]) & np.isfinite(lc["TIME"])
bad_bits=[1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 13, 14, 16, 17]
qual = lc["QUALITY"]
for b in bad_bits:
m &= qual & 2 ** (b - 1) == 0
x = lc["TIME"][m]
y = lc["FLUX"][m]
mu = np.median(y)
y = (y / mu - 1) * 1e3
smooth = savgol_filter(y, 501, polyorder=5)
resid = y - smooth
sigma = np.sqrt(np.mean(resid**2))
m = resid < 3 * sigma
plt.plot(x, y, "k", label="data")
plt.plot(x, smooth, label="smoothed")
plt.plot(x[~m], y[~m], "xr", label="outliers")
plt.legend(fontsize=12)
plt.xlabel("time")
plt.ylabel("flux")
x_ref = np.min(x[m])
x = np.ascontiguousarray(x[m], dtype=np.float64)
x -= x_ref
y = np.ascontiguousarray(y[m], dtype=np.float64)
smooth = savgol_filter(y, 501, polyorder=5)
m = np.zeros(len(x), dtype=bool)
period_grid = np.exp(np.linspace(np.log(10), np.log(50), 50000))
bls_results = []
periods = []
t0s = []
depths = []
for i in range(2):
bls = BoxLeastSquares(x[~m], y[~m] - smooth[~m])
bls_power = bls.power(period_grid, 0.1, oversample=20)
bls_results.append(bls_power)
index = np.argmax(bls_power.power)
periods.append(bls_power.period[index])
t0s.append(bls_power.transit_time[index])
depths.append(bls_power.depth[index])
m |= bls.transit_mask(x, periods[-1], 0.5, t0s[-1])
fig, axes = plt.subplots(len(bls_results), 1, figsize=(10, 6), sharex=True)
for i in range(len(bls_results)):
ax = axes[i]
ax.axvline(np.log10(periods[i]), color="C1", lw=5, alpha=0.8)
ax.plot(np.log10(bls_results[i].period), bls_results[i].power, "k")
ax.set_ylabel("bls power")
ax.set_yticks([])
ax = axes[-1]
ax.set_xlim(np.log10(period_grid.min()), np.log10(period_grid.max()))
ax.set_xlabel("$\log_{10}$ period");
fig, axes = plt.subplots(len(bls_results), 1, figsize=(10, 6), sharex=True)
for i in range(len(bls_results)):
ax = axes[i]
p = periods[i]
x_fold = (x - t0s[i] + 0.5*p) % p - 0.5*p
ax.plot(x_fold, y - smooth, ".k")
bins = np.linspace(-0.4, 0.4, 20)
denom, _ = np.histogram(x_fold, bins)
num, _ = np.histogram(x_fold, bins, weights=y - smooth)
denom[num == 0] = 1.0
ax.plot(0.5*(bins[1:] + bins[:-1]), num / denom, color="C1")
ax = axes[-1]
ax.set_xlim(-0.4, 0.4)
ax.set_xlabel("time since transit");
N_planets = len(periods)
def build_model(mask=None):
if mask is None:
mask = np.ones(len(x), dtype=bool)
with pm.Model() as model:
# Stellar properties
mean = pm.Normal("mean", mu=0.0, sd=1.0)
u = distributions.Triangle("u", shape=2)
mstar = pm.Bound(pm.Normal, lower=0.0)("mstar", mu=1.12, sd=0.05)
rstar = pm.Bound(pm.Normal, lower=0.0)("rstar", mu=1.21, sd=0.11)
t0 = pm.Normal("t0", mu=t0s, sd=0.1, shape=N_planets)
logperiod = pm.Normal("logperiod", mu=np.log(periods), sd=1e-4, shape=N_planets)
period = pm.Deterministic("period", tt.exp(logperiod))
# Radius/impact parameter joint distribution
rb_test = 0.5 + np.zeros((2, N_planets))
rb_test[0, :] = np.sqrt(1e-3)*np.sqrt(depths)
rb = distributions.RadiusImpactParameter(
"rb", min_radius=0.01, max_radius=0.1,
shape=(2, N_planets), testval=rb_test)
ror = pm.Deterministic("ror", rb[0])
r = pm.Deterministic("r", ror * rstar)
b = pm.Deterministic("b", rb[1])
# Log-uniform prior on r
pm.Potential("logrprior", -tt.log(r))
ecc = pm.Beta("ecc", alpha=0.867, beta=3.03, shape=N_planets,
testval=np.array([0.1, 0.1]))
omega = distributions.Angle("omega", shape=N_planets)
# Set up the orbit
orbit = exoplanet.orbits.KeplerianOrbit(
period=period,
t0=t0,
ecc=ecc,
omega=omega,
r_star=rstar,
m_star=mstar,
b=b,
)
lc = exoplanet.StarryLightCurve(u, r_star=rstar)
transits = lc.get_light_curve(r, orbit, x[mask], texp=texp, oversample=5, order=2)
transit = tt.sum(transits, axis=-1) * 1e3
pm.Deterministic("transit", transit)
# Likelihood
logs2 = pm.Normal("logs2", mu=-5.0, sd=5.0)
logw0 = pm.Normal("logw0", mu=0.0, sd=5.0)
logQ = pm.Normal("logQ", mu=1.0, sd=5.0)
delta = pm.Normal("delta", mu=0.0, sd=10.0)
logS0 = pm.Deterministic("logS0", delta - logw0 - logQ)
kernel = terms.SHOTerm(S0=tt.exp(logS0), w0=tt.exp(logw0), Q=tt.exp(logQ))
gp = GP(kernel, x[mask], tt.exp(logs2) + np.zeros_like(x[mask]), J=2)
resid = y[mask] - transit - mean
pm.Potential("obs", gp.log_likelihood(resid))
gp_pred = gp.predict()
pm.Deterministic("gp", mean + gp_pred)
pm.Deterministic("resid", resid - gp_pred)
soln = pm.find_MAP(start=model.test_point, vars=[logs2, logS0, logw0, logQ, mean])
soln = pm.find_MAP(start=soln, vars=[logperiod, t0])
soln = pm.find_MAP(start=soln)
return model, soln
model0, soln0 = build_model()
logp = 2,010.2, ||grad|| = 0.17778: 100%|██████████| 35/35 [00:00<00:00, 63.94it/s]
logp = 2,204.1, ||grad|| = 0.22811: 100%|██████████| 42/42 [00:00<00:00, 49.04it/s]
logp = 4,841.4, ||grad|| = 978.65: 100%|██████████| 1249/1249 [00:32<00:00, 37.93it/s]
with model0:
result, func, args = exoplanet.utils.eval_in_model(model0.obs, return_func=True, profile=True)
func(*args)
array(68.33818922)
func.profile.summary()
Function profiling
==================
Message: /Users/dforeman/research/projects/dfm/exoplanet/exoplanet/utils.py:17
Time in 2 calls to Function.__call__: 1.670885e-02s
Time in Function.fn.__call__: 1.658106e-02s (99.235%)
Time in thunks: 1.629901e-02s (97.547%)
Total compile time: 1.678204e+00s
Number of Apply nodes: 144
Theano Optimizer time: 1.381569e+00s
Theano validate time: 3.431535e-02s
Theano Linker time (includes C, CUDA code generation/compiling): 1.467879e-01s
Import time 5.626464e-02s
Node make_thunk time 1.410229e-01s
Node sigmoid(InplaceDimShuffle{x,x,0}.0) time 1.939797e-02s
Node Elemwise{Composite{(i0 + (i1 * i2))}}(TensorConstant{(1, 1, 1) of 1.0}, sigmoid.0, Elemwise{cos,no_inplace}.0) time 1.499081e-02s
Node Elemwise{Composite{(i0 - (i1 * (Composite{(i0 * arctan2((sqrt((i1 - i2)) * i3), (sqrt((i1 + i2)) * (i1 + i4))))}(i2, i3, i4, i5, i6) - (scalar_sigmoid(i7) * sin(Composite{(i0 * arctan2((sqrt((i1 - i2)) * i3), (sqrt((i1 + i2)) * (i1 + i4))))}(i2, i3, i4, i5, i6)))) * i8))}}(InplaceDimShuffle{x,x,0}.0, TensorConstant{(1, 1, 1) ..4309189535}, TensorConstant{(1, 1, 1) of 2.0}, TensorConstant{(1, 1, 1) of 1.0}, sigmoid.0, Elemwise{cos,no_inplace}.0, Elemwise{Sin}[(0, 0)].0, InplaceDimShuffle{x,x,0}.0, Elemwise{exp,no_inplace}.0) time 1.413298e-02s
Node Elemwise{Composite{(i0 - sqr(i1))}}[(0, 1)](TensorConstant{(1, 1, 1) of 1.0}, sigmoid.0) time 1.326609e-02s
Node Elemwise{Mul}[(0, 0)](Elemwise{Composite{exp(((i0 - i1) - i2))}}.0, Elemwise{exp,no_inplace}.0, Elemwise{exp,no_inplace}.0) time 1.629114e-03s
Time in all call to theano.grad() 3.060994e+01s
Time since theano import 171.645s
Class
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>
49.8% 49.8% 0.008s 4.06e-03s C 2 1 exoplanet.theano_ops.kepler.solver.KeplerOp
33.8% 83.5% 0.006s 4.66e-05s C 118 67 theano.tensor.elemwise.Elemwise
6.6% 90.2% 0.001s 5.41e-04s C 2 1 exoplanet.theano_ops.starry.limbdark.LimbDarkOp
3.7% 93.8% 0.001s 4.99e-05s C 12 6 theano.tensor.basic.Join
1.4% 95.3% 0.000s 5.78e-05s C 4 2 theano.tensor.basic.Alloc
1.2% 96.5% 0.000s 8.83e-06s Py 22 6 theano.ifelse.IfElse
1.2% 97.6% 0.000s 9.61e-05s C 2 1 exoplanet.theano_ops.celerite.solve.SolveOp
1.0% 98.6% 0.000s 7.89e-05s C 2 1 exoplanet.theano_ops.celerite.factor.FactorOp
0.7% 99.3% 0.000s 1.38e-05s C 8 4 theano.tensor.elemwise.Sum
0.4% 99.7% 0.000s 1.19e-06s C 60 31 theano.tensor.elemwise.DimShuffle
0.2% 99.9% 0.000s 1.50e-06s C 22 11 theano.tensor.subtensor.Subtensor
0.0% 100.0% 0.000s 7.15e-07s C 8 4 theano.tensor.basic.Reshape
0.0% 100.0% 0.000s 1.07e-06s C 2 1 theano.compile.ops.Shape_i
0.0% 100.0% 0.000s 2.38e-07s C 8 4 theano.compile.ops.Rebroadcast
0.0% 100.0% 0.000s 9.54e-07s C 2 3 theano.tensor.opt.MakeVector
0.0% 100.0% 0.000s 9.54e-07s C 2 1 exoplanet.theano_ops.starry.get_cl.GetClOp
... (remaining 0 Classes account for 0.00%(0.00s) of the runtime)
Ops
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>
49.8% 49.8% 0.008s 4.06e-03s C 2 1 KeplerOp{tol=1e-08, maxiter=2000}
11.4% 61.2% 0.002s 9.28e-04s C 2 1 Elemwise{Composite{(sqrt((sqr(((((i0 * i1 * i2 * i3 * i4) / i5) - ((i0 * i6 * i2 * i3 * i7) / i5)) - i8)) + sqr(((i9 * i10) - i11)))) / i12)}}[(0, 4)]
6.6% 67.8% 0.001s 5.41e-04s C 2 1 LimbDarkOp
5.1% 72.9% 0.001s 1.03e-04s C 8 4 Elemwise{cos,no_inplace}
4.6% 77.4% 0.001s 9.27e-05s C 8 4 Elemwise{Sin}[(0, 0)]
3.9% 81.3% 0.001s 3.16e-04s C 2 1 Elemwise{Composite{(((i0 * i1 * i2 * i3 * i4) / i5) + ((i0 * i6 * i2 * i3 * i7) / i5))}}
3.7% 85.0% 0.001s 4.99e-05s C 12 6 Join
2.2% 87.1% 0.000s 1.78e-04s C 2 1 Elemwise{Composite{(((i0 * i1) - i2) / i3)}}[(0, 1)]
1.8% 89.0% 0.000s 4.93e-05s C 6 3 Elemwise{Composite{(i0 + (i1 * i2))}}
1.8% 90.8% 0.000s 1.46e-04s C 2 1 Elemwise{Composite{((i0 * (i1 - i2)) / i3)}}
1.4% 92.2% 0.000s 5.78e-05s C 4 2 Alloc
1.2% 93.4% 0.000s 8.83e-06s Py 22 6 if{inplace}
1.2% 94.5% 0.000s 9.61e-05s C 2 1 SolveOp{J=2, n_rhs=1}
1.0% 95.5% 0.000s 7.89e-05s C 2 1 FactorOp{J=2, n_rhs=-1}
0.9% 96.4% 0.000s 7.34e-05s C 2 1 Elemwise{Mul}[(0, 1)]
0.7% 97.1% 0.000s 1.78e-05s C 6 3 Elemwise{second,no_inplace}
0.6% 97.7% 0.000s 4.91e-05s C 2 1 Sum{axis=[1, 2], acc_dtype=float64}
0.3% 98.0% 0.000s 1.28e-05s C 4 2 Elemwise{Composite{exp((i0 * i1 * i2))}}
0.3% 98.2% 0.000s 2.05e-05s C 2 1 Elemwise{Composite{((i0 * i1) + log(i2))}}[(0, 0)]
0.2% 98.4% 0.000s 1.44e-06s C 20 10 Subtensor{int64}
... (remaining 49 Ops account for 1.59%(0.00s) of the runtime)
Apply
------
<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>
49.8% 49.8% 0.008s 4.06e-03s 2 62 KeplerOp{tol=1e-08, maxiter=2000}(Elemwise{Composite{((i0 * (i1 - i2)) / i3)}}.0, Alloc.0)
11.4% 61.2% 0.002s 9.28e-04s 2 100 Elemwise{Composite{(sqrt((sqr(((((i0 * i1 * i2 * i3 * i4) / i5) - ((i0 * i6 * i2 * i3 * i7) / i5)) - i8)) + sqr(((i9 * i10) - i11)))) / i12)}}[(0, 4)](TensorConstant{(1, 1, 1) of -1.0}, Elemwise{cos,no_inplace}.0, Elemwise{Composite{((i0 * i1 * sqr(i2)) ** i3)}}[(0, 2)].0, Elemwise{Composite{(i0 - sqr(i1))}}[(0, 1)].0, Elemwise{cos,no_inplace}.0, Elemwise{Composite{(i0 + (i1 * i2))}}.0, Elemwise{Sin}[(0, 0)].0, Elemwise{Sin}[(0, 0)].0, TensorConstan
6.6% 67.8% 0.001s 5.41e-04s 2 109 LimbDarkOp(Elemwise{Composite{((i0 * i1) / i2)}}[(0, 1)].0, Elemwise{Composite{(sqrt((sqr(((((i0 * i1 * i2 * i3 * i4) / i5) - ((i0 * i6 * i2 * i3 * i7) / i5)) - i8)) + sqr(((i9 * i10) - i11)))) / i12)}}[(0, 4)].0, Elemwise{second,no_inplace}.0, Elemwise{Composite{(((i0 * i1) - i2) / i3)}}[(0, 1)].0)
4.6% 72.3% 0.001s 3.71e-04s 2 67 Elemwise{cos,no_inplace}(KeplerOp{tol=1e-08, maxiter=2000}.1)
4.1% 76.4% 0.001s 3.31e-04s 2 71 Elemwise{Sin}[(0, 0)](KeplerOp{tol=1e-08, maxiter=2000}.1)
3.9% 80.3% 0.001s 3.16e-04s 2 88 Elemwise{Composite{(((i0 * i1 * i2 * i3 * i4) / i5) + ((i0 * i6 * i2 * i3 * i7) / i5))}}(TensorConstant{(1, 1, 1) of -1.0}, Elemwise{Sin}[(0, 0)].0, Elemwise{Composite{((i0 * i1 * sqr(i2)) ** i3)}}[(0, 2)].0, Elemwise{Composite{(i0 - sqr(i1))}}[(0, 1)].0, Elemwise{cos,no_inplace}.0, Elemwise{Composite{(i0 + (i1 * i2))}}.0, Elemwise{cos,no_inplace}.0, Elemwise{Sin}[(0, 0)].0)
2.2% 82.5% 0.000s 1.78e-04s 2 104 Elemwise{Composite{(((i0 * i1) - i2) / i3)}}[(0, 1)](Elemwise{Sin}[(0, 0)].0, Elemwise{Composite{(((i0 * i1 * i2 * i3 * i4) / i5) + ((i0 * i6 * i2 * i3 * i7) / i5))}}.0, TensorConstant{(1, 1, 1) of 0.0}, Elemwise{exp,no_inplace}.0)
2.0% 84.4% 0.000s 1.61e-04s 2 137 Join(TensorConstant{1}, Elemwise{second,no_inplace}.0, Elemwise{Composite{((i0 * i1) + (i2 * i3))}}.0, Elemwise{Composite{((i0 * i1) - (i2 * i3))}}[(0, 1)].0)
1.8% 86.2% 0.000s 1.46e-04s 2 55 Elemwise{Composite{((i0 * (i1 - i2)) / i3)}}(TensorConstant{(1, 1, 1) ..5307179586}, TensorConstant{[[[-1.0216..270e+01]]]}, Elemwise{Composite{(i0 - (i1 * (Composite{(i0 * arctan2((sqrt((i1 - i2)) * i3), (sqrt((i1 + i2)) * (i1 + i4))))}(i2, i3, i4, i5, i6) - (scalar_sigmoid(i7) * sin(Composite{(i0 * arctan2((sqrt((i1 - i2)) * i3), (sqrt((i1 + i2)) * (i1 + i4))))}(i2, i3, i4, i5, i6)))) * i8))}}.0, Elemwise{exp,no_inplace}.0)
1.8% 88.0% 0.000s 1.46e-04s 2 76 Elemwise{Composite{(i0 + (i1 * i2))}}(TensorConstant{(1, 1, 1) of 1.0}, sigmoid.0, Elemwise{cos,no_inplace}.0)
1.4% 89.4% 0.000s 1.13e-04s 2 134 Join(TensorConstant{1}, Elemwise{second,no_inplace}.0, Elemwise{cos,no_inplace}.0, Elemwise{Sin}[(0, 0)].0)
1.4% 90.8% 0.000s 1.11e-04s 2 30 Alloc(sigmoid.0, TensorConstant{3515}, TensorConstant{5}, Shape_i{0}.0)
1.2% 92.0% 0.000s 9.61e-05s 2 139 SolveOp{J=2, n_rhs=1}(Join.0, Join.0, FactorOp{J=2, n_rhs=-1}.0, FactorOp{J=2, n_rhs=-1}.1, InplaceDimShuffle{0,x}.0)
1.0% 92.9% 0.000s 7.89e-05s 2 138 FactorOp{J=2, n_rhs=-1}(Alloc.0, Join.0, Join.0, Join.0)
0.9% 93.8% 0.000s 7.34e-05s 2 114 Elemwise{Mul}[(0, 1)](TensorConstant{[[[0.08333..8333333]]]}, LimbDarkOp.0)
0.6% 94.5% 0.000s 4.99e-05s 2 89 Elemwise{second,no_inplace}(TensorConstant{(1, 3515, .. 1) of 0.0}, InplaceDimShuffle{x,x,x,0}.0)
0.6% 95.1% 0.000s 4.91e-05s 2 119 Sum{axis=[1, 2], acc_dtype=float64}(Elemwise{Mul}[(0, 1)].0)
0.5% 95.5% 0.000s 3.90e-05s 2 133 Elemwise{Sin}[(0, 0)](Elemwise{mul,no_inplace}.0)
0.5% 96.0% 0.000s 3.90e-05s 2 131 Elemwise{cos,no_inplace}(Elemwise{mul,no_inplace}.0)
0.3% 96.4% 0.000s 1.41e-05s 4 74 if{inplace}(Elemwise{lt,no_inplace}.0, Elemwise{Mul}[(0, 4)].0, TensorConstant{[]})
... (remaining 124 Apply instances account for 3.64%(0.00s) of the runtime)
Here are tips to potentially make your code run faster
(if you think of new ones, suggest them on the mailing list).
Test them first, as they are not guaranteed to always provide a speedup.
- Try the Theano flag floatX=float32
We don't know if amdlibm will accelerate this scalar op. arctan2
We don't know if amdlibm will accelerate this scalar op. arccos
- Try installing amdlibm and set the Theano flag lib.amdlibm=True. This speeds up only some Elemwise operation.
with model0:
# Sigma clipping
resid = exoplanet.utils.eval_in_model(model0.resid, soln0)
sigma = np.sqrt(np.mean(resid**2))
mask = np.abs(resid) < 4 * sigma
plt.plot(x, resid)
plt.plot(x[~mask], resid[~mask], "xr")
model, soln = build_model(mask)
logp = 2,607.9, ||grad|| = 0.0023071: 100%|██████████| 42/42 [00:00<00:00, 59.56it/s]
logp = 3,120, ||grad|| = 663.87: 100%|██████████| 35/35 [00:00<00:00, 50.02it/s]
logp = 5,311.1, ||grad|| = 545.25: 100%|██████████| 471/471 [00:11<00:00, 39.94it/s]
with model:
plt.plot(x[mask], y[mask], ".k")
plt.plot(x[mask], exoplanet.utils.eval_in_model(model.transit, soln))
plt.plot(x[mask], exoplanet.utils.eval_in_model(model.gp, soln))
plt.figure()
# plt.plot(x[mask], exoplanet.utils.eval_in_model(model.resid))
plt.plot(x[mask], exoplanet.utils.eval_in_model(model.resid, soln))
schedule = exoplanet.sampling.TuningSchedule(start=100, window=50)
with model:
burnin = schedule.tune(tune=2000, start=soln, step_kwargs=dict(regular_window=0))
Only 2 samples in chain.
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [logQ, logw0, logS0, logs2, omega, ecc, rb, logperiod, t0, rstar, mstar, u, mean]
Sampling 2 chains: 100%|██████████| 204/204 [06:37<00:00, 6.30s/draws]
The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.
The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.
Only 2 samples in chain.
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [logQ, logw0, logS0, logs2, omega, ecc, rb, logperiod, t0, rstar, mstar, u, mean]
Sampling 2 chains: 100%|██████████| 104/104 [03:49<00:00, 5.69s/draws]
The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.
The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.
Only 2 samples in chain.
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [logQ, logw0, logS0, logs2, omega, ecc, rb, logperiod, t0, rstar, mstar, u, mean]
Sampling 2 chains: 100%|██████████| 204/204 [06:01<00:00, 1.87s/draws]
Only 2 samples in chain.
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [logQ, logw0, logS0, logs2, omega, ecc, rb, logperiod, t0, rstar, mstar, u, mean]
Sampling 2 chains: 100%|██████████| 404/404 [12:19<00:00, 8.14s/draws]
The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.
Only 2 samples in chain.
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [logQ, logw0, logS0, logs2, omega, ecc, rb, logperiod, t0, rstar, mstar, u, mean]
Sampling 2 chains: 100%|██████████| 804/804 [41:33<00:00, 8.78s/draws]
Only 2 samples in chain.
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [logQ, logw0, logS0, logs2, omega, ecc, rb, logperiod, t0, rstar, mstar, u, mean]
Sampling 2 chains: 100%|██████████| 2304/2304 [2:14:21<00:00, 6.06s/draws]
plt.plot(burnin["r"][:, 1])
[<matplotlib.lines.Line2D at 0x1c3cfefeb8>]
samples = pm.trace_to_dataframe(burnin, varnames=["u", "rstar", "mstar"])
samples.columns = [k.replace("_", " ") for k in samples.columns]
corner.corner(samples);
samples = pm.trace_to_dataframe(burnin, varnames=["t0", "period", "r", "b", "ecc", "omega"])
samples.columns = [k.replace("_", " ") for k in samples.columns]
corner.corner(samples);
samples = pm.trace_to_dataframe(burnin, varnames=["mean", "logs2", "logS0", "logw0", "logQ"])
samples.columns = [k.replace("_", " ") for k in samples.columns]
corner.corner(samples);
0.5*self.S0*self.w0*Q*tt.stack([1.0+1.0/f, 1.0-1.0/f]),
samples = pm.trace_to_dataframe(burnin, varnames=["logS0", "logw0", "logQ"])
samples["delta"] = samples.logS0 + samples.logw0 + samples.logQ
corner.corner(samples);
with model:
trace = schedule.sample(draws=2000)
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [logQ, logw0, logS0, logs2, omega, ecc, rb, logperiod, t0, rstar, mstar, u, mean]
Sampling 2 chains: 80%|████████ | 3285/4100 [5:24:51<3:42:30, 16.38s/draws]
ERROR:root:Internal Python error in the inspect module.
Below is the traceback from this internal error.
Traceback (most recent call last):
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2881, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "<ipython-input-14-6b9cae1d1b36>", line 2, in <module>
trace = schedule.sample(draws=2000)
File "/Users/dforeman/research/projects/dfm/exoplanet/exoplanet/sampling.py", line 160, in sample
return pm.sample(start=start, step=step, **kwargs)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/pymc3/sampling.py", line 472, in sample
trace = trace[discard:]
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/pymc3/backends/base.py", line 292, in __getitem__
return self._slice(idx)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/pymc3/backends/base.py", line 496, in _slice
new_traces = [trace._slice(slice) for trace in self._straces.values()]
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/pymc3/backends/base.py", line 496, in <listcomp>
new_traces = [trace._slice(slice) for trace in self._straces.values()]
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/pymc3/backends/ndarray.py", line 287, in _slice
sliced = NDArray(model=self.model, vars=self.vars)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/pymc3/backends/ndarray.py", line 160, in __init__
super(NDArray, self).__init__(name, model, vars, test_point)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/pymc3/backends/base.py", line 50, in __init__
self.fn = model.fastfn(vars)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/pymc3/model.py", line 939, in fastfn
f = self.makefn(outs, mode, *args, **kwargs)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/pymc3/model.py", line 909, in makefn
mode=mode, *args, **kwargs)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/compile/function.py", line 317, in function
output_keys=output_keys)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/compile/pfunc.py", line 486, in pfunc
output_keys=output_keys)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/compile/function_module.py", line 1839, in orig_function
name=name)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/compile/function_module.py", line 1519, in __init__
optimizer_profile = optimizer(fgraph)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/gof/opt.py", line 108, in __call__
return self.optimize(fgraph)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/gof/opt.py", line 97, in optimize
ret = self.apply(fgraph, *args, **kwargs)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/gof/opt.py", line 251, in apply
sub_prof = optimizer.optimize(fgraph)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/gof/opt.py", line 97, in optimize
ret = self.apply(fgraph, *args, **kwargs)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/gof/opt.py", line 2513, in apply
lopt_change = self.process_node(fgraph, node, lopt)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/gof/opt.py", line 2074, in process_node
remove=remove)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/gof/toolbox.py", line 569, in replace_all_validate_remove
chk = fgraph.replace_all_validate(replacements, reason)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/gof/toolbox.py", line 518, in replace_all_validate
fgraph.replace(r, new_r, reason=reason, verbose=False)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/gof/fg.py", line 514, in replace
self.change_input(node, i, new_r, reason=reason)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/gof/fg.py", line 452, in change_input
r, new_r, reason=reason)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/gof/fg.py", line 594, in execute_callbacks
fn(self, *args, **kwargs)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/tensor/opt.py", line 1380, in on_change_input
self.update_shape(new_r, r)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/tensor/opt.py", line 1244, in update_shape
for i in xrange(r.ndim)])
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/theano/tensor/opt.py", line 1244, in <listcomp>
for i in xrange(r.ndim)])
KeyboardInterrupt
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 1821, in showtraceback
stb = value._render_traceback_()
AttributeError: 'KeyboardInterrupt' object has no attribute '_render_traceback_'
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/IPython/core/ultratb.py", line 1132, in get_records
return _fixed_getinnerframes(etb, number_of_lines_of_context, tb_offset)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/IPython/core/ultratb.py", line 313, in wrapped
return f(*args, **kwargs)
File "/Users/dforeman/anaconda/lib/python3.6/site-packages/IPython/core/ultratb.py", line 358, in _fixed_getinnerframes
records = fix_frame_records_filenames(inspect.getinnerframes(etb, context))
File "/Users/dforeman/anaconda/lib/python3.6/inspect.py", line 1453, in getinnerframes
frameinfo = (tb.tb_frame,) + getframeinfo(tb, context)
File "/Users/dforeman/anaconda/lib/python3.6/inspect.py", line 1411, in getframeinfo
filename = getsourcefile(frame) or getfile(frame)
File "/Users/dforeman/anaconda/lib/python3.6/inspect.py", line 666, in getsourcefile
if getattr(getmodule(object, filename), '__loader__', None) is not None:
File "/Users/dforeman/anaconda/lib/python3.6/inspect.py", line 709, in getmodule
f = getabsfile(module)
File "/Users/dforeman/anaconda/lib/python3.6/inspect.py", line 678, in getabsfile
_filename = getsourcefile(object) or getfile(object)
File "/Users/dforeman/anaconda/lib/python3.6/inspect.py", line 663, in getsourcefile
if os.path.exists(filename):
File "/Users/dforeman/anaconda/lib/python3.6/genericpath.py", line 19, in exists
os.stat(path)
KeyboardInterrupt
---------------------------------------------------------------------------
samples = pm.trace_to_dataframe(trace, varnames=["u", "rstar", "mstar"])
samples.columns = [k.replace("_", " ") for k in samples.columns]
corner.corner(samples);
samples = pm.trace_to_dataframe(trace, varnames=["t0", "period", "r", "b", "ecc", "omega"])
samples.columns = [k.replace("_", " ") for k in samples.columns]
corner.corner(samples);
samples = pm.trace_to_dataframe(trace, varnames=["mean", "logs2", "logS0", "logw0", "logQ"])
samples.columns = [k.replace("_", " ") for k in samples.columns]
corner.corner(samples);
with model:
plt.plot(x[mask], y[mask] - np.median(trace["gp"], axis=0), ".k")
plt.plot(x[mask], np.median(trace["transit"], axis=0))
# plt.plot(x[mask], exoplanet.utils.eval_in_model(model.transit, soln))
# plt.plot(x[mask], exoplanet.utils.eval_in_model(model.gp, soln))
# plt.figure()
# # plt.plot(x[mask], exoplanet.utils.eval_in_model(model.resid))
# plt.plot(x[mask], exoplanet.utils.eval_in_model(model.resid, soln))