Internal documentation

eigen_safe(x)

Compute the eigendecomposition of x. If x is a StaticMatrix, it is converted to a Matrix before the decomposition.

estimated_sectional_curvature(M::AbstractManifold, p, X, Y; r::Real=1e-3, N::Int=10000)

Approximate sectional curvature of manifold M in the plane spanned by vectors X and Y from tangent space at p using a circle on M of radius r divided into N segments.

The approximation is derived from the Bertrand–Diguet–Puiseux theorem which states that

\[\kappa_p(X, Y) = \lim_{r \to 0^+} 3\frac{2\pi r-C(r)}{\pi r^3},\]

where $C(r)$ is the circumference of the circle of radius $r$ around p in submanifold of M spanned by X and Y. The circumference calculation method has a tendency to return curvature values larger than the exact ones. X and Y are orthonormalized internally, so the estimate depends only on the plane they span.

estimated_sectional_curvature_matrix(M::AbstractManifold, p, B::AbstractBasis; r::Real=1e-3, N::Int=10000)

Estimate the matrix of sectional curvatures of manifold M at point p using estimated_sectional_curvature. Entry (i, j) corresponds to sectional curvature of the surface spanned by vectors i and j from basis B.

get_parameter_type(M::AbstractManifold)

Get parameter argument of the constructor of manifold M. Returns either :field or :type.

See also

get_parameter, TypeParameter

isnormal(x; kwargs...) -> Bool

Check if the matrix or number x is normal, that is, if it commutes with its adjoint:

\[x x^\mathrm{H} = x^\mathrm{H} x.\]

By default, this is an equality check. Provide kwargs for isapprox to perform an approximate check.

log_safe(x)

Compute the matrix logarithm of x. If x is a StaticMatrix, it is converted to a Matrix before computing the log.

log_safe!(y, x)

Compute the matrix logarithm of x. If the eltype of y is real, then the imaginary part of x is ignored, and a DomainError is raised if real(x) has no real logarithm.

mul!_safe(Y, A, B) -> Y

Call mul! safely, that is, A and/or B are permitted to alias with Y.

normal_tvector_distribution(M::AbstractManifold, p, σ)

Normal distribution in ambient space with standard deviation σ projected to tangent space at p.

nzsign(z[, absz])

Compute a modified sign(z) that is always nonzero, i.e. where

\[\operatorname{nzsign}(z) = \begin{cases} 1 & \text{if } |z| < ε^2,\\ \frac{z}{|z|} & \text{otherwise.} \end{cases}\]

Here $ε$ is eps of the floating point type of $|z|$, and a given absz replaces $|z|$.

Note that the condition absz == 0 would be incorrectly handled by ForwardDiff.jl.

projected_distribution(M::AbstractManifold, d, [p=rand(d)])

Wrap the standard distribution d into a manifold-valued distribution. Generated points will be of similar type to p. By default, the type is not changed.

realify(X::AbstractMatrix{T𝔽}, 𝔽::AbstractNumbers) -> Y::AbstractMatrix{<:Real}

Given a matrix $X ∈ 𝔽^{n×n}$, compute $Y ∈ ℝ^{m×m}$, where $m = n \operatorname{dim}_𝔽$, and $\operatorname{dim}_𝔽$ is the real_dimension of the number field $𝔽$, using the map $ϕ \colon X ↦ Y$, that preserves the matrix product, so that for all $C,D ∈ 𝔽^{n×n}$,

\[ϕ(C) ϕ(D) = ϕ(CD).\]

See realify! for an in-place version, and unrealify! to compute the inverse of $ϕ$.

realify!(Y::AbstractMatrix{<:Real}, X::AbstractMatrix{T𝔽}, 𝔽::AbstractNumbers)

In-place version of realify.

realify!(Y::AbstractMatrix{<:Real}, X::AbstractMatrix{<:Complex}, ::typeof(ℂ))

Given a complex matrix $X = A + iB ∈ ℂ^{n×n}$, compute its realified matrix $Y ∈ ℝ^{2n×2n}$, written where

\[Y = \begin{pmatrix}A & -B \\ B & A \end{pmatrix}.\]

symmetrize(X)

