Skip to content

Commit

Permalink
Merge pull request #76 from tronsgaard/fix-supersampling
Browse files Browse the repository at this point in the history 
Fix missing LC supersampling
  • Loading branch information
@nespinoza
nespinoza authored Mar 16, 2022
2 parents 63d6f17 + 737280e commit 162f197
Showing 1 changed file with 38 additions and 45 deletions.
83 changes: 38 additions & 45 deletions juliet/fit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1282,11 +1282,11 @@ def __init__(self, data, sampler = 'multinest', n_live_points = 500, nwalkers =
self.sampler_prefix = self.sampler+'_'

# Before starting, check if force_dynesty or force_pymultinest is on; change options accordingly:
if force_dynesty and (self.sampler is 'multinest'):
if force_dynesty and (self.sampler == 'multinest'):
print('PyMultinest installation not detected. Forcing dynesty as the sampler.')
self.sampler = 'dynesty'
self.sampler_prefix = '_dynesty_NS_'
if force_pymultinest and (self.sampler is 'dynesty'):
if force_pymultinest and (self.sampler == 'dynesty'):
print('dynesty installation not detected. Forcing PyMultinest as the sampler.')
self.sampler = 'multinest'
self.sampler_prefix = ''
Expand Down Expand Up @@ -1827,8 +1827,8 @@ def get_GP_plus_deterministic_model(self, parameter_values, instrument = None):
else:
return self.model[instrument]['deterministic']

def evaluate_model(self, instrument = None, parameter_values = None, resampling = None, nresampling = None, etresampling = None, \
all_samples = False, nsamples = 1000, return_samples = False, t = None, GPregressors = None, LMregressors = None, \
def evaluate_model(self, instrument = None, parameter_values = None,
all_samples = False, nsamples = 1000, return_samples = False, t = None, GPregressors = None, LMregressors = None,
return_err = False, alpha = 0.68, return_components = False, evaluate_transit = False):
"""
This function evaluates the current lc or rv model given a set of posterior distribution samples and/or parameter values. Example usage:
Expand All @@ -1851,15 +1851,6 @@ def evaluate_model(self, instrument = None, parameter_values = None, resampling
'q1_TESS', etc.), and inside each of those keys an array of N samples is expected (i.e., parameter_values['p_p1'] is an array of length N). The
indexes have to be consistent between different parameters.
:param resampling: (optional, boolean)
Boolean indicating if the model needs to be resampled or not. Only works for lightcurves.
:param nresampling: (optional, int)
Number of points to resample for a given time-stamp. Only used if resampling = True. Only applicable to lightcurves.
:param etresampling: (optional, double)
Exposure time of the resampling (same unit as times). Only used if resampling = True. Only applicable to lightcurves.
:param all_samples: (optional, boolean)
If True, all posterior samples will be used to evaluate the model. Default is False.
Expand Down Expand Up @@ -1915,13 +1906,11 @@ def evaluate_model(self, instrument = None, parameter_values = None, resampling
if not self.global_model:
raise Exception("Input error: an instrument has to be defined for non-global models in order to evaluate the model.")

if (resampling is not None) and (self.modeltype == 'lc') and (instrument is not None):
if resampling:
self.model[instrument]['params'], self.model[instrument]['m'] = init_batman(self.times[instrument], self.dictionary[instrument]['ldlaw'],\
nresampling = nresampling,\
etresampling = etresampling)
else:
self.model[instrument]['params'], self.model[instrument]['m'] = init_batman(self.times[instrument], self.dictionary[instrument]['ldlaw'])
if self.modeltype == 'lc':
nresampling = self.dictionary[instrument].get('nresampling')
etresampling = self.dictionary[instrument].get('exptimeresampling')
self.model[instrument]['params'], self.model[instrument]['m'] = init_batman(self.times[instrument], self.dictionary[instrument]['ldlaw'],
nresampling=nresampling, etresampling=etresampling)

# Save the original inames in the case of non-global models, and set self.inames to the input model. This is because if the model
# is not global, in general we don't care about generating the models for the other instruments (and in the lightcurve and RV evaluation part,
Expand Down Expand Up @@ -1957,8 +1946,8 @@ def evaluate_model(self, instrument = None, parameter_values = None, resampling
idx_samples = np.random.choice(np.arange(nsampled),np.min([nsamples,nsampled]),replace=False)
idx_samples = idx_samples[np.argsort(idx_samples)]

# Create the output_model arrays: these will save on each iteration the full model (lc/rv + GP, output_model_samples),
# the GP-only model (GP, output_modelGP_samples) and the lc/rv-only model (lc/rv, output_modelDET_samples) --- the latter ones
# Create the output_model arrays: these will save on each iteration the full model (lc/rv + GP, output_model_samples),
# the GP-only model (GP, output_modelGP_samples) and the lc/rv-only model (lc/rv, output_modelDET_samples) --- the latter ones
# will make sense only if there is a GP model. If not, it will be a zero-array throughout the evaluation process:
if t is None:
# If user did not give input times, then output samples follow the times on which the model was fitted:
Expand Down Expand Up @@ -1997,22 +1986,36 @@ def evaluate_model(self, instrument = None, parameter_values = None, resampling
# those in dictionaries:
for ginstrument in instruments:
if self.dictionary[ginstrument]['TransitFit'] or self.dictionary[ginstrument]['TransitFitCatwoman']:
nresampling = self.dictionary[ginstrument].get('nresampling')
etresampling = self.dictionary[ginstrument].get('exptimeresampling')
supersample_params, supersample_m = {}, {}
sample_params, sample_m = {}, {}
if not self.dictionary[ginstrument]['TransitFitCatwoman']:
supersample_params[ginstrument],supersample_m[ginstrument] = init_batman(t, self.dictionary[ginstrument]['ldlaw'])
sample_params[ginstrument],sample_m[ginstrument] = init_batman(self.times[ginstrument], self.dictionary[ginstrument]['ldlaw'])
supersample_params[ginstrument],supersample_m[ginstrument] = init_batman(t, self.dictionary[ginstrument]['ldlaw'],
nresampling=nresampling, etresampling=etresampling)
sample_params[ginstrument],sample_m[ginstrument] = init_batman(self.times[ginstrument], self.dictionary[ginstrument]['ldlaw'],
nresampling=nresampling, etresampling=etresampling)
else:
supersample_params[ginstrument],supersample_m[ginstrument] = init_catwoman(t, self.dictionary[ginstrument]['ldlaw'])
sample_params[ginstrument],sample_m[ginstrument] = init_catwoman(self.times[ginstrument], self.dictionary[ginstrument]['ldlaw'])
supersample_params[ginstrument],supersample_m[ginstrument] = init_catwoman(t, self.dictionary[ginstrument]['ldlaw'],
nresampling=nresampling, etresampling=etresampling)
sample_params[ginstrument],sample_m[ginstrument] = init_catwoman(self.times[ginstrument], self.dictionary[ginstrument]['ldlaw'],
nresampling=nresampling, etresampling=etresampling)
else:
# If model is not global, the variables saved are not dictionaries but simply the objects, as we are just going to evaluate the
# model for one dataset (the one of the input instrument):
if self.dictionary[instrument]['TransitFit'] or self.dictionary[instrument]['TransitFitCatwoman']:
nresampling = self.dictionary[instrument].get('nresampling')
etresampling = self.dictionary[instrument].get('exptimeresampling')
if not self.dictionary[instrument]['TransitFitCatwoman']:
supersample_params,supersample_m = init_batman(t, self.dictionary[instrument]['ldlaw'])
sample_params,sample_m = init_batman(self.times[instrument], self.dictionary[instrument]['ldlaw'])
supersample_params,supersample_m = init_batman(t, self.dictionary[instrument]['ldlaw'],
nresampling=nresampling, etresampling=etresampling)
sample_params,sample_m = init_batman(self.times[instrument], self.dictionary[instrument]['ldlaw'],
nresampling=nresampling, etresampling=etresampling)
else:
supersample_params,supersample_m = init_catwoman(t, self.dictionary[instrument]['ldlaw'])
sample_params,sample_m = init_catwoman(self.times[instrument], self.dictionary[instrument]['ldlaw'])
supersample_params,supersample_m = init_catwoman(t, self.dictionary[instrument]['ldlaw'],
nresampling=nresampling, etresampling=etresampling)
sample_params,sample_m = init_catwoman(self.times[instrument], self.dictionary[instrument]['ldlaw'],
nresampling=nresampling, etresampling=etresampling)
else:
# If we are trying to evaluate radial-velocities, we don't need to generate objects because radvel receives the times as inputs
# on each call. In this case then we save the original times (self.t has *all* the times of all the instruments) and instrument
Expand Down Expand Up @@ -2386,11 +2389,10 @@ def evaluate_model(self, instrument = None, parameter_values = None, resampling
if self.lm_boolean[instrument]:
components['lm'] = self.model[instrument]['LM']
else:

x = self.evaluate_model(instrument = instrument, parameter_values = self.posteriors, resampling = resampling, \
nresampling = nresampling, etresampling = etresampling, all_samples = all_samples, \
nsamples = nsamples, return_samples = return_samples, t = t, GPregressors = GPregressors, \
LMregressors = LMregressors, return_err = return_err, return_components = return_components, alpha = alpha, \

x = self.evaluate_model(instrument = instrument, parameter_values = self.posteriors, all_samples = all_samples,
nsamples = nsamples, return_samples = return_samples, t = t, GPregressors = GPregressors,
LMregressors = LMregressors, return_err = return_err, return_components = return_components, alpha = alpha,
evaluate_transit = evaluate_transit)
if return_samples:
if return_err:
Expand All @@ -2415,15 +2417,6 @@ def evaluate_model(self, instrument = None, parameter_values = None, resampling
else:
output_model = x

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

if not self.global_model:
# Return original inames back in case of non-global models:
self.inames = original_inames
Expand Down Expand Up @@ -2710,7 +2703,7 @@ def set_posterior_samples(self, posterior_samples):
self.posteriors = posterior_samples
self.median_posterior_samples = {}
for parameter in self.posteriors.keys():
if parameter is not 'unnamed':
if parameter != 'unnamed':
self.median_posterior_samples[parameter] = np.median(self.posteriors[parameter])
for parameter in self.priors:
if self.priors[parameter]['distribution'] == 'fixed':
Expand Down

0 comments on commit 162f197

Please sign in to comment.