Covariance Routines

In the original paper, the kernel was specified using the velocity distance between two wavelengths (\(\lambda_i\), \(\lambda_j\)). In subsequent versions, we choose to simply use the absolute magnitude of the distance between the natural log of the wavelengths, which is functionally equivalent to the former definition, but makes the computation easier. We let

\[\zeta_i = \ln \lambda_i\]

\[r_{ij} = |\zeta_i - \zeta_j|\]

Of course, when we are fitting for two spectral components, f and g, then there will be two sets of input wavelength vectors and thus two distances

\[\begin{split}r_{f,ij} = |\zeta_{f,i} - \zeta_{f,j}| \\ r_{g,ij} = |\zeta_{g,i} - \zeta_{g,j}|\end{split}\]

With two stars, to evaluate the kernel for a specific element in the matrix actually requires four inputs

\[k(\zeta_{f,i}, \zeta_{f,j}, \zeta_{g,i}, \zeta_{g,j} | a_f, l_f, a_g, l_g) = a_f^2 \exp \left ( - \frac{|\zeta_{f,i} - \zeta_{f,j}|^2}{2 l_f^2} \right) + a_g^2 \exp \left ( - \frac{|\zeta_{g,i} - \zeta_{g,j}|^2}{2 l_g^2} \right)\]

psoap.covariance.cycle_calibration(wl, fl, sigma, amp_f, l_f, ncycles, order=1, limit_array=3, mu_GP=1.0, soften=1.0)

Given a chunk of spectra, cycle n_cycles amongst all spectra and return the spectra with inferred calibration adjustments.

order : what order of Chebyshev polynomials to use. 1st order = line.

Only use limit_array number of spectra to save memory.

psoap.covariance.cycle_calibration_chunk(chunk, amp_f, l_f, n_cycles, order=1, limit_array=3, mu_GP=1.0, soften=1.0)

Do the calibration on a chunk at at time, incorporating the masks.

psoap.covariance.lnlike_f(V11, wl_f, fl, sigma, amp_f, l_f, mu_GP=1.0)

Calculate the log-likelihood for a single-lined spectrum.

This function takes a pre-allocated array and fills out the covariance matrices and evaluates the likelihood function for a single-lined spectrum, assuming a squared-exponential kernel (does not celerite).

Parameters:

• V11 (numpy 2D array) – Description of arg1
• wl_f (numpy 1D array) – Description of arg2
• fl (numpy 1D array) – ae
• amp_f (float) – amplitude of GP
• l_f (float) – length scale of GP
• mu_GP (float) – mean of GP

Returns:

The log-likelihood value

Return type:

float

psoap.covariance.lnlike_f_g(V11, wl_f, wl_g, fl, sigma, amp_f, l_f, amp_g, l_g, mu_GP=1.0)

V11 is a matrix to be allocated. wl_known, fl_known, and sigma_known are flattened 1D arrays.

psoap.covariance.lnlike_f_g_george(lwl_f, lwl_g, fl, sigma, amp_f, l_f, amp_g, l_g, mu_GP=1.0)

Evaluate the joint likelihood for f and g using George.

psoap.covariance.lnlike_f_g_h(V11, wl_f, wl_g, wl_h, fl, sigma, amp_f, l_f, amp_g, l_g, amp_h, l_h, mu_GP=1.0)

V11 is a matrix to be allocated. wl_known, fl_known, and sigma_known are flattened 1D arrays.

psoap.covariance.optimize_GP_f(wl_known, fl_known, sigma_known, amp_f, l_f, mu_GP=1.0)

Optimize the GP hyperparameters for the given slice of data. Amp and lv are starting guesses.

psoap.covariance.optimize_calibration(lwl0, lwl1, lwl_cal, fl_cal, fl_fixed, A, B, C, order=1, mu_GP=1.0)

Determine the calibration parameters for this epoch of observations. This is a more general method than psoap.covariance.optimize_calibration_static(), since it allows arbitrary covariance matrices, which should be used when there is orbital motion. Assumes that covariance matrices are appropriately filled out.

Parameters:

• lwl0 (float) – left side evaluation point for Chebyshev
• lwl1 (float) – right side evaluation point for Chebyshev
• lwl_cal (np.array) – the wavelengths corresponding to the epoch we want to calibrate
• fl_cal (np.array) – the fluxes corresponding to the epoch we want to calibrate
• fl_fixed (np.array) – the remaining epochs of data to calibrate in reference to.
• A (2D np.array) – matrix_functions.fill_V11_f(A, lwl_cal, amp, l_f) with sigma_cal already added to the diagonal
• B (2D np.array) – matrix_functions.fill_V11_f(B, lwl_fixed, amp, l_f) with sigma_fixed already added to the diagonal
• C (2D np.array) – matrix_functions.fill_V12_f(C, lwl_cal, lwl_fixed, amp, l_f) cross matrix (with no sigma added, since these are independent measurements).
• order (int) – the degree polynomial to use. order = 1 is a line, order = 2 is a line + parabola

