Kernel Building¶
The class of kernels that can be used in celerite are sums of terms of the following form:
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, andbeta_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:
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 Modeling Protocol, 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()
andterms.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()
andterms.Term.get_complex_coefficients()
but it will always return a tuple of length 6, even ifalpha_complex_imag
was omitted fromget_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
, andbeta_complex_imag
as described above.Raises: ValueError
– For invalid dimensions for the coefficients.
-
get_complex_coefficients
()¶ Get the arrays
alpha_complex_*
andbeta_complex_*
This method should be overloaded by subclasses to return the arrays
alpha_complex_real
,alpha_complex_imag
,beta_complex_real
, andbeta_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
, andbeta_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 asomega
.
-
get_real_coefficients
()¶ Get the arrays
alpha_real
andbeta_real
This method should be overloaded by subclasses to return the arrays
alpha_real
andbeta_real
given the current parameter settings. By default, this term is empty.Returns: alpha_real
andbeta_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 astau
.
-
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
andlog_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:
-
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
, andlog_d
. The parameterlog_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:
-
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
, andlog_omega0
.Parameters:
-
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
andlog_rho
. The parametereps
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: