C++ interface¶
All of the computationally expensive parts of celerite2 are implemented in C++. This interface was not designed to be as user friendly as the Python interfaces and it isn’t as thoroughly documented, but these details should be enough to get you started if you’re interested in using or contributing to the low-level code directly. Please feel free to open an issue if you have questions!
One thing to note when using these functions is that the performance will be
greatly improved (especially for small values of J
) by using fixed size
matrices because Eigen can work its magic. If
you use the Terms interface defined below to construct your low rank
matrices, the results will be right (using the Width parameter), but you might
need to enforce this manually if you’re building your own matrices. Much of the
code in the pybind11 interface is dedicated
to making sure that small matrices correctly typed.
It is also worth noting that the notation in the equations below (and throughout
these pages) is a little sloppy since it generally neglects the P
matrix.
Take a look at the original celerite paper
for more details about the numerical stability introduced by using this matrix,
and how it affects the algorithms.
Basic interface¶
If you don’t need access to the backpropagation functions, you can use the following functions to compute the standard celerite operations. If you need support for backpropagation, use the functions defined in the Autodiff Interface section, but know that they will require a larger memory footprint.
-
template<typename
Input
, typenameCoeffs
, typenameDiag
, typenameLowRank
, typenameDense
>
voidcelerite2::core
::
to_dense
(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<Diag> &a, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, Eigen::MatrixBase<Dense> const &K_out)¶ Get the dense representation of a celerite matrix.
- Parameters
t
: (N,): The input coordinates (must be sorted)c
: (J,): The transport coefficientsa
: (N,): The diagonal componentU
: (N, J): The first low rank matrixV
: (N, J): The second low rank matrixK_out
: (N, N): The dense matrix
-
template<typename
Input
, typenameCoeffs
, typenameDiag
, typenameLowRank
, typenameDiagOut
, typenameLowRankOut
>
Eigen::Indexcelerite2::core
::
factor
(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<Diag> &a, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, Eigen::MatrixBase<DiagOut> const &d_out, Eigen::MatrixBase<LowRankOut> const &W_out)¶ Compute the Cholesky factorization of the system.
This computes
d
andW
such that:diag(a) + tril(U*V^T) + triu(V*U^T) = L*diag(d)*L^T
where
L = 1 + tril(U*W^T)
This can be safely applied in place:
d_out
can point toa
andW_out
can point toV
, and the memory will be reused.- Parameters
t
: (N,): The input coordinates (must be sorted)c
: (J,): The transport coefficientsa
: (N,): The diagonal componentU
: (N, J): The first low rank matrixV
: (N, J): The second low rank matrixd_out
: (N,): The diagonal component of the Cholesky factorW_out
: (N, J): The second low rank component of the Cholesky factor
-
template<typename
Input
, typenameCoeffs
, typenameLowRank
, typenameRightHandSide
, typenameRightHandSideOut
>
voidcelerite2::core
::
solve_lower
(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)¶ Compute the solution of a lower triangular linear equation.
This computes
Z
such that:Y = L * Y
where
L = 1 + tril(U*W^T)
This can be safely applied in place.
- Parameters
t
: (N,): The input coordinates (must be sorted)c
: (J,): The transport coefficientsU
: (N, J): The first low rank matrixW
: (N, J): The second low rank matrixY
: (N, Nrhs): The right hand sideZ_out
: (N, Nrhs): The solution of this equation
-
template<typename
Input
, typenameCoeffs
, typenameLowRank
, typenameRightHandSide
, typenameRightHandSideOut
>
voidcelerite2::core
::
solve_upper
(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)¶ Compute the solution of a upper triangular linear equation.
This computes
Z
such that:Y = L^T * Y
where
L = 1 + tril(U*W^T)
This can be safely applied in place.
- Parameters
t
: (N,): The input coordinates (must be sorted)c
: (J,): The transport coefficientsU
: (N, J): The first low rank matrixW
: (N, J): The second low rank matrixY
: (N, Nrhs): The right hand sideZ_out
: (N, Nrhs): The solution of this equation
-
template<typename
Input
, typenameCoeffs
, typenameLowRank
, typenameRightHandSide
, typenameRightHandSideOut
>
voidcelerite2::core
::
matmul_lower
(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)¶ Apply a strictly lower matrix multiply.
This computes:
Z += tril(U * V^T) * Y
where
tril
is the strictly lower triangular function.Note that this will update the value of
Z
.- Parameters
t
: (N,): The input coordinates (must be sorted)c
: (J,): The transport coefficientsU
: (N, J): The first low rank matrixV
: (N, J): The second low rank matrixY
: (N, Nrhs): The matrix to be multipliedZ_out
: (N, Nrhs): The matrix to be updated
-
template<typename
Input
, typenameCoeffs
, typenameLowRank
, typenameRightHandSide
, typenameRightHandSideOut
>
voidcelerite2::core
::
matmul_upper
(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)¶ Apply a strictly upper matrix multiply.
This computes:
Z += triu(V * U^T) * Y
where
triu
is the strictly lower triangular function.Note that this will update the value of
Z
.- Parameters
t
: (N,): The input coordinates (must be sorted)c
: (J,): The transport coefficientsU
: (N, J): The first low rank matrixV
: (N, J): The second low rank matrixY
: (N, Nrhs): The matrix to be multipliedZ_out
: (N, Nrhs): The matrix to be updated
-
template<typename
Input
, typenameCoeffs
, typenameLowRank
, typenameRightHandSide
, typenameRightHandSideOut
>
voidcelerite2::core
::
general_matmul_lower
(const Eigen::MatrixBase<Input> &t1, const Eigen::MatrixBase<Input> &t2, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)¶ The general lower-triangular dot product of a rectangular celerite system.
- Parameters
t1
: (N,): The left input coordinates (must be sorted)t2
: (M,): The right input coordinates (must be sorted)c
: (J,): The transport coefficientsU
: (N, J): The first low rank matrixV
: (M, J): The second low rank matrixY
: (M, Nrhs): The matrix that will be multipliedZ_out
: (N, Nrhs): The result of the operation
-
template<typename
Input
, typenameCoeffs
, typenameLowRank
, typenameRightHandSide
, typenameRightHandSideOut
>
voidcelerite2::core
::
general_matmul_upper
(const Eigen::MatrixBase<Input> &t1, const Eigen::MatrixBase<Input> &t2, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)¶ The general upper-triangular dot product of a rectangular celerite system.
- Parameters
t1
: (N,): The left input coordinates (must be sorted)t2
: (M,): The right input coordinates (must be sorted)c
: (J,): The transport coefficientsU
: (N, J): The first low rank matrixV
: (M, J): The second low rank matrixY
: (M, Nrhs): The matrix that will be multipliedZ_out
: (N, Nrhs): The result of the operation
Autodiff Interface¶
If you don’t need access to the backpropagation functions, you should use the functions defined in the Basic interface section above. The functions defined in this section track the extra matrices that are required for the reverse pass in reverse-mode automatic differentiation. See Foreman-Mackey (2018) for more information about the reverse pass calculations.
Each of these functions follows the same argument naming conventions. As an
example, imagine that we have a function called func
that has two input
arguments, x
(lowercase variables are vectors) and Y
(and uppercase
variables are matrices), and and two outputs A
and b
. Let’s say that the
reverse pass also requires a cached variable (or “workspace”) called S
. In
this case, the signature of the forward pass would be something like:
template <typename VectorIn, typename MatrixIn, typename MatrixOut,
typename VectorOut, typename Workspace>
void func(
const Eigen::MatrixBase<VectorIn> &x,
const Eigen::MatrixBase<MatrixIn> &Y,
Eigen::MatrixBase<MatrixOut> const &A_out,
Eigen::MatrixBase<VectorOut> const &b_out,
Eigen::MatrixBase<Workspace> const &S_out
)
Note that, in general, the _out
parameters don’t need to have the right
shape because they will be resized in place. The reverse pass of this function,
would be implemented in another function func_rev
with the following
signature:
template <typename VectorIn, typename MatrixIn, typename Workspace,
typename VectorOut, typename MatrixOut>
void func_rev(
// Original inputs
const Eigen::MatrixBase<VectorIn> &x,
const Eigen::MatrixBase<MatrixIn> &Y,
// Original outputs
const Eigen::MatrixBase<MatrixIn> &A,
const Eigen::MatrixBase<VectorIn> &b,
const Eigen::MatrixBase<Workspace> &S,
// The sensitivities of the outputs, note: S is not included
const Eigen::MatrixBase<MatrixIn> &bA,
const Eigen::MatrixBase<VectorIn> &bb,
// The (resulting) sensitivities of the inputs
Eigen::MatrixBase<VectorOut> const &bx_out,
Eigen::MatrixBase<MatrixOut> const &bY_out
)
where the b
prefix before a parameter indicates the overbar from the
notation in Foreman-Mackey (2018):
The forward methods for each celerite operation are documented below, and the reverse passes are all implemented following the rules listed above. Take a look at the source code to see the signatures.
-
template<bool
update_workspace
= true, typenameInput
, typenameCoeffs
, typenameDiag
, typenameLowRank
, typenameDiagOut
, typenameLowRankOut
, typenameWork
>
Eigen::Indexcelerite2::core
::
factor
(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<Diag> &a, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, Eigen::MatrixBase<DiagOut> const &d_out, Eigen::MatrixBase<LowRankOut> const &W_out, Eigen::MatrixBase<Work> const &S_out)¶ Compute the Cholesky factorization of the system.
This computes
d
andW
such that:K = L*diag(d)*L^T
where
K
is the celerite matrix andL = 1 + tril(U*W^T)
This can be safely applied in place:
d_out
can point toa
andW_out
can point toV
, and the memory will be reused. In this particular case, thecelerite2::core::factor_rev
function doesn’t usea
andV
, but this won’t be true for all_rev
functions.- Parameters
t
: (N,): The input coordinates (must be sorted)c
: (J,): The transport coefficientsa
: (N,): The diagonal componentU
: (N, J): The first low rank matrixV
: (N, J): The second low rank matrixd_out
: (N,): The diagonal component of the Cholesky factorW_out
: (N, J): The second low rank component of the Cholesky factorS_out
: (N, J*J): The cached value of the S matrix at each step
-
template<bool
update_workspace
= true, typenameInput
, typenameCoeffs
, typenameLowRank
, typenameRightHandSide
, typenameRightHandSideOut
, typenameWork
>
voidcelerite2::core
::
solve_lower
(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out, Eigen::MatrixBase<Work> const &F_out)¶ Apply a strictly lower matrix multiply.
This computes:
Z += tril(U * V^T) * Y
where
tril
is the strictly lower triangular function.Note that this will update the value of
Z
.- Parameters
t
: (N,): The input coordinates (must be sorted)c
: (J,): The transport coefficientsU
: (N, J): The first low rank matrixW
: (N, J): The second low rank matrixY
: (N, Nrhs): The matrix to be multipliedZ_out
: (N, Nrhs): The matrix to be updatedF_out
: (N, J*Nrhs): The workspace
-
template<bool
update_workspace
= true, typenameInput
, typenameCoeffs
, typenameLowRank
, typenameRightHandSide
, typenameRightHandSideOut
, typenameWork
>
voidcelerite2::core
::
solve_upper
(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out, Eigen::MatrixBase<Work> const &F_out)¶ Compute the solution of a upper triangular linear equation.
This computes
Z
such that:Y = L^T * Y
where
L = 1 + tril(U*W^T)
This can be safely applied in place.
- Parameters
t
: (N,): The input coordinates (must be sorted)c
: (J,): The transport coefficientsU
: (N, J): The first low rank matrixW
: (N, J): The second low rank matrixY
: (N, Nrhs): The right hand sideZ_out
: (N, Nrhs): The solution of this equationF_out
: (N, J*Nrhs): The workspace
-
template<bool
update_workspace
= true, typenameInput
, typenameCoeffs
, typenameLowRank
, typenameRightHandSide
, typenameRightHandSideOut
, typenameWork
>
voidcelerite2::core
::
matmul_lower
(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out, Eigen::MatrixBase<Work> const &F_out)¶ Apply a strictly lower matrix multiply.
This computes:
Z += tril(U * V^T) * Y
where
tril
is the strictly lower triangular function.Note that this will update the value of
Z
.- Parameters
t
: (N,): The input coordinates (must be sorted)c
: (J,): The transport coefficientsU
: (N, J): The first low rank matrixV
: (N, J): The second low rank matrixY
: (N, Nrhs): The matrix to be multipliedZ_out
: (N, Nrhs): The matrix to be updatedF_out
: (N, J*Nrhs): The workspace
-
template<bool
update_workspace
= true, typenameInput
, typenameCoeffs
, typenameLowRank
, typenameRightHandSide
, typenameRightHandSideOut
, typenameWork
>
voidcelerite2::core
::
matmul_upper
(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out, Eigen::MatrixBase<Work> const &F_out)¶ Apply a strictly upper matrix multiply.
This computes:
Z += triu(V * U^T) * Y
where
triu
is the strictly lower triangular function.Note that this will update the value of
Z
.- Parameters
t
: (N,): The input coordinates (must be sorted)c
: (J,): The transport coefficientsU
: (N, J): The first low rank matrixV
: (N, J): The second low rank matrixY
: (N, Nrhs): The matrix to be multipliedZ_out
: (N, Nrhs): The matrix to be updatedF_out
: (N, J*Nrhs): The workspace
Terms¶
The only the most basic terms are implemented in C++ and they are mostly used for testing purposes, but it would be possible to generalize them to other use cases.
-
template<typename
T
, intJ_
= Eigen::Dynamic>
classcelerite2
::
Term
¶ The abstract base class from which terms should inherit
Public Types
-
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Width, Order>
LowRank
¶ The
Eigen
type for the low-rank matrices used internally
-
typedef Eigen::Matrix<Scalar, Width, 1>
CoeffVector
¶ The
Eigen
type for a fixed width vector of coefficients
-
typedef std::tuple<Vector, Vector, Vector, Vector, Vector, Vector>
Coeffs
¶ A tuple of vectors giving the coefficients for the celerite model
-
typedef std::tuple<CoeffVector, Vector, LowRank, LowRank>
Matrices
¶ A tuple of matrices representing this celerite process
Public Functions
-
Term
()¶
-
int
get_width
() const¶
-
void
set_coefficients
(const Vector &ar, const Vector &cr, const Vector &ac, const Vector &bc, const Vector &cc, const Vector &dc)¶ Set the coefficients of the term
- Parameters
ar
: (J_real,): The real amplitudes.cr
: (J_real,): The real exponential.ac
: (J_comp,): The complex even amplitude.bc
: (J_comp,): The complex odd amplitude.cc
: (J_comp,): The complex exponential.dc
: (J_comp,): The complex frequency.
-
Matrices
get_celerite_matrices
(const Vector &x, const Vector &diag) const¶ Get the matrices required to represent the celerite process
- Parameters
x
: (N,): The independent coordinates of the data.diag
: (N,): The diagonal variance of the process.
-
template<typename
Other
>
Term<typename std::common_type<Scalar, typename Other::Scalar>::type, sum_width<Width, Other::Width>::value>operator+
(const Other &other) const¶ Adding two terms builds a new term where the coefficients have been concatenated
- Parameters
other
: (Term): The term to add to this one.
-
typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Width, Order>
-
template<typename
T
>
classcelerite2
::
SHOTerm
: public celerite2::Term<T, 2>¶ A term representing a stochastically-driven, damped harmonic oscillator
- Parameters
S0
: The power atomega = 0
.w0
: The undamped angular frequency.Q
: The quality factor.eps
: A regularization parameter used for numerical stability.
Public Types
Public Functions
Public Static Attributes
-
constexpr int
Width
= 2¶
-
template<typename
T
>
classcelerite2
::
RealTerm
: public celerite2::Term<T, 1>¶ The simplest celerite model
- Parameters
a
: The amplitude of the term.c
: The exponent of the term.
Public Types
Public Static Attributes
-
constexpr int
Width
= 1¶
-
template<typename
T
>
classcelerite2
::
ComplexTerm
: public celerite2::Term<T, 2>¶ A general celerite model
- Parameters
a
: The real part of the amplitude.b
: The complex part of the amplitude.c
: The real part of the exponent.d
: The complex part of the exponent.
Public Types
Public Static Attributes
-
constexpr int
Width
= 2¶