Returns:

a tuple of two data products. The first is the fl_cal vector, now calibrated. The second is the array of the Chebyshev coefficients, in case one wants to re-evaluate the calibration polynomials.

Return type:

(np.array, np.array)

psoap.covariance.optimize_calibration_ST1(lwl0, lwl1, lwl_cal, fl_cal, fl_fixed, gp, A, C, mu_GP=1.0, order=1)

Determine the calibration parameters for this epoch of observations.

lwl0, lwl1: set the points for the Chebyshev.

This is a more general method than optimize_calibration_static, since it allows arbitrary covariance matrices, which should be used when there is orbital motion.

lwl_cal: the wavelengths corresponding to the epoch we want to calibrate fl_cal: the fluxes corresponding to the epoch we want to calibrate

fl_fixed: the remaining epochs of data to calibrate in reference to.

gp: the celerite GP

order: the degree polynomial to use. order = 1 is a line, order = 2 is a line + parabola

Assumes that covariance matrices are appropriately filled out.

psoap.covariance.optimize_calibration_static(wl0, wl1, wl_cal, fl_cal, sigma_cal, wl_fixed, fl_fixed, sigma_fixed, amp, l_f, order=1, mu_GP=1.0)

Determine the calibration parameters for this epoch of observations. Assumes all wl, fl arrays are 1D, and that the relative velocities between all epochs are zero.

Parameters:

• wl0 (float) – left wl point to evaluate the Chebyshev
• wl1 (float) – right wl point to evaluate the Chebyshev
• wl_cal (np.array) – the wavelengths of the epoch to calibrate
• fl_cal (np.array) – the fluxes of the epoch to calibrate
• sigma_cal (np.array) – the sigmas of the epoch to calibrate
• wl_fixed (np.array) – the 1D (flattened) array of the reference wavelengths
• fl_fixed (np.array) – the 1D (flattened) array of the reference fluxes
• sigma_fixed (np.array) – the 1D (flattened) array of the reference sigmas
• amp (float) – the GP amplitude
• l_f (float) – the GP length
• order (int) – the Chebyshev order to use
• mu_GP (optional) – the mean of the GP to assume.

Returns:

a tuple of two data products. The first is the fl_cal vector, now calibrated. The second is the array of the Chebyshev coefficients, in case one wants to re-evaluate the calibration polynomials.

Return type:

(np.array, np.array)

psoap.covariance.optimize_epoch_velocity_f(lwl_epoch, fl_epoch, sigma_epoch, lwl_fixed, fl_fixed, sigma_fixed, gp)

Optimize the wavelengths of the chosen epoch relative to the fixed wavelengths. Identify the velocity required to redshift the chosen epoch.

psoap.covariance.predict_f(lwl_known, fl_known, sigma_known, lwl_predict, amp_f, l_f, mu_GP=1.0)

wl_known are known wavelengths. wl_predict are the prediction wavelengths. Assumes all inputs are 1D arrays.

psoap.covariance.predict_f_g(lwl_f, lwl_g, fl_fg, sigma_fg, lwl_f_predict, lwl_g_predict, mu_f, amp_f, l_f, mu_g, amp_g, l_g, get_Sigma=True)

Given that f + g is the flux that we’re modeling, jointly predict the components.

psoap.covariance.predict_f_g_h(lwl_f, lwl_g, lwl_h, fl_fgh, sigma_fgh, lwl_f_predict, lwl_g_predict, lwl_h_predict, mu_f, mu_g, mu_h, amp_f, l_f, amp_g, l_g, amp_h, l_h)

Given that f + g + h is the flux that we’re modeling, jointly predict the components.

psoap.covariance.predict_f_g_h_sum(lwl_f, lwl_g, lwl_h, fl_fgh, sigma_fgh, lwl_f_predict, lwl_g_predict, lwl_h_predict, mu_fgh, amp_f, l_f, amp_g, l_g, amp_h, l_h)

Given that f + g + h is the flux that we’re modeling, predict the joint sum.

psoap.covariance.predict_python(wl_known, fl_known, sigma_known, wl_predict, amp_f, l_f, mu_GP=1.0)

wl_known are known wavelengths. wl_predict are the prediction wavelengths.

This is the covariance module.