<jupyter_start><jupyter_text>The Sample<jupyter_code>muirhead_data = pd.read_csv("datafiles/Muirhead2013_isochrones/muirhead_data_incmissing.txt", sep=" ")

# ALL Kepler planets from exo archive
planets = pd.read_csv('datafiles/exoplanetarchive/cumulative_kois.csv')

# Take the Kepler planet archive entries for the planets in Muirhead et al. 2013 sample
spectplanets = pd.read_csv('datafiles2/spectplanets.csv')

# Kepler-Gaia Data
kpgaia = Table.read('datafiles/Kepler-Gaia/kepler_dr2_4arcsec.fits', format='fits').to_pandas();

# Combined spectroscopy data + Gaia/Kepler data for our sample
muirhead_comb = pd.read_csv('datafiles2/muirhead_comb.csv')

# Only targets from table above with published luminosities from Gaia
muirhead_comb_lums = pd.read_csv('datafiles2/muirhead_comb_lums.csv')<jupyter_output>/Users/ssagear/anaconda3/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3418: TableReplaceWarning: converted column 'r_result_flag' from integer to float
  exec(code_obj, self.user_global_ns, self.user_ns)
/Users/ssagear/anaconda3/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3418: TableReplaceWarning: converted column 'r_modality_flag' from integer to float
  exec(code_obj, self.user_global_ns, self.user_ns)
/Users/ssagear/anaconda3/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3418: TableReplaceWarning: converted column 'teff_err1' from integer to float
  exec(code_obj, self.user_global_ns, self.user_ns)
/Users/ssagear/anaconda3/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3418: TableReplaceWarning: converted column 'teff_err2' from integer to float
  exec(code_obj, self.user_global_ns, self.user_ns)<jupyter_text>Defining a "test planet"<jupyter_code># Kepler ID for Kepler-1582 b
kepid = 5868793
kepname = spectplanets.loc[spectplanets['kepid'] == kepid].kepler_name.values[0]
kepoi_name = spectplanets.loc[spectplanets['kepid'] == kepid].kepoi_name.values[0]

kp1582b = muirhead_comb.loc[muirhead_comb['KIC'] == kepid]
isodf = pd.read_csv("datafiles/isochrones/iso_lums_" + str(kepid) + ".csv")
mstar = isodf["mstar"].mean()
mstar_err = isodf["mstar"].std()

rstar = isodf["radius"].mean()
rstar_err = isodf["radius"].std()
rho_star, mass, radius = ph.find_density_dist_symmetric(mstar, mstar_err, rstar, rstar_err, arrlen)<jupyter_output><empty_output><jupyter_text>Creating a fake light curve based on a real planet<jupyter_code>period, period_uerr, period_lerr, rprs, rprs_uerr, rprs_lerr, a_arc, a_uerr_arc, a_lerr_arc, inc, e_arc, w_arc = ph.planet_params_from_archive(spectplanets, kepoi_name)

# We calculate a_rs to ensure that it's consistent with the spec/Gaia stellar density.
a_rs = ph.calc_a(period*86400.0, mstar*smass_kg, rstar*srad_m)
a_rs_err = np.mean((a_uerr_arc, abs(a_lerr_arc)))

print('Stellar mass (Msun): ', mstar, 'Stellar radius (Rsun): ', rstar)
print('Period (Days): ', period, 'Rp/Rs: ', rprs)
print('a/Rs: ', a_rs)
print('i (deg): ', inc)
a_lerr_arc
perioderr = np.nanmean((abs(period_uerr), abs(period_lerr)))
rprserr = np.nanmean((abs(rprs_uerr), abs(rprs_lerr)))<jupyter_output><empty_output><jupyter_text>--------------------------------------------------------------------------------------  $e = 0.0$, $\omega = 90.0$  --------------------------------------------------------------------------------------<jupyter_code># 30 minute cadence
cadence = 0.02142857142857143

time = np.arange(-25, 25, cadence)
# Calculate flux from transit model
e = 0.0
w = 90.0
flux = ph.integratedlc(time, period, rprs, a_rs, e, inc, w, 0)

# Adding some gaussian noise
noise = np.random.normal(0,0.00006,len(time))
nflux = flux+noise

flux_err = np.array([0.00006]*len(nflux))<jupyter_output><empty_output><jupyter_text>Fitting the transit<jupyter_code>periodPDF = ph.get_period_dist(time, nflux, 4, 6, arrlen)
print('Period fit: ', ph.mode(periodPDF))
pdist = periodPDF
per_guess = ph.mode(pdist)
ttimes = np.concatenate((-np.arange(0, time[-1], period)[1:], np.arange(0, time[-1], period)))
ttimes = np.sort(ttimes)

time, flux, flux_err = ph.get_transit_cutout_full(ttimes, 4, time, nflux, flux_err)

