Here is a short readable (and copy-and-pastable) example as to how to use Namaste 2 to fit a single transit.
In [1]:
from namaste import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
#%matplotlib inline
%reload_ext autoreload
%autoreload 2
In [3]:
#Using 203311200 (from Osborn 2016) as a test case
st=namaste.Star('203311200',namaste.Settings(nsteps=1000,nthreads=4,mission='K2'))
In [4]:
#Can add data directly here:
st.addRad(2.024,0.176)#solar
st.addMass(1.33,0.35)#solar
st.addTeff(6088,250)#K
#Caculates density:
st.addDens()
In [9]:
#Or you can use a csv file, for example:
st.csvfile_dat('ExampleData/EPIC203311200_beststelpars.csv')
#st.addLightcurve('/Volumes/CHEOPS/K2_MONOS/BestLCs/EPIC203311200_bestlc_ev.npy')
In [5]:
#Saving the star
st.SaveAll(overwrite=True)
In [7]:
#We need to use Teff and Logg information to initialise Limb Darkening (default: Kepler)
st.initLD()
To change settings use the st.settings approach
In [8]:
st.settings.printall()
In [9]:
#eg change number of steps in mcmc:
st.settings.update(nsteps=8000)
#or the "sigma" at which to cut anomalies:
st.settings.update(anomcut=4.5)
#Not doing a GP!
st.settings.update(GP=True)
st.settings.update(kernel='Real')
In [10]:
st.addLightcurve('ExampleData/EPIC203311200_bestlc_ev.npy')
Out[10]:
In [11]:
plt.plot(st.Lcurve.lc[:,0],st.Lcurve.lc[:,1],'.')
plt.plot(st.Lcurve.lc[st.Lcurve.lcmask,0],st.Lcurve.lc[st.Lcurve.lcmask,1],',')
Out[11]:
In [12]:
# Can set 'kernel' in settings as "quasi" or "Real"
st.settings.update(kernel='Real')
st.addGP()
#This then optimizes the hyperparameters using the lightcurve (and the out-of-transit regions, if a transit has been specified)
st.Optimize_GP()
So we have set up our star, importantly with the density included. Now we must add a monotransit. This is set up such that multiple monotransits can be added (although I have yet to check this feature).
The format is $t_{\rm cen}$, $b$, and $R_p/R_s$
In [13]:
#Adds this to mono1
st.AddMonotransit(2121.02,0.5,0.0032,name='mono1')
In [14]:
#Accessing pure transit fit:
st.BuildMeanModel()
model=st.meanmodel_comb.get_value(np.arange(2120.6,2121.4,0.005))
plt.plot(st.Lcurve.lc[abs(st.Lcurve.lc[:,0]-st.mono1.tcen)<1.5,0],st.Lcurve.lc[abs(st.Lcurve.lc[:,0]-st.mono1.tcen)<1.5,1],'.')
plt.plot(np.arange(2120.6,2121.4,0.005),model,'--')
Out[14]:
In [18]:
#Now we have our mean-functions (eg transit models) added, we need to build the priors and the model:
st.BuildMeanPriors()
st.BuildMeanModel()
In [15]:
# Then we need the combine the kernel we've optimised with the monotransit mean model into a Celerite GP:
import celerite
st.gp=celerite.GP(kernel=st.kern,mean=st.meanmodel_comb,fit_mean=True)
In [16]:
# Now we need to build all the priors (not just the mean functions)
st.BuildAllPriors(st.gp.get_parameter_names())
In [17]:
# Making walkers initial positions
chx=np.random.choice(st.settings.npdf,st.settings.nwalkers,replace=False)
#Tidying up so that the Celerite priors match our prior names
st.fitdict.update({'mean:'+nm:getattr(st.mono1,nm+'s')[chx] for nm in ['tcen','b','vel','RpRs']})
st.fitdict.update({'mean:'+nm:getattr(st,nm+'s')[chx] for nm in ['LD1','LD2']})
In [18]:
#Removing medians:
for row in st.fitdict:
st.fitdict[row][np.isnan(st.fitdict[row])]=np.nanmedian(np.isnan(st.fitdict[row]))
dists=[st.fitdict[cname] for cname in st.gp.get_parameter_names()]
st.init_mcmc_params=np.column_stack(dists)
In [19]:
#Masking an arbitrary region around the transit for analysis (2.75d in this case)
mask=abs(st.Lcurve.lc[:,0]-st.gp.get_parameter('mean:tcen'))<2.75
In [21]:
# Doing an initial fit with PlotModel:
_=namaste.PlotModel(st.Lcurve.lc[mask,:], st.gp, np.median(st.init_mcmc_params,axis=0), fname=st.settings.outfilesloc+st.objname+'_initfit.png',GP=True)
In [34]:
#RUNNING THE MCMC WITH EMCEE:
import emcee
st.sampler = emcee.EnsembleSampler(st.settings.nwalkers, len(st.gp.get_parameter_vector()), namaste.MonoLogProb, args=(st.Lcurve.lc,st.priors,st.gp), threads=st.settings.nthreads)
If in the following line you see:
PicklingError: Can't pickle <class 'namaste.namaste.MonotransitModel'>: it's not the same object as namaste.namaste.MonotransitModel
then try re-running this script!
In [35]:
#Just checking if the MCMC works/progresses:
st.sampler.run_mcmc(st.init_mcmc_params, 1, rstate0=np.random.get_state())
Out[35]:
Running the MCMC
In [36]:
_ = st.sampler.run_mcmc(st.init_mcmc_params, int(st.settings.nsteps/12), rstate0=np.random.get_state())
Out[36]:
In [38]:
# let's make an area of the samples and their names
st.samples = st.sampler.chain.reshape((-1, len(dists)))
st.samples = np.column_stack((st.samples,st.sampler.lnprobability.reshape(-1)))
st.sampleheaders=list(st.gp.get_parameter_names())+['logprob']
In [35]:
st.RunMCMC()
Now let's save the MCMC samples and plot the MCMC result:
In [39]:
st.SaveMCMC()
In [45]:
st.PlotMCMC()
In [77]:
st_nogp=namaste.Star('203311200',namaste.Settings(nsteps=1000,nthreads=4,mission='K2'))
st_nogp.settings.update(GP=False)
st_nogp.settings.update(verbose=False)
In [78]:
#Or you can use a csv file, for example:
st_nogp.csvfile_dat('ExampleData/EPIC203311200_beststelpars.csv')
#Initialising LDs:
st_nogp.initLD()
In [79]:
#Adding a lightcurve:
st_nogp.addLightcurve('/Volumes/CHEOPS/K2_MONOS/BestLCs/EPIC203311200_bestlc_ev.npy')
Out[79]:
In [80]:
#When we're not doing GPs, lets flatten the LC:
st_nogp.Lcurve.lc=st_nogp.Lcurve.flatten()
In [81]:
plt.plot(st_nogp.Lcurve.lc[:,0],st_nogp.Lcurve.lc[:,1],'.')
plt.plot(st_nogp.Lcurve.lc[st_nogp.Lcurve.lcmask,0],st_nogp.Lcurve.lc[st_nogp.Lcurve.lcmask,1],',')
Out[81]:
In [88]:
#Adds this to mono1
st_nogp.AddMonotransit(2121.02,0.5,0.0032,name='mono1')
#Now we have our mean-functions (eg transit models) added, we need to build the priors and the model:
st_nogp.BuildMeanPriors()
In [90]:
#Accessing pure transit fit:
st_nogp.BuildMeanModel()
model=st_nogp.meanmodel_comb.get_value(np.arange(2120.6,2121.4,0.005))
plt.plot(st_nogp.Lcurve.lc[abs(st_nogp.Lcurve.lc[:,0]-st_nogp.mono1.tcen)<1.5,0],st_nogp.Lcurve.lc[abs(st_nogp.Lcurve.lc[:,0]-st_nogp.mono1.tcen)<1.5,1],'.')
plt.plot(np.arange(2120.6,2121.4,0.005),model,'--')
Out[90]:
In [91]:
st_nogp.RunMCMC()
In [94]:
st_nogp.SaveMCMC(suffix='NoGP',overwrite=True)
In [99]:
st_nogp.PlotMCMC(suffix='NoGP')