Given a square matrix X compute 1/2 .* (X' + X).

symmetrize!(Y, X)

Given a square matrix X compute 1/2 .* (X' + X) in place of Y.

unrealify!(X::AbstractMatrix{T𝔽}, Y::AbstractMatrix{<:Real}, 𝔽::AbstractNumbers[, n])

Given a real matrix $Y ∈ ℝ^{m×m}$, where $m = n \operatorname{dim}_𝔽$, and $\operatorname{dim}_𝔽$ is the real_dimension of the number field $𝔽$, compute in-place its equivalent matrix $X ∈ 𝔽^{n×n}$. Note that this function does not check that $Y$ has a valid structure to be un-realified.

See realify! for the inverse of this function.

usinc(θ::Real)

Unnormalized version of sinc function, i.e. $\operatorname{usinc}(θ) = \frac{\sin(θ)}{θ}$. This is equivalent to sinc(θ/π).

Note that ForwardDiff.jl would return wrong answer at θ=0 if a simple equality was used.

usinc_from_cos(x::Real)

Unnormalized version of sinc function, i.e. $\operatorname{usinc}(θ) = \frac{\sin(θ)}{θ}$, computed from $x = cos(θ)$.

vec2skew!(X, v)
vec2skew(v, k)

Create a skew symmetric matrix X of size $k×k$ from a vector v (in place of X), for example for v=[1,2,3] and k=3 this yields

[  0  1  2;  -1  0  3;  -2 -3  0]

(int_term::IntegratorTerminatorNearChartBoundary)(u, t, integrator)

Terminate integration when integrator goes too closely to chart boundary. Closeness is determined by ϵ value of IntegratorTerminatorNearChartBoundary int_term.

Arguments:

• int_term: object containing keyword arguments for check_chart_switch, such as the desired maximum distance to boundary,
• u: parameters of a point at which the integrator is solving a differential equation.
• t: time parameter of the integrator
• integrator: state of the integrator. Internal parameters are expected to contained the manifold on which the equation is solved, the atlas and the current chart index.

StitchedChartSolution{Prob,TM<:AbstractManifold,TA<:AbstractAtlas,TChart}

Solution of an ODE on a manifold M in charts of an AbstractAtlas A.

When StitchedChartSolution{:Exp} is used as a function with a number t as an argument, a pair (p, X) is returned such that $p\in \mathcal{M}$ is the point at time t of the geodesic and $X \in T_p \mathcal{M}$ is the velocity of the geodesic at that point. Similarly, StitchedChartSolution{:PT} called with number t returns a triple (p, X, Y) where (p, X) corresponds to the geodesic along which the vector is transported and $Y\in T_p\mathcal{M}$ is the vector transported to p.

_adjoint_coordinate_map(M::AbstractManifold, A::AbstractAtlas,
    i_from, a_from, L, i_to, a_to, Yc)

Apply the Riemannian adjoint of a linear map represented in chart-induced bases. If L represents a map from the tangent space at the point with coordinates a_from in chart i_from to the tangent space at the point with coordinates a_to in chart i_to, this function returns the coordinates of its adjoint applied to Yc. The adjoint is computed using the local metric matrices as

\[L^* = G_{\mathrm{from}}^{-1}L^\mathsf{T}G_{\mathrm{to}}.\]

_jacobi_exp_matrix(M::AbstractManifold, a, Xc, A::AbstractAtlas, i0; kwargs...)

Compute the coordinate matrix of the differential of the exponential map in a chart. It is obtained by solving the chart-coordinate geodesic and a matrix-valued Jacobi equation, with respect to either the argument (if wrt is set to :argument) or the basepoint (if wrt is set to :basepoint).

The geodesic coordinates satisfy

\[\dot a^k = X^k, \qquad \dot X^k = -\Gamma^k_{ij}(a)X^iX^j.\]

For the matrices $Y$ and $dY$, whose columns are Jacobi fields and their covariant derivatives, respectively, the system is

\[\begin{aligned} \dot Y^k{}_r &= dY^k{}_r - \Gamma^k_{ij}(a)X^iY^j{}_r, \\ \dot{dY}^k{}_r &= -\Gamma^k_{ij}(a)X^i dY^j{}_r - R^k_{\ell ij}(a)Y^\ell{}_rX^iX^j. \end{aligned}\]

The initial conditions are $a(0) = a$, $X(0)$ is set to Xc. If the keyword argument wrt is set to :basepoint, then $Y(0)$ is set to 0, and $dY(0) = I$. If the keyword argument wrt is set to :argument, then $Y(0) = I$, and $dY(0) = 0$. Thus, the returned matrix $Y(1)$ represents $D_X\exp_p(X)$ in the chart-induced bases if wrt is set to :argument and $D_p\exp_p(X)$ if wrt is set to :basepoint. The function also returns the final chart index and the final point coordinates.

solve_chart_adjoint_differential_exp_argument(
    M::AbstractManifold, a, Xc, A::AbstractAtlas, i0, Yc; kwargs...
)

Compute the chart coordinates at the base point of the adjoint of $D_X\exp_p(X)$ applied to Yc. Here Yc contains coordinates in the induced basis of the final chart reached by the geodesic. p is the point with coordinates a in chart i0. X is the tangent vector with coordinates Xc in the induced basis of chart i0 at p.

solve_chart_adjoint_differential_exp_basepoint(
    M::AbstractManifold, a, Xc, A::AbstractAtlas, i0, Yc; kwargs...
)

Compute the chart coordinates at the base point of the adjoint of $D_p\exp_p(X)$ applied to Yc. Here Yc contains coordinates in the induced basis of the final chart reached by the geodesic. p is the point with coordinates a in chart i0. X is the tangent vector with coordinates Xc in the induced basis of chart i0 at p.

solve_chart_adjoint_differential_log_argument(
    M::AbstractManifold, a, Xc, A::AbstractAtlas, i0, Yc; kwargs...
)

Compute the chart coordinates at q = exp_p(X) of the adjoint of $D_q\log_p(q)$ applied to Yc. The input uses the induced basis of the initial chart; the output uses that of the final chart reached by the geodesic. p is the point with coordinates a in chart i0. X is the tangent vector with coordinates Xc in the induced basis of chart i0 at p.

solve_chart_adjoint_differential_log_basepoint(
    M::AbstractManifold, a, Xc, A::AbstractAtlas, i0, Yc; kwargs...
)

Compute the chart coordinates at the base point of the adjoint of $D_p\log_p(q)$, where $q = \exp_p(X)$. Both the input and output use the induced basis of the initial chart. p is the point with coordinates a in chart i0. X is the tangent vector with coordinates Xc in the induced basis of chart i0 at p.

solve_chart_differential_exp_argument(
    M::AbstractManifold, a, Xc, A::AbstractAtlas, i0, Yc; kwargs...
)

Solve the Jacobi equation for $D_X\exp_p(X)[Y]$. The coordinate vectors Xc and Yc are represented in the chart-induced basis at p.

solve_chart_differential_exp_basepoint(
    M::AbstractManifold, a, Xc, A::AbstractAtlas, i0, Yc; kwargs...
)

Solve the Jacobi equation for $D_p\exp_p(X)[Y]$. The coordinate vectors Xc and Yc are represented in the chart-induced basis at p.

The ODE is solved by solve_chart_jacobi_field with dYc set to 0.

solve_chart_differential_log_argument(
    M::AbstractManifold, a, Xc, A::AbstractAtlas, i0, Yc; kwargs...
)

Solve the Jacobi equation for $D_q\log_p(q)[Y]$, where $q = \exp_p(X)$. The coordinate vector Yc is represented in the chart-induced basis at q; the differential is the covariant derivative in solution(0)[4].

solve_chart_differential_log_basepoint(
    M::AbstractManifold, a, Xc, A::AbstractAtlas, i0, Yc; kwargs...
)

Solve the Jacobi equation for $D_p\log_p(q)[Y]$, where $q = \exp_p(X)$. The coordinate vector Yc is represented in the chart-induced basis at p; the differential is the covariant derivative in solution(0)[4].

solve_chart_exp_ode(
    M::AbstractManifold, a, Xc, A::AbstractAtlas, i0;
    solver=AutoVern9(Rodas5P()),
    final_time::Real=1.0,
    check_chart_switch_kwargs=NamedTuple(),
    kwargs...,
)

Solve geodesic ODE on a manifold M from point of coordinates a in chart i0 from an AbstractAtlas A in direction of coordinates Xc in the induced basis. The geodesic is solved up to time final_time (by default equal to 1).

Chart switching

If the solution exceeds the domain of chart i0 (which is detected using the check_chart_switch function with additional keyword arguments check_chart_switch_kwargs), a new chart is selected using get_chart_index on the final point in the old chart.

Returned value

The function returns an object of type StitchedChartSolution{:Exp} to represent the geodesic.

solve_chart_jacobi_field(
    M::AbstractManifold, a, Xc, A::AbstractAtlas, i0, Yc, dYc;
    solver=AutoVern9(Rodas5P()), final_time::Real = 1.0,
    check_chart_switch_kwargs = NamedTuple(), kwargs...
)

Solve the Jacobi equation along the geodesic starting at parameters a in chart i0 with initial velocity coordinates Xc. Yc and dYc are, respectively, the coordinates of the initial Jacobi field and its initial covariant derivative in the induced basis of the chart.

In chart coordinates, this solves the system

\[\begin{aligned} \dot a^k &= X^k, & \dot X^k &= -\Gamma^k_{ij}(a)X^iX^j, \\ \dot Y^k &= dY^k - \Gamma^k_{ij}(a)X^iY^j, & \dot{dY}^k &= -\Gamma^k_{ij}(a)X^i dY^j - R^k_{\ell ij}(a)Y^\ell X^iX^j, \end{aligned}\]

with initial conditions $a(0) = a$, $X(0)$ is equal to Xc, $Y(0)$ is equal to Yc, and $dY(0) = dYc$. Here, dY represents the coordinates of $\nabla_{\dot\gamma}Y$. $\Gamma^k_{ij}$ are the Christoffel symbols of the affine connection calculated using the mutating variant of affine_connection and $R^k_{\ell ij}$ are the components of the Riemann curvature tensor calculated using the mutating variant of riemann_tensor.

The returned StitchedChartSolution{:Jacobi} returns (p, X, Y, dY) at time t, where p is the point on the geodesic, X its velocity, Y the Jacobi field, and dY its covariant derivative.

solve_chart_parallel_transport_ode(
    M::AbstractManifold, a, Xc, A::AbstractAtlas, i0, Yc;
    solver=AutoVern9(Rodas5P()), check_chart_switch_kwargs=NamedTuple(), final_time::Real=1.0,
    kwargs...
)

Parallel transport vector with coordinates Yc along geodesic on a manifold M from point of coordinates a in a chart i0 from an AbstractAtlas A in direction of coordinates Xc in the induced basis.

solve_chart_volume_density(
    M::AbstractManifold, a, Xc, A::AbstractAtlas, i0; kwargs...
)

Compute the volume density of the exponential map in chart coordinates. The coordinates a and Xc are represented in the induced basis of chart i0 from atlas A.

christoffel_symbols_first(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend = default_differential_backend(),
)

Compute the Christoffel symbols of the first kind in local coordinates of basis B. The Christoffel symbols are (in Einstein summation convention)

\[Γ_{ijk} = \frac{1}{2} \Bigl[g_{kj,i} + g_{ik,j} - g_{ij,k}\Bigr],\]

where $g_{ij,k}=\frac{∂}{∂ p^k} g_{ij}$ is the coordinate derivative of the local representation of the metric tensor. The dimensions of the resulting multi-dimensional array are ordered $(i,j,k)$.

christoffel_symbols_second(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend = default_differential_backend(),
)

Compute the Christoffel symbols of the second kind in local coordinates of basis B. For affine connection manifold the Christoffel symbols need to be explicitly implemented while, for a MetricManifold they are computed as (in Einstein summation convention)

\[Γ^{l}_{ij} = g^{kl} Γ_{ijk},\]

where $Γ_{ijk}$ are the Christoffel symbols of the first kind (see christoffel_symbols_first), and $g^{kl}$ is the inverse of the local representation of the metric tensor. The dimensions of the resulting multi-dimensional array are ordered $(l,i,j)$.

christoffel_symbols_second_jacobian(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend = default_differential_backend(),
)

Get partial derivatives of the Christoffel symbols of the second kind for manifold M at p with respect to the coordinates of B, i.e.

\[\frac{∂}{∂ p^l} Γ^{k}_{ij} = Γ^{k}_{ij,l}.\]

The dimensions of the resulting multi-dimensional array are ordered $(k,i,j,l)$.

det_local_metric(M::AbstractManifold, p, B::AbstractBasis)

Return the determinant of local matrix representation of the metric tensor $g$, i.e. of the matrix $G(p)$ representing the metric in the tangent space at $p$ with as a matrix.

See also local_metric

einstein_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend())

Compute the Einstein tensor of the manifold M at the point p, see https://en.wikipedia.org/wiki/Einstein_tensor

gaussian_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend())

