Kernel Building¶

The class of kernels that can be used in celerite are sums of terms of the following form:

\[\begin{split}k(\tau) &= \sum_{j=1}^J \frac{1}{2}\left[ (a_j + i\,b_j)\,e^{-(c_j+i\,d_j)\,\tau} + (a_j - i\,b_j)\,e^{-(c_j-i\,d_j)\,\tau} \right] \\ &= \sum_{j=1}^J \frac{1}{2}\left[ \alpha_j\,e^{-\beta_j\,\tau} + {\alpha_j}^*\,e^{-{\beta_j}^*\,\tau} \right]\end{split}\]

At a basic level, the interface involves providing the solver with lists of \(\alpha_j\) and \(\beta_j\). It’s useful to know if \(b_j = 0\) and \(d_j = 0\) for any terms because we can use that fact to gain computational efficiency. Therefore, there are 6 arrays that should be provided to the code:

alpha_real: an array of \(a_j\) for the terms where \(b_j = 0\) and \(d_j = 0\),
beta_real: an array of \(c_j\) for the terms where \(b_j = 0\) and \(d_j = 0\),
alpha_complex_real: an array of \(a_j\) for the other terms,
alpha_complex_imag: an array of \(b_j\) for the other terms,
beta_complex_real: an array of \(c_j\) for the other terms, and
beta_complex_imag: an array of \(d_j\) for the other terms.

In practice, this is further abstracted and users should use the provided term objects described below or subclass Term to implement a custom kernel. Subclasses should overload the Term.get_real_coefficients() and Term.get_complex_coefficients() methods to provide custom \(\alpha\) and \(\beta\) arrays.

The following gives a simple example for a custom term of the form:

\[k(\tau) = \frac{a}{2+b}\,e^{-c\,\tau}\left[ \cos\left(\frac{2\,\pi\,\tau}{P}\right) + (1+b) \right]\]

with the parameters log_a, log_b, log_c, and log_P.

class CustomTerm(terms.Term):
    parameter_names = ("log_a", "log_b", "log_c", "log_P")

    def get_real_coefficients(self):
        b = np.exp(self.log_b)
        return (
            np.exp(self.log_a) * (1.0 + b) / (2.0 + b),
            np.exp(self.log_c),
        )

    def get_complex_coefficients(self):
        b = np.exp(self.log_b)
        return (
            np.exp(self.log_a) / (2.0 + b),
            0.0,
            np.exp(self.log_c),
            2*np.pi*np.exp(-self.log_P),
        )

In this example, all of the returned coefficients are scalars but they can also be returned as arrays.

Finally, terms are combined by adding them:

term1 = terms.RealTerm(0.5, -5.0)
term2 = terms.ComplexTerm(1.0, 0.1, 0.6, -0.3)
kernel = term1 + term2

This sum is also a Term and it provides the same interface, concatenating the parameters and coefficients in the correct order.

Since Term objects implement the Python: Modeling, for the kernel described above, we can do things like the following:

print(kernel.get_parameter_dict())

With the output:

OrderedDict([('terms[0]:log_a', 0.5),
             ('terms[0]:log_c', -5.0),
             ('terms[1]:log_a', 1.0),
             ('terms[1]:log_b', 0.1),
             ('terms[1]:log_c', 0.6),
             ('terms[1]:log_d', -0.2)])

The following are the Term objects provided by celerite.

class celerite.terms.Term(*args, **kwargs)¶

The abstract base “term” that is the superclass of all other terms

Subclasses should overload the terms.Term.get_real_coefficients() and terms.Term.get_complex_coefficients() methods.

check_parameters()¶

Check for negative power in the PSD using Sturm’s theorem

Returns:	`True` for valid parameters.

coefficients¶

All of the coefficient arrays

This property is the concatenation of the results from terms.Term.get_real_coefficients() and terms.Term.get_complex_coefficients() but it will always return a tuple of length 6, even if alpha_complex_imag was omitted from get_complex_coefficients.

Returns:	(array[j_real], array[j_real], array[j_complex], array[j_complex], array[j_complex], array[j_complex]): `alpha_real`, `beta_real`, `alpha_complex_real`, `alpha_complex_imag`, `beta_complex_real`, and `beta_complex_imag` as described above.
Raises:	`ValueError` – For invalid dimensions for the coefficients.

get_complex_coefficients()¶

Get the arrays alpha_complex_* and beta_complex_*