t0 = 0.0<jupyter_output><empty_output><jupyter_text>Dynesty Transit Fit Likelihood Functions<jupyter_code>u = np.array([0,1])
a_rs_err
u*10*a_rs_err#+a_rs-((10*a_rs_err)/2.)
ndim = 5

def tfit_loglike(theta):
    """
    Transit fit dynesty function

    model = ph.integratedlc_fitter()
    gerr = sigma of g distribution

    """

    perl, rprsl, a_rsl, incl, t0l = theta

    model = ph.integratedlc_fitter(time, perl, rprsl, a_rsl, incl, t0l)
    sigma2 = flux_err ** 2

    return -0.5 * np.sum((flux - model) ** 2 / sigma2 + np.log(sigma2))


def tfit_prior_transform(utheta):
    """Transforms samples `u` drawn from the unit cube to samples to those
    from our uniform prior within [-10., 10.) for each variable."""

    uper, urp, ua, uinc, ut0 = utheta

    perl = uper*10*perioderr+period-((10*perioderr)/2.)
    rprsl = urp*10*rprserr+rprs-((10*rprserr)/2.)
    a_rsl = ua*10*a_rs_err+a_rs-((10*a_rs_err)/2.)
    incl = uinc*0.5+inc-0.25
    t0l = ut0*0.01-0.005

    return perl, rprsl, a_rsl, incl, t0l<jupyter_output><empty_output><jupyter_text>Run dynesty transit fit<jupyter_code>dsampler = dynesty.DynamicNestedSampler(tfit_loglike, tfit_prior_transform, ndim=ndim, bound='balls', sample='rwalk', nlive=1500, update_interval=700)
dsampler.run_nested(dlogz_init=0.05)
dres = dsampler.results
# Corner Plot of Transit Fit
fig, axes = dyplot.cornerplot(dres, labels=["period", "Rp/Rs", "a/Rs", "i", "t0"])
# plotting the static run
fig, axes = dyplot.traceplot(dres, truths=np.zeros(5),
                             show_titles=True, trace_cmap='plasma',
                             quantiles=None)
# Thinning distributions to size arrlen

pdist = random.choices(dres.samples[:,0], k=arrlen)
rdist = random.choices(dres.samples[:,1], k=arrlen)
adist = random.choices(dres.samples[:,2], k=arrlen)
idist = random.choices(dres.samples[:,3], k=arrlen)
t0dist = random.choices(dres.samples[:,4], k=arrlen)

per_f = ph.mode(pdist)
rprs_f = ph.mode(rdist)
a_f = ph.mode(adist)
i_f = ph.mode(idist)
t0_f = ph.mode(t0dist)
# Create a light curve with the fit parameters
fit1 = ph.integratedlc_fitter(time, per_f, rprs_f, a_f, i_f, t0)
t0_f
plt.errorbar(time, flux, yerr=flux_err, c='blue', alpha=0.5, label='Original LC')
plt.plot(time, fit1, c='red', alpha=1.0, label='Fit LC')
#plt.xlim(-0.1, 0.1)
plt.legend()

print('Stellar mass (Msun): ', mstar, 'Stellar radius (Rsun): ', rstar)

print('\n')

print('Input params:')
print('Rp/Rs: ', rprs)
print('a/Rs: ', a_rs)
print('i (deg): ', inc)

print('\n')

print('Fit params:')
print('Rp/Rs: ', rprs_f)
print('a/Rs: ', a_f)
print('i (deg): ', i_f)<jupyter_output>Stellar mass (Msun):  0.17438311970562165 Stellar radius (Rsun):  0.19884797856314


Input params:
Rp/Rs:  0.036066
a/Rs:  33.79155851141583
i (deg):  89.98


Fit params:
Rp/Rs:  0.03629227534998896
a/Rs:  33.87705752650823
i (deg):  90.06616024528012<jupyter_text>Determining T14 and T23<jupyter_code>T14dist = ph.get_T14(pdist, rdist, adist, idist)
T14errs = ph.get_sigmas(T14dist)

T23dist = ph.get_T23(pdist, rdist, adist, idist)
T23errs = ph.get_sigmas(T23dist)<jupyter_output><empty_output><jupyter_text>Get $g$<jupyter_code>gs, rho_c = ph.get_g_distribution(rho_star, pdist, rdist, T14dist, T23dist)

g_mean = ph.mode(gs)
g_sigma = np.mean(np.abs(ph.get_sigmas(gs)))
g_mean
g_sigma
g = g_mean
gerr = g_sigma<jupyter_output><empty_output><jupyter_text>Dynesty (e,w) Fitting LIkelihood Functions<jupyter_code>ndim = 2

def loglike(theta):
    """The log-likelihood function."""

    w, e = theta
    
    model = (1+e*np.sin(w*(np.pi/180.)))/np.sqrt(1-e**2)
    sigma2 = gerr ** 2
    
    return -0.5 * np.sum((g - model) ** 2 / sigma2 + np.log(sigma2))

