Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
if rv_dictionary['GPType'] == 'CeleriteQPKernel':
# Now that we know the type of GP, we extract the "instrument" corresponding to each GP
# parameter. For example, for the ExpSineSquaredKernel, it might happen the user wants
# to have a common GP_Prot parameter shared along many instruments, e.g., GP_Prot_TESS_K2_RV,
# which means the user wants a common Prot for TESS and K2 photometry, and also for the RVs. However,
# the same GP might have a different Gamma, i.e., there might be a GP_Gamma_TESS, GP_Gamma_K2 and GP_Gamma_RV.
# The idea here is to, e.g., in the case of TESS photometry, gather lc_dictionary[instrument]['GP_Prot'] =
# 'TESS_K2_RV', and lc_dictionary[instrument]['GP_Gamma'] = 'TESS':
for GPvariable in ['B','C','L','Prot']:
for pnames in priors.keys():
vec = pnames.split('_')
if (vec[0] == 'GP') and (GPvariable in vec[1]) and ('rv' in vec[-1].lower()):
rv_dictionary['GP_'+GPvariable] = '_'.join(vec[2:])
#for instrument in inames_rv:
rot_kernel = terms.TermSum(RotationTerm(log_amp=np.log(10.),\
log_timescale=np.log(10.0),\
log_period=np.log(3.0),\
log_factor=np.log(1.0)))
# Jitter term; dont add it, jitters will be added directly on the log-like (see Espinoza+2018).
#kernel_jitter = terms.JitterTerm(np.log(100*1e-6))
# Wrap GP object to compute likelihood:
kernel = rot_kernel #+ kernel_jitter
rv_dictionary['GPObject'] = celerite.GP(kernel, mean=0.0)
# Note order of GP Vector: logB, logL, logProt, logC, logJitter
rv_dictionary['GPVector'] = np.zeros(5)
rv_dictionary['GPObject'].compute(rv_dictionary['X'],yerr=rverr_rv)
if rv_dictionary['GPType'] == 'ExpSineSquaredSEKernel':
for GPvariable in ['sigma','alpha','Gamma','Prot']:
for pnames in priors.keys():
lc_dictionary[instrument]['GPType'] = 'CeleriteQPKernel'
break
print('\t Detected ',lc_dictionary[instrument]['GPType'],'for the GP')
# For each instrument for which there are external parameters, activate GP:
lc_dictionary[instrument]['GPDetrend'] = True
# Save variables, standarize them if GP is SEKernel:
lc_dictionary[instrument]['X'] = GPDict[instrument]['variables']
lc_dictionary[instrument]['nX'] = lc_dictionary[instrument]['X'].shape[1]
print('\t (',lc_dictionary[instrument]['nX'],'external parameters)')
if lc_dictionary[instrument]['GPType'] == 'SEKernel':
for i in range(lc_dictionary[instrument]['nX']):
lc_dictionary[instrument]['X'][:,i] = (lc_dictionary[instrument]['X'][:,i] - \
np.mean(lc_dictionary[instrument]['X'][:,i]))/\
np.sqrt(np.var(lc_dictionary[instrument]['X'][:,i]))
if lc_dictionary[instrument]['GPType'] == 'CeleriteQPKernel':
rot_kernel = terms.TermSum(RotationTerm(log_amp=np.log(10.),\
log_timescale=np.log(10.0),\
log_period=np.log(3.0),\
log_factor=np.log(1.0)))
# Now that we know the type of GP, we extract the "instrument" corresponding to each GP
# parameter. For example, for the ExpSineSquaredKernel, it might happen the user wants
# to have a common GP_Prot parameter shared along many instruments, e.g., GP_Prot_TESS_K2_RV,
# which means the user wants a common Prot for TESS and K2 photometry, and also for the RVs. However,
# the same GP might have a different Gamma, i.e., there might be a GP_Gamma_TESS, GP_Gamma_K2 and GP_Gamma_RV.
# The idea here is to, e.g., in the case of TESS photometry, gather lc_dictionary[instrument]['GP_Prot'] =
# 'TESS_K2_RV', and lc_dictionary[instrument]['GP_Gamma'] = 'TESS':
for GPvariable in ['B','C','L','Prot']:
for pnames in priors.keys():
vec = pnames.split('_')
if (vec[0] == 'GP') and (GPvariable in vec[1]) and (instrument in vec):
lc_dictionary[instrument]['GP_'+GPvariable] = '_'.join(vec[2:])
# actual kernel:
self.variables = self.all_kernel_variables[self.kernel_name]
phantomvariable = 0
if self.kernel_name == 'SEKernel':
# Generate GPExpSquared base kernel:
self.kernel = 1.*george.kernels.ExpSquaredKernel(np.ones(self.nX),ndim = self.nX, axes = range(self.nX))
# (Note no jitter kernel is given, as with george one defines this in the george.GP call):
elif self.kernel_name == 'ExpSineSquaredSEKernel':
# Generate the kernels:
K1 = 1.*george.kernels.ExpSquaredKernel(metric = 1.0)
K2 = george.kernels.ExpSine2Kernel(gamma=1.0,log_period=1.0)
self.kernel = K1*K2
# (Note no jitter kernel is given, as with george one defines this in the george.GP call):
elif self.kernel_name == 'CeleriteQPKernel':
# Generate rotational kernel:
rot_kernel = terms.TermSum(RotationTerm(log_amp=np.log(10.),\
log_timescale=np.log(10.0),\
log_period=np.log(3.0),\
log_factor=np.log(1.0)))
# Jitter term:
kernel_jitter = terms.JitterTerm(np.log(100*1e-6))
# Wrap GP kernel and object:
if self.instrument in ['rv','lc']:
self.kernel = rot_kernel
else:
self.kernel = rot_kernel + kernel_jitter
# We are using celerite:
self.use_celerite = True
elif self.kernel_name == 'CeleriteExpKernel':
# Generate exponential kernel:
exp_kernel = terms.RealTerm(log_a=np.log(10.), log_c=np.log(10.))
# Jitter term:
def __add__(self, b):
return TermSum(self, b)
def __init__(self, *terms):
models = []
for term in terms:
models += term.terms
super(TermSum, self).__init__([("terms[{0}]".format(i), t)
for i, t in enumerate(models)])
def __radd__(self, b):
return TermSum(b, self)