GP API

[petigura]

GP example API using Kepler-78 as an Example

# Load in data
data_hires = load_data()
data_harps = load_data()

# RV orbital parameters
basis_input  = 'per tc e w k'
basis_fit = 'per tc secosw sesinw logk' 
params = radvel.Parameters(1,basis=basis_input)

# Initialize values
params['per1'] = radvel.Parameter(value=1201.1 + 0.4)
params['tc1'] = radvel.Parameter(value=2456778 + 1)
params['e1'] = radvel.Parameter(value=0.01)
params['w1'] = radvel.Parameter(value=0.01)
params['logk1'] = radvel.Parameter(value=1)
params['dvdt'] = radvel.Parameter(value=0)
params['curv'] = radvel.Parameter(value=0)

# params = radvel.basis.covert(params, basis_fit) 

params = params.basis.to_any_basis(basis_fit)
params['dvdt'].vary = False
params['curv'].vary = False

# Add in noise parameters
params['gamma_hires'] = radvel.Parameter(value=1.0)
params['gamma_harps'] = radvel.Parameter(value=1.0)

# Define GP parameters 
hyperparam_names = [
    'gp_amp_hires', # GP variability amplitude for HIRES
    'gp_per_share', # GP variability period (shared between all instruments)
    'gp_decay_share' # GP variability decay period  (shared between all instruments)
    gp_jit_hires', # GP white noise component (jitter) for HIRES 
]
params['gp_amp_hires'] = radvel.Parameter(value=1.0) # period
params['gp_per_share'] = radvel.Parameter(value=1.0) # 
params['gp_decay_share'] = radvel.Parameter(value=1.0) # 
params['gp_jit_hires'] = radvel.Parameter(value=1.0) # 

# GP should probably be implemented as their own object for modularity and to implement different compute backends
gp = QuasiPeriodicGP(names, backend='cholesky') # other backends could be george, etch
gp.mean(params, x, y, yerr) # compute mean function given best fit parameters
gp.std(params, x, y, yerr) # compute 1 sigma bands
gp.loglike(params, x, y, yerr) # compute loglike

# Define HIRES likelihood 
like_hires = GPLikelihood(data_hires['t'], data_hires['rv'], data_hires['rv_err'], model, gp)
like_hires.loglike(params) # compute loglikelihood
like_hires.mean(params) # compute mean function
like_hires.mean(params) # compute 1sig bands

# Define GP parameters for HARPS
hyperparam_names = [
    'gp_amp_harps', # GP variability amplitude for HARPS
    'gp_per_share', # GP variability period (shared between all instruments)
    'gp_decay_share', # GP variability decay period  (shared between all instruments)
    'gp_jit_harps', # GP white noise component (jitter) for HARPS
]
params['gp_amp_harps'] = radvel.Parameter(value=1.0) # period
params['gp_per_share'] = radvel.Parameter(value=1.0) # 
params['gp_decay_share'] = radvel.Parameter(value=1.0) # 
params['gp_jit_harps'] = radvel.Parameter(value=1.0) # 
gp = QuasiPeriodicGP(names, backend='cholesky')

like_harps = GPLikelihood(data_harps['t'], data_harps['rv'], data_harps['rv_err'], model, gp)

# Define composite likelihood
like = radvel.likelihood.CompositeLikelihood([like_hires,like_harps])

# Define posterior (likelihood+priors)
post = radvel.posterior.Posterior(like)
post0 = copy.deepcopy(post)
post = radvel.posterior.Posterior(like)
post.priors += [radvel.prior.EccentricityPrior( 1 )] # Keeps eccentricity < 1

[petigura]

@sblunt, I made some tweaks to our draft API to make it applicable to the Kepler 78 dataset

[petigura]

@sblunt, small tweaks to the draft API per our discussion. The GP object might more appropriately be a "kernel" object

# Load in data
data_hires = load_data()
data_harps = load_data()

# RV orbital parameters
basis_input  = 'per tc e w k'
basis_fit = 'per tc secosw sesinw logk' 
params = radvel.Parameters(1,basis=basis_input)

# Initialize values
params['per1'] = radvel.Parameter(value=1201.1 + 0.4)
params['tc1'] = radvel.Parameter(value=2456778 + 1,vary=False)
params['e1'] = radvel.Parameter(value=0.01,vary=False)
params['w1'] = radvel.Parameter(value=0.01,vary=False)
params['logk1'] = radvel.Parameter(value=1,vary=True)
params['dvdt'] = radvel.Parameter(value=0, vary=True)
params['curv'] = radvel.Parameter(value=0,vary=False)

params = params.basis.to_any_basis(basis_fit)

# Add in noise parameters
params['gamma_hires'] = radvel.Parameter(value=1.0,vary=True)
params['gamma_harps'] = radvel.Parameter(value=1.0,vary=True)

# Define GP parameters for HIRES
hyperparam_names = [
    'gp_amp_hires', # GP variability amplitude for HIRES
    'gp_per_share', # GP variability period (shared between all instruments)
    'gp_decay_share', # GP variability decay period  (shared between all instruments)
    'gp_jit_hires', # GP white noise component (jitter) for HIRES 
]
params['gp_amp_hires'] = radvel.Parameter(value=1.0) # period
params['gp_per_share'] = radvel.Parameter(value=1.0) # 
params['gp_decay_share'] = radvel.Parameter(value=1.0) # 
params['gp_jit_hires'] = radvel.Parameter(value=1.0) # 

# GP should probably be implemented as their own object for modularity and to implement different compute backends
gp = QuasiPeriodicGP(hyperparam_names, backend='cholesky') # other backends could be george, etc

# Define HIRES likelihood 
like_hires = GPLikelihood(data_hires['t'], data_hires['rv'], data_hires['rv_err'], model, gp)
like_hires.loglike(params) # compute loglikelihood
like_hires.mean(params) # compute mean function
like_hires.std(params) # compute 1sig bands

# Define GP parameters for HARPS
hyperparam_names = [
    'gp_amp_harps', # GP variability amplitude for HARPS
    'gp_per_share', # GP variability period (shared between all instruments)
    'gp_decay_share', # GP variability decay period  (shared between all instruments)
    'gp_jit_harps', # GP white noise component (jitter) for HARPS
]

# Add in the parameters needed for HARPS (no need to re-initialize the shared parameters)
params['gp_amp_harps'] = radvel.Parameter(value=1.0) # period
params['gp_jit_harps'] = radvel.Parameter(value=1.0) # 
gp = QuasiPeriodicGP(hyperparam_names, backend='cholesky')

like_harps = GPLikelihood(data_harps['t'], data_harps['rv'], data_harps['rv_err'], model, gp)

# Define composite likelihood
like = radvel.likelihood.CompositeLikelihood([like_hires,like_harps])

# Define posterior (likelihood+priors)
post = radvel.posterior.Posterior(like)
post0 = copy.deepcopy(post)
post = radvel.posterior.Posterior(like)
post.priors += [radvel.prior.EccentricityPrior( 1 )] # Keeps eccentricity < 1

[bjfultn]

@petigura This doesn't seem like an issue with a well-defined resolution and closing criteria but more like notes. Not sure exactly where to store this. Maybe @sblunt should transfer this information over to Confluence?

[sblunt]

I can do that. We could make a "Notes on RadVel Meetings" page or something like that.

[bjfultn]

OK, or you could try to set up a "project" here on GitHub. Whichever you prefer.

[sblunt]

I moved all of this information to Confluence. Under "projects" on the main page, click on "RadVel" to access these notes.