Source code for geomstats.geometry.connection

"""Affine connections."""

from scipy.optimize import minimize

import geomstats.backend as gs
import geomstats.errors
import geomstats.vectorization
from geomstats.integrator import integrate


N_STEPS = 10
EPSILON = 1e-3


[docs]class Connection: r"""Class for affine connections. Parameters ---------- dim : int Dimension of the underlying manifold. default_point_type : str, {\'vector\', \'matrix\'} Point type. Optional, default: \'vector\'. default_coords_type : str, {\'intrinsic\', \'extrinsic\', etc} Coordinate type. Optional, default: \'intrinsic\'. """ def __init__( self, dim, default_point_type='vector', default_coords_type='intrinsic'): geomstats.errors.check_integer(dim, 'dim') geomstats.errors.check_parameter_accepted_values( default_point_type, 'default_point_type', ['vector', 'matrix']) self.dim = dim self.default_point_type = default_point_type self.default_coords_type = default_coords_type
[docs] def christoffels(self, base_point): """Christoffel symbols associated with the connection. Parameters ---------- base_point : array-like, shape=[..., dim] Point on the manifold. Returns ------- gamma : array-like, shape=[..., dim, dim, dim] Christoffel symbols, with the covariant index on the first dimension. """ raise NotImplementedError( 'The Christoffel symbols are not implemented.')
[docs] def connection(self, tangent_vec_a, tangent_vec_b, base_point): """Covariant derivative. Connection applied to `tangent_vector_b` in the direction of `tangent_vector_a`, both tangent at `base_point`. Parameters ---------- tangent_vec_a : array-like, shape=[..., dim] Tangent vector at base point. tangent_vec_b : array-like, shape=[..., dim] Tangent vector at base point. base_point : array-like, shape=[..., dim] Point on the manifold. """ raise NotImplementedError( 'connection is not implemented.')
[docs] def geodesic_equation(self, position, velocity): """Compute the geodesic ODE associated with the connection. Parameters ---------- velocity : array-like, shape=[..., dim] Tangent vector at the position. position : array-like, shape=[..., dim] Point on the manifold, the position at which to compute the geodesic ODE. Returns ------- geodesic_ode : array-like, shape=[..., dim] Value of the vector field to be integrated at position. """ gamma = self.christoffels(position) equation = gs.einsum( '...kij,...i->...kj', gamma, velocity) equation = - gs.einsum( '...kj,...j->...k', equation, velocity) return equation
[docs] def exp(self, tangent_vec, base_point, n_steps=N_STEPS, step='euler', point_type=None): """Exponential map associated to the affine connection. Exponential map at base_point of tangent_vec computed by integration of the geodesic equation (initial value problem), using the christoffel symbols Parameters ---------- tangent_vec : array-like, shape=[..., dim] Tangent vector at the base point. base_point : array-like, shape=[..., dim] Point on the manifold. n_steps : int Number of discrete time steps to take in the integration. Optional, default: N_STEPS. step : str, {'euler', 'rk4'} The numerical scheme to use for integration. Optional, default: 'euler'. point_type : str, {'vector', 'matrix'} Type of representation used for points. Optional, default: None. Returns ------- exp : array-like, shape=[..., dim] Point on the manifold. """ initial_state = (base_point, tangent_vec) flow, _ = integrate( self.geodesic_equation, initial_state, n_steps=n_steps, step=step) exp = flow[-1] return exp
[docs] def log(self, point, base_point, n_steps=N_STEPS, step='euler'): """Compute logarithm map associated to the affine connection. Solve the boundary value problem associated to the geodesic equation using the Christoffel symbols and conjugate gradient descent. Parameters ---------- point : array-like, shape=[..., dim] Point on the manifold. base_point : array-like, shape=[..., dim] Point on the manifold. n_steps : int Number of discrete time steps to take in the integration. Optional, default: N_STEPS. step : str, {'euler', 'rk4'} Numerical scheme to use for integration. Optional, default: 'euler'. Returns ------- tangent_vec : array-like, shape=[..., dim] Tangent vector at the base point. """ def objective(velocity): """Define the objective function.""" velocity = gs.array(velocity) velocity = gs.cast(velocity, dtype=base_point.dtype) velocity = gs.reshape(velocity, base_point.shape) delta = self.exp(velocity, base_point, n_steps, step) - point loss = 1. / self.dim * gs.sum(delta ** 2, axis=-1) return gs.sum(loss) objective_with_grad = gs.autograd.value_and_grad(objective) tangent_vec = gs.random.rand(*gs.flatten(base_point).shape) res = minimize( objective_with_grad, tangent_vec, method='L-BFGS-B', jac=True, options={'disp': False, 'maxiter': 25}) tangent_vec = gs.array(res.x) tangent_vec = gs.reshape(tangent_vec, base_point.shape) return tangent_vec
def _pole_ladder_step(self, base_point, next_point, base_shoot, return_geodesics=False): """Compute one Pole Ladder step. One step of pole ladder scheme [LP2013a]_ using the geodesic to transport along as main_geodesic of the parallelogram. Parameters ---------- base_point : array-like, shape=[..., dim] Point on the manifold, from which to transport. next_point : array-like, shape=[..., dim] Point on the manifold, to transport to. base_shoot : array-like, shape=[..., dim] Point on the manifold, end point of the geodesics starting from the base point with initial speed to be transported. return_geodesics : bool, optional (defaults to False) Whether to return the geodesics of the construction. Returns ------- next_step : dict of array-like and callable with following keys: next_tangent_vec : array-like, shape=[..., dim] Tangent vector at end point. end_point : array-like, shape=[..., dim] Point on the manifold, closes the geodesic parallelogram of the construction. geodesics : list of callable, len=3 (only if `return_geodesics=True`) Three geodesics of the construction. References ---------- .. [LP2013a] Marco Lorenzi, Xavier Pennec. Efficient Parallel Transport of Deformations in Time Series of Images: from Schild's to Pole Ladder. Journal of Mathematical Imaging and Vision, Springer Verlag, 2013,50 (1-2), pp.5-17. ⟨10.1007/s10851-013-0470-3⟩ """ mid_tangent_vector_to_shoot = 1. / 2. * self.log( base_point=base_point, point=next_point) mid_point = self.exp( base_point=base_point, tangent_vec=mid_tangent_vector_to_shoot) tangent_vector_to_shoot = - self.log( base_point=mid_point, point=base_shoot) end_shoot = self.exp( base_point=mid_point, tangent_vec=tangent_vector_to_shoot) next_tangent_vec = - self.log( base_point=next_point, point=end_shoot) end_point = self.exp( base_point=next_point, tangent_vec=next_tangent_vec) geodesics = [] if return_geodesics: main_geodesic = self.geodesic( initial_point=base_point, end_point=next_point) diagonal = self.geodesic( initial_point=mid_point, end_point=base_shoot) final_geodesic = self.geodesic( initial_point=next_point, end_point=end_shoot) geodesics = [main_geodesic, diagonal, final_geodesic] return {'next_tangent_vec': next_tangent_vec, 'geodesics': geodesics, 'end_point': end_point} def _schild_ladder_step(self, base_point, next_point, base_shoot, return_geodesics=False): """Compute one Schild's Ladder step. One step of the Schild's ladder scheme [LP2013a]_ using the geodesic to transport along as one side of the parallelogram. Parameters ---------- base_point : array-like, shape=[..., dim] Point on the manifold, from which to transport. next_point : array-like, shape=[..., dim] Point on the manifold, to transport to. base_shoot : array-like, shape=[..., dim] Point on the manifold, end point of the geodesics starting from the base point with initial speed to be transported. return_geodesics : bool Whether to return points computed along each geodesic of the construction. Optional, default: False. Returns ------- transported_tangent_vector : array-like, shape=[..., dim] Tangent vector at end point. end_point : array-like, shape=[..., dim] Point on the manifold, closes the geodesic parallelogram of the construction. References ---------- .. [LP2013a] Marco Lorenzi, Xavier Pennec. Efficient Parallel Transport of Deformations in Time Series of Images: from Schild's to Pole Ladder. Journal of Mathematical Imaging and Vision, Springer Verlag, 2013,50 (1-2), pp.5-17. ⟨10.1007/s10851-013-0470-3⟩ """ mid_tangent_vector_to_shoot = 1. / 2. * self.log( base_point=base_shoot, point=next_point) mid_point = self.exp( base_point=base_shoot, tangent_vec=mid_tangent_vector_to_shoot) tangent_vector_to_shoot = - self.log( base_point=mid_point, point=base_point) end_shoot = self.exp( base_point=mid_point, tangent_vec=tangent_vector_to_shoot) next_tangent_vec = self.log( base_point=next_point, point=end_shoot) geodesics = [] if return_geodesics: main_geodesic = self.geodesic( initial_point=base_point, end_point=next_point) diagonal = self.geodesic( initial_point=base_point, end_point=end_shoot) second_diagonal = self.geodesic( initial_point=base_shoot, end_point=next_point) final_geodesic = self.geodesic( initial_point=next_point, end_point=end_shoot) geodesics = [ main_geodesic, diagonal, second_diagonal, final_geodesic] return {'next_tangent_vec': next_tangent_vec, 'geodesics': geodesics, 'end_point': end_shoot}
[docs] def ladder_parallel_transport( self, tangent_vec_a, tangent_vec_b, base_point, n_steps=1, step='pole', **single_step_kwargs): """Approximate parallel transport using the pole ladder scheme. Approximate Parallel transport using either the pole ladder or the Schild's ladder scheme [LP2013b]_. Pole ladder is exact in symmetric spaces [GJSP2019]_ while Schild's ladder is a first order approximation. Both schemes are available any affine connection manifolds whose exponential and logarithm maps are implemented. `tangent_vec_a` is transported along the geodesic starting at the base_point with initial tangent vector `tangent_vec_b`. Parameters ---------- tangent_vec_a : array-like, shape=[..., dim] Tangent vector at base point to transport. tangent_vec_b : array-like, shape=[..., dim] Tangent vector at base point, initial speed of the geodesic along which to transport. base_point : array-like, shape=[..., dim] Point on the manifold, initial position of the geodesic along which to transport. n_steps : int The number of pole ladder steps. Optional, default: 1. step : str, {'pole', 'schild'} The scheme to use for the construction of the ladder at each step. Optoinal, default: 'pole'. **single_step_kwargs : keyword arguments for the step functions Returns ------- ladder : dict of array-like and callable with following keys transported_tangent_vector : array-like, shape=[..., dim] Approximation of the parallel transport of tangent vector a. trajectory : list of list of callable, len=n_steps List of lists containing the geodesics of the construction, only if `return_geodesics=True` in the step function. The geodesics are methods of the class connection. References ---------- .. [LP2013b] Marco Lorenzi, Xavier Pennec. Efficient Parallel Transpor of Deformations in Time Series of Images: from Schild's to Pole Ladder.Journal of Mathematical Imaging and Vision, Springer Verlag, 2013, 50 (1-2), pp.5-17. ⟨10.1007/s10851-013-0470-3⟩ .. [GJSP2019] N. Guigui, Shuman Jia, Maxime Sermesant, Xavier Pennec. Symmetric Algorithmic Components for Shape Analysis with Diffeomorphisms. GSI 2019, Aug 2019, Toulouse, France. pp.10. ⟨hal-02148832⟩ """ current_point = gs.copy(base_point) next_tangent_vec = gs.copy(tangent_vec_a) / n_steps methods = {'pole': self._pole_ladder_step, 'schild': self._schild_ladder_step} single_step = methods[step] base_shoot = self.exp( base_point=current_point, tangent_vec=next_tangent_vec) trajectory = [] for i_point in range(n_steps): frac_tangent_vector_b = (i_point + 1) / n_steps * tangent_vec_b next_point = self.exp( base_point=base_point, tangent_vec=frac_tangent_vector_b) next_step = single_step( base_point=current_point, next_point=next_point, base_shoot=base_shoot, **single_step_kwargs) current_point = next_point base_shoot = next_step['end_point'] trajectory.append(next_step['geodesics']) transported_tangent_vec = n_steps * next_step['next_tangent_vec'] return {'transported_tangent_vec': transported_tangent_vec, 'trajectory': trajectory}
[docs] def riemannian_curvature(self, base_point): """Compute Riemannian curvature tensor associated with the connection. Parameters ---------- base_point: array-like, shape=[..., dim] Point on the manifold. """ raise NotImplementedError( 'The Riemannian curvature tensor is not implemented.')
[docs] def geodesic(self, initial_point, end_point=None, initial_tangent_vec=None, point_type=None): """Generate parameterized function for the geodesic curve. Geodesic curve defined by either: - an initial point and an initial tangent vector, - an initial point and an end point. Parameters ---------- initial_point : array-like, shape=[..., dim] Point on the manifold, initial point of the geodesic. end_point : array-like, shape=[..., dim], optional Point on the manifold, end point of the geodesic. If None, an initial tangent vector must be given. initial_tangent_vec : array-like, shape=[..., dim], Tangent vector at base point, the initial speed of the geodesics. Optional, default: None. If None, an end point must be given and a logarithm is computed. point_type : str, {'vector', 'matrix'} Point type. Optional, default: 'vector'. Returns ------- path : callable Time parameterized geodesic curve. If a batch of initial conditions is passed, the output array's first dimension represents time, and the second corresponds to the different initial conditions. """ if point_type is None: point_type = self.default_point_type geomstats.errors.check_parameter_accepted_values( point_type, 'point_type', ['vector', 'matrix']) if end_point is None and initial_tangent_vec is None: raise ValueError('Specify an end point or an initial tangent ' 'vector to define the geodesic.') if end_point is not None: shooting_tangent_vec = self.log( point=end_point, base_point=initial_point) if initial_tangent_vec is not None: if not gs.allclose(shooting_tangent_vec, initial_tangent_vec): raise RuntimeError( 'The shooting tangent vector is too' ' far from the input initial tangent vector.') initial_tangent_vec = shooting_tangent_vec if point_type == 'vector': initial_point = gs.to_ndarray(initial_point, to_ndim=2) initial_tangent_vec = gs.to_ndarray( initial_tangent_vec, to_ndim=2) else: initial_point = gs.to_ndarray(initial_point, to_ndim=3) initial_tangent_vec = gs.to_ndarray( initial_tangent_vec, to_ndim=3) n_initial_conditions = initial_tangent_vec.shape[0] def path(t): """Generate parameterized function for geodesic curve. Parameters ---------- t : array-like, shape=[n_points,] Times at which to compute points of the geodesics. """ t = gs.array(t, gs.float32) t = gs.to_ndarray(t, to_ndim=1) if point_type == 'vector': tangent_vecs = gs.einsum( 'i,...k->...ik', t, initial_tangent_vec) else: tangent_vecs = gs.einsum( 'i,...kl->...ikl', t, initial_tangent_vec) points_at_time_t = [ self.exp(tv, pt) for tv, pt in zip(tangent_vecs, initial_point)] points_at_time_t = gs.stack(points_at_time_t, axis=1) return points_at_time_t[:, 0] if n_initial_conditions == 1 else \ points_at_time_t return path
[docs] def torsion(self, base_point): """Compute torsion tensor associated with the connection. Parameters ---------- base_point: array-like, shape=[..., dim] Point on the manifold. """ raise NotImplementedError( 'The torsion tensor is not implemented.')