Closed avivajpeyi closed 2 years ago
Do we know of an example TOI that I can try out?
A good place to start might be to find examples with period < 0 in the catalog. I did start implementing single transits at one point if you want to take a look there as a starting point: https://github.com/exoplanet-dev/tess.world/blob/main/src/tess_world/templates/template.py#L305
From https://exofop.ipac.caltech.edu/tess/view_toi.php
, I found 66/2394 candidates have period < 0
. Going to test PE with TOI 2180
.
single transit --> most likely a short period single transit
Hence incorporate that into the prior --> steep power-law prior in period
Re un with above, period, duration, b, probably very correlated. We can look at corner and see whats going on (Plot log_period)
Plot divergences on top of the corner to see where and which parameters might be causing the issue. Maybe this always happens at particular parameters?
Check if this works now that I have merged in #76
Reran 2180:
failed when testing the model https://github.com/dfm/tess-atlas/blob/81db5436b8f86f9da758e6272eaffa02e7890301/src/tess_atlas/notebook_templates/toi_template.py#L330
There is an issue with the model. The model(testval) has a nan:
dur_interval__ NaN
r_log__ -0.92
b_impact__ -1.39
t0_interval__ -1.19
p_1_lowerbound__ -inf
f0 -3.22
u_quadlimbdark__ -2.77
jitter_log__ 0.11
sigma_log__ 0.11
rho_log__ -0.55
obs -362020.84
The error is arising from the transit duration dur
and the orbital period p
.
{'TOI': 2180.01,
'Period (days)': 0.0,
'Epoch (TBJD)': 1830.7900000000373,
'Depth (ppt)': 2.4,
'Duration (days)': 1.3089166666666667,
'Single Transit': True}
def build_planet_transit_model(tic_entry):
t = tic_entry.lightcurve.time
y = tic_entry.lightcurve.flux
yerr = tic_entry.lightcurve.flux_err
n = tic_entry.planet_count
t0s = np.array([planet.t0 for planet in tic_entry.candidates])
depths = np.array([planet.depth for planet in tic_entry.candidates])
periods = np.array([planet.period for planet in tic_entry.candidates])
tmaxs = np.array([planet.tmax for planet in tic_entry.candidates])
durations = np.array([planet.duration for planet in tic_entry.candidates])
max_durations = np.array(
[planet.duration_max for planet in tic_entry.candidates]
)
min_durations = np.array(
[planet.duration_min for planet in tic_entry.candidates]
)
with pm.Model() as my_planet_transit_model:
## define planet parameters
# 1) d: transit duration (duration of eclipse)
d_priors = pm.Bound(
pm.Lognormal, lower=min_durations, upper=max_durations
)(
name=DURATION,
mu=np.log(durations),
sigma=np.log(1.2),
shape=n,
testval=durations,
)
...
# 1) t0: the time of the first transit in data (a reference time)
t0_norm = pm.Bound(
pm.Normal, lower=t0s - max_durations, upper=t0s + max_durations
)
t0_priors = t0_norm(
TIME_START, mu=t0s, sigma=0.5 * durations, shape=n, testval=t0s
)
# 2) period: the planets' orbital period
p_params, p_priors_list, tmax_priors_list = [], [], []
for n, planet in enumerate(tic_entry.candidates):
# if only one transit in data we use the period
if planet.has_data_only_for_single_transit:
p_prior = pm.Pareto(
name=f"{ORBITAL_PERIOD}_{planet.index}",
m=planet.period_min,
alpha=2.0 / 3.0,
testval=planet.period_min,
)
p_param = p_prior
tmax_prior = planet.t0
...
p_priors = pm.Deterministic(ORBITAL_PERIOD, tt.stack(p_priors_list))
tmax_priors = pm.Deterministic(TIME_END, tt.stack(tmax_priors_list))
...
LOL
tic_entry.candidates[0].num_periods
>>> -9223372036854775808
With the above changes, I was able to get the model to compile and run on the test val (no nans! woohoo!)
However, got the darn celerite2.backprop.LinAlgError: failed to factorize or solve matrix
error again!
Hmm looking at the logs, I see that there are two 0s in the point being optimized:
{'dur_interval__': array([4.59231013]),
'r_log__': array([-3.01614327]),
'b_impact__': array([-0.09347175]),
't0_interval__': array([0.]),
'p_1_lowerbound__': array(4.99338746),
'f0': array(0.),
'u_quadlimbdark__': array([4.4408921e-16, 0.0000000e+00]),
'jitter_log__': array(-0.70142785),
'sigma_log__': array(19.9134664),
'rho_log__': array(30.25482971)}
The error is being thrown at
theta = pmx.optimize(**kwargs, vars=noise_params)
and
noise_params = [jitter_prior, sigma_prior, rho_prior]
run_toi(1847)
results with an error at the start of the sampling step, because p_1_lowerbound -inf
An error occurred while executing the following cell:
------------------
if tic_entry.inference_data is None:
inference_data = run_inference(planet_transit_model)
else:
logger.info("Using cached run")
inference_data = tic_entry.inference_data
inference_data
------------------
---------------------------------------------------------------------------
SamplingError Traceback (most recent call last)
/tmp/ipykernel_176089/141900171.py in <module>
1 if tic_entry.inference_data is None:
----> 2 inference_data = run_inference(planet_transit_model)
3 else:
4 logger.info("Using cached run")
5 inference_data = tic_entry.inference_data
/tmp/ipykernel_176089/3349958539.py in run_inference(model)
4 sampling_kwargs = dict(tune=2000, draws=2000, chains=2, cores=2)
5 logger.info(f"Run sampler with kwargs: {sampling_kwargs}")
----> 6 inference_data = pmx.sample(
7 **sampling_kwargs, start=init_params, return_inferencedata=True
8 )
/fred/oz200/avajpeyi/envs/tess/lib/python3.8/site-packages/pymc3_ext/sampling/sampling.py in sample(draws, tune, model, step_kwargs, warmup_window, adapt_window, cooldown_window, initial_accept, target_accept, gamma, k, t0, parameter_groups, adapt_type, recompute_interval, regularization_steps, regularization_variance, **kwargs)
95 )
96
---> 97 return pm.sample(draws=draws, tune=tune, model=model, step=step, **kwargs)
/fred/oz200/avajpeyi/envs/tess/lib/python3.8/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, jitter_max_retries, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs)
433 for chain_start_vals in start:
434 update_start_vals(chain_start_vals, model.test_point, model)
--> 435 check_start_vals(start, model)
436
437 if cores is None:
/fred/oz200/avajpeyi/envs/tess/lib/python3.8/site-packages/pymc3/util.py in check_start_vals(start, model)
235
236 if not np.all(np.isfinite(initial_eval)):
--> 237 raise SamplingError(
238 "Initial evaluation of model at starting point failed!\n"
239 "Starting values:\n{}\n\n"
SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'dur_interval__': array([-1.75828217, 0.66137893]), 'r_log__': array([-2.3121765, -2.9045715]), 'b_impact__': array([-0.18072936, -0.10394958]), 't0_interval__': array([0., 0.]), 'p_1_lowerbound__': array(-inf), 'tmax_2_interval__': array(0.), 'f0': array(0.), 'u_quadlimbdark__': array([4.4408921e-16, 0.0000000e+00]), 'jitter_log__': array(0.5556663), 'sigma_log__': array(0.5556663), 'rho_log__': array(0.05250777), 'dur': array([0.149375 , 0.08883333]), 'r': array([0.09904544, 0.05477226]), 'b': array([0.5, 0.5]), 't0': array([1390.70655 , 2130.150082]), 'p_1': array(750.87597302), 'tmax_2': array(2130.150082), 'p': array([750.87597302, nan]), 'tmax': array([1390.70655 , 2130.150082]), 'u': array([0.70710678, 0. ]), 'jitter': array(1.74310202), 'sigma': array(1.74310202), 'rho': array(1.05391075), 'rho_circ': array([128.85573981, nan])}
Initial evaluation results:
dur_interval__ 0.28
r_log__ -1.84
b_impact__ -2.78
t0_interval__ 0.47
p_1_lowerbound__ -inf
tmax_2_interval__ -0.51
f0 -3.22
u_quadlimbdark__ -2.77
jitter_log__ 0.11
sigma_log__ 0.11
rho_log__ -0.55
obs -37873.84
Name: Log-probability of test_point, dtype: float64
SamplingError: Initial evaluation of model at starting point failed!
So the issue might be in the optimise_init_param
step?
From the flux vs time plot, looks like we've just missed the transit?
notebook: toi_1847.ipynb.zip
^ the test_model
function should also check for inf
s, not just isnull
Looks like im missing the transits for several single-transit TOIs!
Here are some examples:
Another set of funky errors -- the transit duration seems wrong?
The 77 analyses for single-transit systems have p_n_lowerbound = -inf
Eg from 1812:
The test-val is wrong -- currently set to prior's lower bound!
So although the periods are super long with even wider uncertainties -- the phase plots look ok?
I wonder if there is any more physics that we can use to constrain the period a bit more for these single-transit systems
Some more diagnostics:
Hhmm why are my ranges going so far out I wonder... i probably samples really far away from the 1,2,3 sigma contours
Finally, the trace plot:
looks like there are a few divergences
For single transits:
period -- are we just getting back the prior?
This could be a cool system to talk about + followup + maybe highlight this in the paper
There might be a way to deal with single transits now!
assert np.all(periods > 0), "We haven't implemented single transits yet"