nespinoza / juliet

A versatile modelling tool for transiting and non-transiting (single and multiple) exoplanetary systems
MIT License
55 stars 30 forks source link

GP-only fit on photometry not working #104

Closed melissa-hobson closed 1 year ago

melissa-hobson commented 1 year ago

I'm trying to fit a GP (only) on photometric data, as done in the docs, and it crashes with the following error:

imagen imagen

Any idea why? I know this used to work, I can fit GP+transit simultaneously, and GP-only fits on RV data - it's just the GP-only on photometric data that's failing.

nespinoza commented 1 year ago

Hi @melissa-hobson,

Can you please provide a minimum working example that reproduces the error (perhaps internally)? That could help me pinpoint the source. From the looks of it, it seems as if the sigma_w_TESS was missing --- but it could also be an issue I introduced with the latest juliet updates.

N.

melissa-hobson commented 1 year ago

A missing sigma_w_TESS was my first thought, but it's definitely there in the priors, both in my own code and in the tutorial (which gives me the same error when just copy-pasted into a script and run, so I've copied that below):

import juliet
import numpy as np
import matplotlib.pyplot as plt

# First, get arrays of times, normalized-fluxes and errors for HATS-46
#from Sector 1 from MAST:
t, f, ferr  = juliet.get_TESS_data('https://archive.stsci.edu/hlsps/'+\
                                   'tess-data-alerts/hlsp_tess-data-'+\
                                   'alerts_tess_phot_00281541555-s01_'+\
                                   'tess_v1_lc.fits')
# Plot the data:
plt.errorbar(t,f,yerr=ferr,fmt='.')
plt.xlim([np.min(t),np.max(t)])
plt.xlabel('Time (BJD - 2457000)')
plt.ylabel('Relative flux')

# Period and t0 from Anderson et al. (201X):
P,t0 =  4.7423729 ,  2457376.68539 - 2457000
# Get phases --- identify out-of-transit (oot) times by phasing the data
# and selecting all points at absolute phases larger than 0.02:
phases = juliet.utils.get_phases(t, P, t0)
idx_oot = np.where(np.abs(phases)>0.02)[0]
# Save the out-of-transit data into dictionaries so we can feed them to juliet:
times, fluxes, fluxes_error = {},{},{}
times['TESS'], fluxes['TESS'], fluxes_error['TESS'] = t[idx_oot],f[idx_oot],ferr[idx_oot]

priors = {}
# Set the priors:
params =  ['mdilution_TESS', 'mflux_TESS', 'sigma_w_TESS', 'GP_sigma_TESS', \
           'GP_rho_TESS']
dists =   ['fixed',          'normal',     'loguniform',   'loguniform',\
           'loguniform']
hyperps = [1., [0.,0.1], [1e-6, 1e6], [1e-6, 1e6],\
           [1e-3,1e3]]

for param, dist, hyperp in zip(params, dists, hyperps):
    priors[param] = {}
    priors[param]['distribution'], priors[param]['hyperparameters'] = dist, hyperp

# Perform the juliet fit. Load dataset first (note the GP regressor will be the times):
dataset = juliet.load(priors=priors, t_lc = times, y_lc = fluxes, \
                      yerr_lc = fluxes_error, GP_regressors_lc = times, \
                      out_folder = 'hats46_detrending')
# # Fit:
results = dataset.fit()
nespinoza commented 1 year ago

Thanks for that @melissa-hobson! I found the bug and this should work on the latest version I just released, v.2.2.4 --- can you do pip install juliet --upgrade and check this is the case?

Thanks!

N.

melissa-hobson commented 1 year ago

Yes, that works! Thank you :)

nespinoza commented 1 year ago

Solved in v.2.2.4