Skip to content

Commit

Permalink
Fixed code from crashing for CeleriteMaternExpKernel. It did because …
Browse files Browse the repository at this point in the history 
…amplitude of matern is set to 1.
  • Loading branch information
@nespinoza
nespinoza committed Sep 4, 2019
1 parent 3b7f2d8 commit 85ae88b
Showing 1 changed file with 74 additions and 33 deletions.
107 changes: 74 additions & 33 deletions juliet/fit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -270,13 +270,17 @@ def save_regressors(self,fname, GP_arguments, global_model):
fout = open(fname,'w')
for GP_instrument in GP_arguments.keys():
GP_regressors = GP_arguments[GP_instrument]
for i in range(GP_regressors.shape[0]):
for j in range(GP_regressors.shape[1]):
fout.write('{0:.10f} '.format(GP_regressors[i,j]))
if GP_instrument is not 'global_model':
fout.write('{0:}\n'.format(GP_instrument))
else:
fout.write('\n')
multi_dimensional = False
if len(GP_regressors.shape) == 2:
multi_dimensional = True
if multi_dimensional:
for i in range(GP_regressors.shape[0]):
for j in range(GP_regressors.shape[1]):
fout.write('{0:.10f} '.format(GP_regressors[i,j]))
fout.write('{0:}\n'.format(GP_instrument))
else:
for i in range(GP_regressors.shape[0]):
fout.write('{0:.10f} {1:}\n'.format(GP_regressors[i],GP_instrument))
fout.close()