def unife_prior_transform(utheta):
    """Uniform eccentricity prior"""
    
    uw, ue = utheta
    w = 360.*uw-90.
    e = 1. * ue
    
    return w, e


def betae_prior_transform(utheta):
    """Uniform eccentricity prior"""
    
    uw, ue = utheta
    w = 360.*uw-90.
    
    a, b = 0.867, 3.03
    e = scipy.stats.beta.ppf(ue, a, b)
    
    return w, e<jupyter_output><empty_output><jupyter_text>Uniform Prior<jupyter_code>dsampler = dynesty.DynamicNestedSampler(loglike, unife_prior_transform, ndim=2, bound='multi', sample='rstagger')
dsampler.run_nested()
dres = dsampler.results
truths = [w, e]

fig, axes = dyplot.cornerplot(dres, truths=truths, show_titles=True, title_kwargs={'y': 1.04}, labels=["w", "e"],
                              fig=plt.subplots(2, 2, figsize=(8, 8)))<jupyter_output><empty_output><jupyter_text>Beta-distribution Prior<jupyter_code>dsampler = dynesty.DynamicNestedSampler(loglike, betae_prior_transform, ndim=2, bound='multi', sample='rstagger')
dsampler.run_nested()
dres = dsampler.results
truths = [w, e]
fig, axes = dyplot.cornerplot(dres, truths=truths, show_titles=True, title_kwargs={'y': 1.04}, labels=["w", "e"],
                              fig=plt.subplots(2, 2, figsize=(8, 8)))<jupyter_output><empty_output><jupyter_text>--------------------------------------------------------------------------------------  $e = 0.9$, $\omega = -85.0$  --------------------------------------------------------------------------------------<jupyter_code># 30 minute cadence
cadence = 0.02142857142857143

time = np.arange(-25, 25, cadence)
# Calculate flux from transit model
e = 0.9
w = -85.0
flux = ph.integratedlc(time, period, rprs, a_rs, e, i, w)

# Adding some gaussian noise
noise = np.random.normal(0,0.00006,len(time))
nflux = flux+noise

flux_err = np.array([0.00006]*len(nflux))<jupyter_output><empty_output><jupyter_text>Fitting the transit<jupyter_code>periodPDF = ph.get_period_dist(time, nflux, 4, 6, arrlen)
print('Period fit: ', ph.mode(periodPDF))
pdist = periodPDF
per_guess = ph.mode(pdist)
ttimes = np.concatenate((-np.arange(0, time[-1], period)[1:], np.arange(0, time[-1], period)))
ttimes = np.sort(ttimes)

time, flux, flux_err = ph.get_transit_cutout_full(ttimes, 4, time, nflux, flux_err)

t0 = 0.0<jupyter_output><empty_output><jupyter_text>Dynesty Transit Fit Likelihood Functions<jupyter_code>ndim = 5

def tfit_loglike(theta):
    """
    Transit fit dynesty function

    model = ph.integratedlc_fitter()
    gerr = sigma of g distribution

    """

    per, rprs, a_rs, inc, t0 = theta

    model = ph.integratedlc_fitter(time, per, rprs, a_rs, inc, t0)
    sigma2 = flux_err ** 2

    return -0.5 * np.sum((flux - model) ** 2 / sigma2 + np.log(sigma2))


def tfit_prior_transform(utheta):
    """Transforms samples `u` drawn from the unit cube to samples to those
    from our uniform prior within [-10., 10.) for each variable."""

    uper, urp, ua, uinc, ut0 = utheta

    per = 3.*uper+3.
    rprs = urp
    a_rs = ua*15.+5.
    inc = uinc*3.+87.
    t0 = 2.*ut0-1.

    return per, rprs, a_rs, inc, t0<jupyter_output><empty_output><jupyter_text>Run dynesty transit fit<jupyter_code>dsampler = dynesty.DynamicNestedSampler(tfit_loglike, tfit_prior_transform, ndim=ndim, nlive=1500,
                                        bound='balls', sample='rwalk')
dsampler.run_nested()
dres = dsampler.results
# Corner Plot of Transit Fit
fig, axes = dyplot.cornerplot(dres, labels=["period", "Rp/Rs", "a/Rs", "i", "t0"])
# Thinning distributions to size arrlen

pdist = random.choices(dres.samples[:,0], k=arrlen)
rdist = random.choices(dres.samples[:,1], k=arrlen)
adist = random.choices(dres.samples[:,2], k=arrlen)
idist = random.choices(dres.samples[:,3], k=arrlen)
t0dist = random.choices(dres.samples[:,4], k=arrlen)

