Can't plot RV when GPregressors is involved

[LucaNap]

Hello, I am having an issue trying to plot the fit of radial velocities only when (rv) GPs are involved. While trying to replicate the example in Juliet's documentation:

# 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)

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

import matplotlib.pyplot as plt
instruments = ['HARPS','FEROS']
colors = ['red','black']

fig = plt.figure(figsize=(10,4))
for instrument,color in zip (instruments,colors):
    plt.errorbar(dataset.times_rv[instrument]-2454705,dataset.data_rv[instrument] - components['mu'][instrument], \
                 yerr = dataset.errors_rv[instrument], fmt = 'o', label = instrument+' data',mfc='white', mec = color, ecolor = color, \
                 elinewidth=1)

plt.plot(model_times-2454705,full_model - components['mu']['FEROS'],label='Full model',color='black')
plt.plot(model_times-2454705,results.rv.model['deterministic'],label = 'Keplerian component', color = 'steelblue')
plt.plot(model_times-2454705,results.rv.model['GP'], label = 'GP component',color='red')
plt.xlim([3701,3715])
plt.ylabel('Radial velocity (m/s)')
plt.xlabel('Time (BJD - 2454705)')
plt.legend(ncol=2)

I get the following error:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-41-2ae74457e435> in <module>
      9 
     10 # Extract full model and components of the RV model:
---> 11 full_model, components = results.rv.evaluate('FEROS', t = model_times, GPregressors = model_times, return_components = True)
     12 
     13 import matplotlib.pyplot as plt

D:\anaconda3\lib\site-packages\juliet\fit.py in evaluate_model(self, instrument, parameter_values, resampling, nresampling, etresampling, all_samples, nsamples, return_samples, t, GPregressors, LMregressors, return_err, alpha, return_components, evaluate_transit)
   2388         else:
   2389 
-> 2390             x = self.evaluate_model(instrument = instrument, parameter_values = self.posteriors, resampling = resampling, \
   2391                                               nresampling = nresampling, etresampling = etresampling, all_samples = all_samples, \
   2392                                               nsamples = nsamples, return_samples = return_samples, t = t, GPregressors = GPregressors, \

D:\anaconda3\lib\site-packages\juliet\fit.py in evaluate_model(self, instrument, parameter_values, resampling, nresampling, etresampling, all_samples, nsamples, return_samples, t, GPregressors, LMregressors, return_err, alpha, return_components, evaluate_transit)
   2116                         self.generate_lc_model(current_parameter_values, evaluate_lc = True)
   2117                     else:
-> 2118                         self.generate_rv_model(current_parameter_values, evaluate_global_errors = True)
   2119 
   2120                     # Save residuals (and global errors, in the case of global models):

D:\anaconda3\lib\site-packages\juliet\fit.py in generate_rv_model(self, parameter_values, evaluate_global_errors)
   1802                 self.model['global'][self.instrument_indexes[instrument]] = self.model[instrument]['deterministic']
   1803                 if evaluate_global_errors:
-> 1804                     self.model['global_variances'][self.instrument_indexes[instrument]] = self.yerr[self.instrument_indexes[instrument]]**2 + \
   1805                                                                                           parameter_values['sigma_w_'+instrument]**2
   1806 

IndexError: index 59 is out of bounds for axis 0 with size 59

It looks like the error is related to GPregressors, but of course the evaluation can't be done without it.

[LucaNap]

I have found out that, with a GP model, everytime results.rv.evaluate() is called the "results" are updated and this makes it impossible to call the function again (if for example one is trying to plot points from two or more instruments, or if one is trying to make more than one plot using results.rv.evaluate).

P.S. Even more strange is that this error can't be bypassed by making a copy of "results" (dataset.fit) because the copy also gets overwritten after results.rv.evaluate() is called.

[LucaNap]

Update: I've found a workaround, which relies on using

import copy
results_1 = copy.deepcopy(results)

everytime results.rv.evaluate is called (now becomes results_1).

[nespinoza]

Hi @LucaNap,

Yes, this is actually not a desirable output of calling the evaluate method for sure, and is in my to-do to look at. I hope to have it fixed for the next juliet version; will thus leave this open for now!

Thanks for bringing this up with such a detailed report!

Néstor

[nespinoza]

Finally got to this -- fixed on dev. Will be on the next juliet version --- but users that need this/find this problem can go and install the dev version.