Passing parameter_values to evaluate.rv doesn't work

[melissa-hobson]

I'm trying to pass specific parameters to rv.evaluate, and getting the following traceback:

---> 93 full_model, components = results.rv.evaluate(instrument='FEROS', t = model_times, GPregressors = model_times, return_components = True, parameter_values = posterior_params_med)

File ~/juliet/juliet/fit.py:2940, in model.evaluate_model(self, instrument, parameter_values, all_samples, nsamples, return_samples, t, GPregressors, LMregressors, return_err, alpha, return_components, evaluate_transit)
2937 self.variances = self.model['global_variances']
2938 if self.dictionary['global_model']['GPDetrend']:
2939 self.model['deterministic'], self.model[
-> 2940 'GP'], output_model = self.get_GP_plus_deterministic_model(
2941 parameter_values)
2942 else:
2943 output_model = self.get_GP_plus_deterministic_model(
2944 parameter_values)

File ~/juliet/juliet/fit.py:2203, in model.get_GP_plus_deterministic_model(self, parameter_values, instrument)
2198 self.dictionary['global_model'][
2199 'noise_model'].set_parameter_vector(parameter_values)
2200 self.dictionary['global_model']['noise_model'].yerr = np.sqrt(
2201 self.variances)
2202 self.dictionary['global_model']['noise_model'].compute_GP(
-> 2203 X=self.original_GPregressors)
2204 # Return mean signal plus GP model:
2205 self.model['GP'] = self.dictionary['global_model']['noise_model'].GP.predict(self.residuals, self.dictionary['global_model']['noise_model'].X,
2206 return_var=False, return_cov=False)

AttributeError: 'model' object has no attribute 'original_GPregressors'

Here is a minimal working example using the TOI-141 data from the docs:

import numpy as np
import juliet

def get_quantiles(dist, alpha=0.68, method='median'):
    """
    get_quantiles function
    DESCRIPTION
        This function returns, in the default case, the parameter median and the error%
        credibility around it. This assumes you give a non-ordered
        distribution of parameters.
    OUTPUTS
        Median of the parameter,upper credibility bound, lower credibility bound
    """
    ordered_dist = dist[np.argsort(dist)]
    param = 0.0
    # Define the number of samples from posterior
    nsamples = len(dist)
    nsamples_at_each_side = int(nsamples * (alpha / 2.) + 1)
    if (method == 'median'):
        med_idx = 0
        if (nsamples % 2 == 0.0):  # Number of points is even
            med_idx_up = int(nsamples / 2.) + 1
            med_idx_down = med_idx_up - 1
            param = (ordered_dist[med_idx_up] + ordered_dist[med_idx_down]) / 2.
            return param, ordered_dist[med_idx_up + nsamples_at_each_side], \
                ordered_dist[med_idx_down - nsamples_at_each_side]
        else:
            med_idx = int(nsamples / 2.)
            param = ordered_dist[med_idx]
            return param, ordered_dist[med_idx + nsamples_at_each_side], \
                ordered_dist[med_idx - nsamples_at_each_side]


priors = {}

# Name of the parameters to be fit:
params = ['P_p1','t0_p1','mu_CORALIE14', \
          'mu_CORALIE07','mu_HARPS','mu_FEROS',\
          'K_p1', 'ecc_p1', 'omega_p1', 'sigma_w_CORALIE14','sigma_w_CORALIE07',\
           'sigma_w_HARPS','sigma_w_FEROS','GP_sigma_rv','GP_rho_rv']

# Distributions:
dists = ['normal','normal','uniform', \
         'uniform','uniform','uniform',\
         'uniform','fixed', 'fixed', 'loguniform', 'loguniform',\
         'loguniform', 'loguniform','loguniform','loguniform']

# Hyperparameters
hyperps = [[1.007917,0.000073], [2458325.5386,0.0011], [-100,100], \
           [-100,100], [-100,100], [-100,100], \
           [0.,100.], 0., 90., [1e-3, 100.], [1e-3, 100.], \
           [1e-3, 100.], [1e-3, 100.],[0.01,100.],[0.01,100.]]

# Populate the priors dictionary:
for param, dist, hyperp in zip(params, dists, hyperps):
    priors[param] = {}
    priors[param]['distribution'], priors[param]['hyperparameters'] = dist, hyperp

# Add second planet to the prior:
params = params + ['P_p2',   't0_p2',  'K_p2',    'ecc_p2','omega_p2']
dists = dists +   ['uniform','uniform','uniform', 'fixed', 'fixed']
hyperps = hyperps + [[1.,10.],[2458325.,2458330.],[0.,100.], 0., 90.]

# Repopulate priors dictionary:
priors = {}

for param, dist, hyperp in zip(params, dists, hyperps):
    priors[param] = {}
    priors[param]['distribution'], priors[param]['hyperparameters'] = dist, hyperp

dataset = juliet.load(priors = priors, rvfilename='rvs_toi141.dat', out_folder = 'toi141_rvs-global', \
                      GPrveparamfile='GP_regressors_rv_toi141.dat')

results = dataset.fit(n_live_points = 300)

# Define minimum and maximum times to evaluate the model on:
min_time, max_time = np.min(dataset.times_rv['FEROS'])-30,\
                 np.max(dataset.times_rv['FEROS'])+30

# Create model times on which we will evaluate the model:
model_times = np.linspace(min_time,max_time,5000)

# evaluate the model on the median posterior parameters
posterior_params_med = {}
for param in params:
    try:
        posterior_params_med[param]= get_quantiles(results.posteriors['posterior_samples'][param])[0]
    except:
        posterior_params_med[param] = priors[param]['hyperparameters']


# Extract full model and components of the RV model:
full_model, components = results.rv.evaluate(instrument='FEROS', t = model_times, GPregressors = model_times, return_components = True, parameter_values = posterior_params_med)

What I think is happening, is that when parameter_values are given, it is at some point "jumping" from treating it as an instrument model to treating it as a global model? Because the original_GPregressors attribute is only defined within evaluate_model for a global model. However, I get the same traceback when not passing an instrument, so letting it be a global model is not the solution either...