per_f = ph.mode(pdist)
rprs_f = ph.mode(rdist)
a_f = ph.mode(adist)
i_f = ph.mode(idist)
t0_f = ph.mode(t0dist)
# Create a light curve with the fit parameters
fit1 = ph.integratedlc_fitter(time1, per_f, rprs_f, a_f, i_f, t0_f)
plt.errorbar(time1, nflux1, yerr=fluxerr1, c='blue', alpha=0.5, label='Original LC')
plt.plot(time1, fit1, c='red', alpha=1.0, label='Fit LC')
plt.xlim(-0.1, 0.1)
plt.legend()

print('Stellar mass (Msun): ', mstar, 'Stellar radius (Rsun): ', rstar)

print('\n')

print('Input params:')
print('Rp/Rs: ', rprs)
print('a/Rs: ', a_rs)
print('i (deg): ', i)

print('\n')

print('Fit params:')
print('Rp/Rs: ', rprs_f)
print('a/Rs: ', a_f)
print('i (deg): ', i_f)<jupyter_output>Stellar mass (Msun):  0.17438311970562165 Stellar radius (Rsun):  0.19884797856314


Input params:
Rp/Rs:  0.036066
a/Rs:  33.79155851141583
i (deg):  89.98


Fit params:
Rp/Rs:  0.03621114696739527
a/Rs:  7.857472745516413
i (deg):  89.06191781175463<jupyter_text>Determining T14 and T23<jupyter_code>T14dist = ph.get_T14(pdist, rdist, adist, idist)
T14errs = ph.get_sigmas(T14dist)

T23dist = ph.get_T23(pdist, rdist, adist, idist)
T23errs = ph.get_sigmas(T23dist)<jupyter_output>/Users/ssagear/Dropbox (UFL)/Research/MetallicityProject/photoeccentric/photoeccentric/photoeccentric.py:113: RuntimeWarning: invalid value encountered in sqrt
  T23[j] = (p[j]/np.pi)*np.arcsin(rs_a[j]*(np.sqrt(((1-rprs[j])**2)-b[j]**2))/np.sin(i[j]*(np.pi/180.0))) #Equation 14 in exoplanet textbook<jupyter_text>Get $g$<jupyter_code>gs, rho_c = ph.get_g_distribution(rho_star, pdist, rdist, T14dist, T23dist)

g_mean = ph.mode(gs)
g_sigma = np.mean(np.abs(ph.get_sigmas(gs)))
g_mean
g_sigma
g = g_mean
gerr = g_sigma<jupyter_output><empty_output><jupyter_text>Dynesty (e,w) Fitting LIkelihood Functions<jupyter_code>ndim = 2

def loglike(theta):
    """The log-likelihood function."""

    w, e = theta
    
    model = (1+e*np.sin(w*(np.pi/180.)))/np.sqrt(1-e**2)
    sigma2 = gerr ** 2
    
    return -0.5 * np.sum((g - model) ** 2 / sigma2 + np.log(sigma2))

def unife_prior_transform(utheta):
    """Uniform eccentricity prior"""
    
    uw, ue = utheta
    w = 360.*uw-90.
    e = 1. * ue
    
    return w, e


def betae_prior_transform(utheta):
    """Uniform eccentricity prior"""
    
    uw, ue = utheta
    w = 360.*uw-90.
    
    a, b = 0.867, 3.03
    e = scipy.stats.beta.ppf(ue, a, b)
    
    return w, e<jupyter_output><empty_output><jupyter_text>Uniform Prior<jupyter_code>dsampler = dynesty.DynamicNestedSampler(loglike, unife_prior_transform, ndim=2, bound='multi', sample='rstagger')
dsampler.run_nested()
dres = dsampler.results
truths = [w, e]

fig, axes = dyplot.cornerplot(dres, truths=truths, show_titles=True, title_kwargs={'y': 1.04}, labels=["w", "e"],
                              fig=plt.subplots(2, 2, figsize=(8, 8)))<jupyter_output><empty_output><jupyter_text>Beta-distribution Prior<jupyter_code>dsampler = dynesty.DynamicNestedSampler(loglike, betae_prior_transform, ndim=2, bound='multi', sample='rstagger')
dsampler.run_nested()
dres = dsampler.results
truths = [w, e]
fig, axes = dyplot.cornerplot(dres, truths=truths, show_titles=True, title_kwargs={'y': 1.04}, labels=["w", "e"],
                              fig=plt.subplots(2, 2, figsize=(8, 8)))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1729903968840239, 'fraction_numerical': 0.04714256933718353, 'mean_word_length': 4.27710843373494, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 81, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '98503', 'license': 'MIT', 'n_tokens_mistral': 6519, 'n_tokens_neox': 5944, 'n_words': 1131, 'path': 'notebooks/.ipynb_checkpoints/dynestynew-checkpoint.ipynb', 'repo_name': 'ssagear/metallicity', 'size': '14891'}