# 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:log_a', 0.5),
('terms:log_c', -5.0),
('terms:log_a', 1.0),
('terms:log_b', 0.1),
('terms:log_c', 0.6),
('terms: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. 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. 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. (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. 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)