Compute the Gaussian curvature of the manifold M at the point p using basis B. This is equal to half of the scalar Ricci curvature, see ricci_curvature.

inverse_local_metric(M::AbstractManifold{𝔽}, p, B::AbstractBasis)

Return the local matrix representation of the inverse metric (cometric) tensor of the tangent space at p on the AbstractManifold M with respect to the AbstractBasis basis B.

The metric tensor (see local_metric) is usually denoted by $G = (g_{ij}) ∈ 𝔽^{d×d}$, where $d$ is the dimension of the manifold.

Then the inverse local metric is denoted by $G^{-1} = g^{ij}$.

local_metric(M::AbstractManifold{𝔽}, p, B::AbstractBasis)

Return the local matrix representation at the point p of the metric tensor $g$ with respect to the AbstractBasis B on the AbstractManifold M. Let $d$denote the dimension of the manifold and $b_1,\ldots,b_d$ the basis vectors. Then the local matrix representation is a matrix $G\in 𝔽^{n×n}$ whose entries are given by $g_{ij} = g_p(b_i,b_j), i,j\in\{1,…,d\}$.

This yields the property for two tangent vectors (using Einstein summation convention) $X = X^ib_i, Y=Y^ib_i \in T_p\mathcal M$ we get $g_p(X, Y) = g_{ij} X^i Y^j$.

