refactor likelihoods

[petigura]

Would be more clear if all parameters were defined at the top of a setup file, as opposed to at several places as the likelihoods are being setup. See comments below.

data = pd.read_csv(os.path.join(radvel.DATADIR,'k2-131.txt'), sep=' ')
t = np.array(data['time']) 
vel = np.array(data['mnvel'])
errvel = np.array(data['errvel'])
telgrps = data.groupby('tel').groups
bjd = 0. 

# Set initial parameter value guesses.
fitting_basis = 'per tc secosw sesinw k'
params = radvel.Parameters(nplanets,basis=fitting_basis, planet_letters=planet_letters)
params['per1'] = radvel.Parameter(value=0.3693038, vary=False) # constrained from transit
params['tc1'] = radvel.Parameter(value=2457582.936, vary=False) # constrained from transit 
params['sesinw1'] = radvel.Parameter(value=0.,vary=False) 
params['secosw1'] = radvel.Parameter(value=0.,vary=False)
params['k1'] = radvel.Parameter(value=6.55)
params['dvdt'] = radvel.Parameter(value=0.,vary=False)
params['curv'] = radvel.Parameter(value=0.,vary=False)
params['gp_B_pfs'] = radvel.Parameter(value=30**2., vary=True)
params['gp_B_harps-n'] = radvel.Parameter(value=1., vary=True) 
params['gp_C'] = radvel.Parameter(value=1., vary=True) 
params['gp_L'] = radvel.Parameter(value=9., vary=True)
params['gp_Prot'] = radvel.Parameter(value=9., vary=True)

time_base = np.median(t)    
model = radvel.model.RVModel(params,time_base=time_base)

# Define PFS likelihood
hnames = ['gp_B_pfs','gp_C','gp_L','gp_Prot']
_data = data.query('tel=="pfs"')
like1 = radvel.likelihood.CeleriteLikelihood(model, _data.time, _data.mnvel, _data.errvel,hnames,suffix='_pfs') # no suffix term
like1.params['gamma_pfs'] = radvel.Parameter(value=np.mean(vel[telgrps['pfs']])) # <- define the linking in the likelihood definition
like1.params['jit_pfs'] = radvel.Parameter(value=3) # <- define the linking in the likelihood definition

# Define HARPS-N likelihood 
hnames = ['gp_B_harps-n','gp_C','gp_L','gp_Prot']
_data = data.query('tel=="harps-n"')
like2 = radvel.likelihood.CeleriteLikelihood(model, _data.time, _data.mnvel, _data.errvel, hnames,suffix='_harps-n') # no suffix term
like2.params['gamma_harps-n'] = radvel.Parameter(value=np.mean(vel[telgrps['harps-n']])) 
like2.params['jit_harps-n'] = radvel.Parameter(value=3) # <- define the linking in the likelihood definition
like = radvel.likelihood.CompositeLikelihood([like1,like2]) # <- define the linking in the likelihood definition

# Define priors  
post = radvel.posterior.Posterior(like) 
priors = [
   radvel.prior.Gaussian('per1', Porb, Porb_unc),
   radvel.prior.Gaussian('tc1', Tc, Tc_unc),
   radvel.prior.Jeffreys('k1', 0.01, 10.),
   radvel.prior.Jeffreys('jit_pfs', 0.01, 10.),
   radvel.prior.Jeffreys('jit_harps-n', 0.01,10.),
   radvel.prior.Gaussian('gp_L',9.5*2,1.), # constraints from photometry on hyperparams (Dai et al. 2017)
   radvel.prior.Gaussian('gp_Prot',9.64,0.12),
   radvel.prior.HardBounds('gp_C',0.,1e10),  # keep other two hyperparams positive
   radvel.prior.HardBounds('gp_B_pfs',0.,1e10),
   radvel.prior.HardBounds('gp_B_harps-n',0.,1e10)
]
post.priors = priors
res = optimize.minimize(
    post.neglogprob_array, post.get_vary_params(), method='Powell',
    options=dict(maxiter=200, maxfev=100000, xatol=1e-8)
)

print(post)

[bjfultn]

@petigura should we wait for this before pulling in #185?

[petigura]

I'd say it's up to you as the benevolent dictator of this project. There is some value in getting out this new capability to the community ASAP as many people are interested in RVs and GPs. The changes I was thinking would potentially break existing capability. Generally, these are reserved for major version bumps. RadVel 2.0?

[bjfultn]

Ah ok, I didn't realize that we were going to break backwards compatibility here. I thought this would just be a change to the documentation and tutorials.

If that is really the case then I don't think we should wait on this to get the Celerite kernel out there.