def save_data(self,fname,t,y,yerr,instruments,lm_boolean,lm_arguments):
Expand Down Expand Up @@ -1268,19 +1272,25 @@ def evaluate_lc_model(self, instrument, parameter_values = None, resampling = No
This functions evaluate the current lc model into a new set of times.
"""
if resampling is not None:
self.model[instrument]['params'], self.model[instrument]['m'] = self.init_batman(self.dictionary[instrument]['ldlaw'],\
nresampling = nresampling,\
etresampling = etresampling)
if resampling:
self.model[instrument]['params'], self.model[instrument]['m'] = self.init_batman(self.dictionary[instrument]['ldlaw'],\
nresampling = nresampling,\
etresampling = etresampling)
else:
self.model[instrument]['params'], self.model[instrument]['m'] = self.init_batman(self.dictionary[instrument]['ldlaw'])
# Check if user gave input parameter_values dictionary. If that's the case, generate again the
# full lightcurve model:
if parameter_values is not None:
self.generate_lc_model(parameter_values)

if resampling is not None:
# get lc, return, then turn all back to normal:
self.model[instrument]['params'], self.model[instrument]['m'] = self.init_batman(self.dictionary[instrument]['ldlaw'],\
nresampling = self.dictionary[instrument]['nresampling'],\
etresampling = self.dictionary[instrument]['exptimeresampling'])
if self.dictionary[instrument]['resampling']:
self.model[instrument]['params'], self.model[instrument]['m'] = self.init_batman(self.dictionary[instrument]['ldlaw'],\
nresampling = self.dictionary[instrument]['nresampling'],\
etresampling = self.dictionary[instrument]['exptimeresampling'])
else:
self.model[instrument]['params'], self.model[instrument]['m'] = self.init_batman(self.dictionary[instrument]['ldlaw'])
return self.model[instrument]['full']

def predict_lc_model(self, t, posteriors = None):
Expand Down Expand Up @@ -1519,6 +1529,17 @@ def __init__(self, data, modeltype, pl = 0.0, pu = 1.0, ecclim = 1., ta = 245846
# Extract number of linear model terms per instrument:
if self.lm_boolean[instrument]:
self.lm_n[instrument] = self.lm_arguments[instrument].shape[1]
# An array of ones to copy around:
self.model[instrument]['ones'] = np.ones(len(self.instrument_indexes[instrument]))
# Generate internal model variables of interest to the user. First, the lightcurve model in the notation of juliet (Mi)
# (full lightcurve plus dilution factors and mflux):
self.model[instrument]['M'] = np.ones(len(self.instrument_indexes[instrument]))
# Linear model (in the notation of juliet, LM):
self.model[instrument]['LM'] = np.zeros(len(self.instrument_indexes[instrument]))
# Now, generate dictionary that will save the final full model (M + LM):
self.model[instrument]['full'] = np.zeros(len(self.instrument_indexes[instrument]))
# Same for the errors:
self.model[instrument]['full_errors'] = np.zeros(len(self.instrument_indexes[instrument]))
if self.dictionary[instrument]['TransitFit']:
# First, take the opportunity to initialize transit lightcurves for each instrument:
if self.dictionary[instrument]['resampling']:
Expand All @@ -1527,20 +1548,9 @@ def __init__(self, data, modeltype, pl = 0.0, pu = 1.0, ecclim = 1., ta = 245846
etresampling = self.dictionary[instrument]['exptimeresampling'])
else:
self.model[instrument]['params'], self.model[instrument]['m'] = self.init_batman(self.dictionary[instrument]['ldlaw'])
# Generate internal model variables of interest to the user. First, the lightcurve model in the notation of juliet (Mi)
# (full lightcurve plus dilution factors and mflux):
self.model[instrument]['M'] = np.ones(len(self.instrument_indexes[instrument]))
# Linear model (in the notation of juliet, LM):
self.model[instrument]['LM'] = np.zeros(len(self.instrument_indexes[instrument]))
# Now, generate dictionary that will save the final full model (M + LM):
self.model[instrument]['full'] = np.zeros(len(self.instrument_indexes[instrument]))
# Same for the errors:
self.model[instrument]['full_errors'] = np.zeros(len(self.instrument_indexes[instrument]))
# Individual transit lightcurves for each planet:
for i in self.numbering:
self.model[instrument]['p'+str(i)] = np.ones(len(self.instrument_indexes[instrument]))
# An array of ones to copy around:
self.model[instrument]['ones'] = np.ones(len(self.t[self.instrument_indexes[instrument]]))
# Now proceed with instrument namings:
for pname in self.priors.keys():
# Check if variable name is a limb-darkening coefficient:
Expand All @@ -1560,6 +1570,17 @@ def __init__(self, data, modeltype, pl = 0.0, pu = 1.0, ecclim = 1., ta = 245846
self.mdilution_iname[instrument] = '_'.join(vec[1:])
else:
self.mdilution_iname[instrument] = vec[1]
else:
# Now proceed with instrument namings:
for pname in self.priors.keys():
# Check if it is a dilution factor:
if pname[0:9] == 'mdilution':
vec = pname.split('_')
if len(vec)>2:
if instrument in vec:
self.mdilution_iname[instrument] = '_'.join(vec[1:])
else:
self.mdilution_iname[instrument] = vec[1]
# Set the model-type to M(t):
self.evaluate = self.evaluate_lc_model
self.predict = self.predict_lc_model
Expand Down Expand Up @@ -1747,9 +1768,9 @@ def set_parameter_vector(self, parameter_values):
elif self.kernel_name == 'CeleriteMaternExpKernel':
self.parameter_vector[0] = np.log(parameter_values['GP_sigma_'+self.input_instrument[0]])
self.parameter_vector[1] = np.log(parameter_values['GP_timescale_'+self.input_instrument[1]])
self.parameter_vector[2] = np.log(parameter_values['GP_rho_'+self.input_instrument[2]])
self.parameter_vector[3] = np.log(parameter_values['GP_rho_'+self.input_instrument[2]])
if not self.global_GP:
self.parameter_vector[3] = np.log(parameter_values['sigma_w_'+self.instrument]*self.sigma_factor)
self.parameter_vector[4] = np.log(parameter_values['sigma_w_'+self.instrument]*self.sigma_factor)
elif self.kernel_name == 'CeleriteSHOKernel':
self.parameter_vector[0] = np.log(parameter_values['GP_S0_'+self.input_instrument[0]])
self.parameter_vector[1] = np.log(parameter_values['GP_Q_'+self.input_instrument[1]])
Expand Down Expand Up @@ -1848,6 +1869,7 @@ def __init__(self, data, model_type, instrument, george_hodlr = True):
# Initialize each kernel on the GP object. First, set the variables to the ones defined above. Then initialize the
# 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))
Expand All @@ -1867,7 +1889,10 @@ def __init__(self, data, model_type, instrument, george_hodlr = True):
# Jitter term:
kernel_jitter = terms.JitterTerm(np.log(100*1e-6))
# Wrap GP kernel and object:
self.kernel = rot_kernel + kernel_jitter
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':
Expand All @@ -1876,7 +1901,10 @@ def __init__(self, data, model_type, instrument, george_hodlr = True):
# Jitter term:
kernel_jitter = terms.JitterTerm(np.log(100*1e-6))
# Wrap GP kernel and object:
self.kernel = exp_kernel + kernel_jitter
if self.instrument in ['rv','lc']:
self.kernel = exp_kernel
else:
self.kernel = exp_kernel + kernel_jitter
# We are using celerite:
self.use_celerite = True
elif self.kernel_name == 'CeleriteMaternKernel':
Expand All @@ -1885,7 +1913,10 @@ def __init__(self, data, model_type, instrument, george_hodlr = True):
# Jitter term:
kernel_jitter = terms.JitterTerm(np.log(100*1e-6))
# Wrap GP kernel and object:
self.kernel = matern_kernel + kernel_jitter
if self.instrument in ['rv','lc']:
self.kernel = matern_kernel
else:
self.kernel = matern_kernel + kernel_jitter
# We are using celerite:
self.use_celerite = True
elif self.kernel_name == 'CeleriteMaternExpKernel':
Expand All @@ -1895,7 +1926,14 @@ def __init__(self, data, model_type, instrument, george_hodlr = True):
# Jitter term:
kernel_jitter = terms.JitterTerm(np.log(100*1e-6))
# Wrap GP kernel and object:
self.kernel = exp_kernel*matern_kernel + kernel_jitter
if self.instrument in ['rv','lc']:
self.kernel = exp_kernel*matern_kernel
else:
self.kernel = exp_kernel*matern_kernel + kernel_jitter
# We add a phantom variable because we want to leave index 2 without value ON PURPOSE: the idea is
# that here, that is always 0 (because this defines the log(sigma) of the matern kernel in the
# multiplication, which we set to 1).
phantomvariable = 1
# We are using celerite:
self.use_celerite = True
elif self.kernel_name == 'CeleriteSHOKernel':
Expand All @@ -1904,7 +1942,10 @@ def __init__(self, data, model_type, instrument, george_hodlr = True):
# Jitter term:
kernel_jitter = terms.JitterTerm(np.log(100*1e-6))
# Wrap GP kernel and object:
self.kernel = sho_kernel + kernel_jitter
if self.instrument in ['rv','lc']:
self.kernel = sho_kernel
else:
self.kernel = sho_kernel + kernel_jitter
# We are using celerite:
self.use_celerite = True
# Check if use_celerite is True; if True, check that the regressor is ordered. If not, don't do the self.init_GP():
Expand All @@ -1924,10 +1965,10 @@ def __init__(self, data, model_type, instrument, george_hodlr = True):
# (e.g., global photometric signal, or global RV signal) that assumes a given
# GP realization for all instruments (but allows different jitters for each
# instrument, added in quadrature to the self.yerr):
self.parameter_vector = np.zeros(len(self.variables))
self.parameter_vector = np.zeros(len(self.variables)+panthomvariable)
self.global_GP = True
else:
# If GP per instrument, then there is one jitter term per instrument directly added in the model:
self.parameter_vector = np.zeros(len(self.variables)+1)
self.parameter_vector = np.zeros(len(self.variables)+1+phantomvariable)
self.global_GP = False
self.set_input_instrument(data.priors)

0 comments on commit 85ae88b

Please sign in to comment.