local_metric_jacobian(M::AbstractManifold, p, B::AbstractBasis;
    backend::AbstractDiffBackend,
)

Get partial derivatives of the local metric of M at p in basis B with respect to the coordinates of p, $\frac{∂}{∂ p^k} g_{ij} = g_{ij,k}$. The dimensions of the resulting multi-dimensional array are ordered $(i,j,k)$.

log_local_metric_density(M::AbstractManifold, p, B::AbstractBasis)

Return the natural logarithm of the metric density $ρ$ of M at p. The density is given by $ρ = \sqrt{|\det [g_{ij}]|}$ for the metric tensor expressed in basis B.

ricci_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend())

Compute the Ricci scalar curvature of the manifold M at the point p using basis B. The curvature is computed as the trace of the Ricci curvature tensor with respect to the metric, that is $R=g^{ij}R_{ij}$ where $R$ is the scalar Ricci curvature at p, $g^{ij}$ is the inverse local metric (see inverse_local_metric) at p and $R_{ij}$ is the Ricci curvature tensor, see ricci_tensor. Both the tensor and inverse local metric are expressed in local coordinates defined by B, and the formula uses the Einstein summation convention.

ricci_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = default_differential_backend())

Compute the Ricci tensor, also known as the Ricci curvature tensor, of the manifold M at the point p using basis B, see https://en.wikipedia.org/wiki/Ricci_curvature#Introduction_and_local_definition.