This method should be overloaded by subclasses to return the arrays alpha_complex_real, alpha_complex_imag, beta_complex_real, and beta_complex_imag given the current parameter settings. By default, this term is empty.

Returns:	(array[j_complex], array[j_complex], array[j_complex], array[j_complex]): `alpha_complex_real`, `alpha_complex_imag`, `beta_complex_real`, and `beta_complex_imag` as described above. `alpha_complex_imag` can be omitted and it will be assumed to be zero.

get_psd(omega)¶

Compute the PSD of the term for an array of angular frequencies

Parameters:	omega (array[..]) – An array of frequencies where the PSD should be evaluated.
Returns:	The value of the PSD for each `omega`. This will have the same shape as `omega`.

get_real_coefficients()¶

Get the arrays alpha_real and beta_real

This method should be overloaded by subclasses to return the arrays alpha_real and beta_real given the current parameter settings. By default, this term is empty.

Returns:	`alpha_real` and `beta_real` as described above.
Return type:	(array[j_real], array[j_real])

get_value(tau)¶

Compute the value of the term for an array of lags

Parameters:	tau (array[..]) – An array of lags where the term should be evaluated.
Returns:	The value of the term for each `tau`. This will have the same shape as `tau`.

terms¶: A list of all the terms included in a sum of terms

class celerite.terms.RealTerm(*args, **kwargs)¶

The simplest celerite term

This term has the form

\[k(\tau) = a_j\,e^{-c_j\,\tau}\]

with the parameters log_a and log_c.

Strictly speaking, for a sum of terms, the parameter a could be allowed to go negative but since it is somewhat subtle to ensure positive definiteness, we require that the amplitude be positive through this interface. Advanced users can build a custom term that has negative coefficients but care should be taken to ensure positivity.

Parameters:	log_a (float) – The log of the amplitude of the term. log_c (float) – The log of the exponent of the term.

class celerite.terms.ComplexTerm(*args, **kwargs)¶

A general celerite term

This term has the form

\[k(\tau) = \frac{1}{2}\,\left[(a_j + b_j)\,e^{-(c_j+d_j)\,\tau} + (a_j - b_j)\,e^{-(c_j-d_j)\,\tau}\right]\]

with the parameters log_a, log_b, log_c, and log_d. The parameter log_b can be omitted and it will be assumed to be zero.

This term will only correspond to a positive definite kernel (on its own) if \(a_j\,c_j \ge b_j\,d_j\) and the log_prior method checks for this constraint.

Parameters:	log_a (float) – The log of the real part of amplitude. log_b (float) – The log of the imaginary part of amplitude. log_c (float) – The log of the real part of the exponent. log_d (float) – The log of the imaginary part of exponent.

class celerite.terms.SHOTerm(*args, **kwargs)¶

A term representing a stochastically-driven, damped harmonic oscillator

The PSD of this term is

\[S(\omega) = \sqrt{\frac{2}{\pi}} \frac{S_0\,\omega_0^4} {(\omega^2-{\omega_0}^2)^2 + {\omega_0}^2\,\omega^2/Q^2}\]

with the parameters log_S0, log_Q, and log_omega0.

Parameters:	log_S0 (float) – The log of the parameter \(S_0\). log_Q (float) – The log of the parameter \(Q\). log_omega0 (float) – The log of the parameter \(\omega_0\).

class celerite.terms.Matern32Term(*args, **kwargs)¶

A term that approximates a Matern-3/2 function

This term is defined as

\[k(\tau) = \sigma^2\,\left[ \left(1+1/\epsilon\right)\,e^{-(1-\epsilon)\sqrt{3}\,\tau/\rho} \left(1-1/\epsilon\right)\,e^{-(1+\epsilon)\sqrt{3}\,\tau/\rho} \right]\]

with the parameters log_sigma and log_rho. The parameter eps controls the quality of the approximation since, in the limit \(\epsilon \to 0\) this becomes the Matern-3/2 function

\[\lim_{\epsilon \to 0} k(\tau) = \sigma^2\,\left(1+ \frac{\sqrt{3}\,\tau}{\rho}\right)\, \exp\left(-\frac{\sqrt{3}\,\tau}{\rho}\right)\]

Parameters:	log_sigma (float) – The log of the parameter \(\sigma\). log_rho (float) – The log of the parameter \(\rho\). eps (Optional[float]) – The value of the parameter \(\epsilon\). (default: 0.01)