solve_exp_ode(
    M::AbstractManifold,
    p,
    X,
    t::Number;
    B::AbstractBasis = DefaultOrthonormalBasis(),
    backend::AbstractDiffBackend = default_differential_backend(),
    solver = AutoVern9(Rodas5P()),
    kwargs...,
)

Approximate the exponential map on the manifold by evaluating the ODE describing the geodesic at 1, assuming the default connection of the given manifold by solving the ordinary differential equation

\[\frac{d^2}{dt^2} p^k + Γ^k_{ij} \frac{d}{dt} p_i \frac{d}{dt} p_j = 0,\]

where $Γ^k_{ij}$ are the Christoffel symbols of the second kind, and the Einstein summation convention is assumed. The argument solver follows the OrdinaryDiffEq conventions. kwargs... specify keyword arguments that will be passed to OrdinaryDiffEq.solve.

Currently, the numerical integration is only accurate when using a single coordinate chart that covers the entire manifold. This excludes coordinates in an embedded space.

using OrdinaryDiffEq

riemann_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend=default_differential_backend())

Compute the Riemann tensor $R^l_{ijk}$, also known as the Riemann curvature tensor, at the point p in local coordinates defined by B. The dimensions of the resulting multi-dimensional array are ordered $(l,i,j,k)$.

The function uses the coordinate expression involving the second Christoffel symbol, see https://en.wikipedia.org/wiki/Riemann_curvature_tensor#Coordinate_expression for details.

See also

christoffel_symbols_second, christoffel_symbols_second_jacobian