# geomstats.numerics package#

## geomstats.numerics.bvp module#

Boundary value problem solvers implementation.

class geomstats.numerics.bvp.ScipySolveBVP(tol=0.001, max_nodes=1000, bc_tol=None, save_result=False)[source]#

Bases: `object`

Wrapper for scipy.integrate.solve_bvp.

integrate(fun, bc, x, y, fun_jac=None, bc_jac=None)[source]#

Solve a boundary value problem for a system of ODEs.

## geomstats.numerics.finite_differences module#

Finite differences machinery.

geomstats.numerics.finite_differences.centered_difference(array, delta=None, axis=-1, endpoints=False)[source]#

Centered difference in a Euclidean space.

Points live in R^m, but are a k-dim embedding (e.g. a curve). Assumes points are in correspondence for cases higher than dim=1.

Parameters:
• array (array-like) – Values of a function.

• delta (float) – Spacing between points.

• axis (int) – Axis in which perform the difference. Must be given backwards.

• endpoints (bool) – If True, endpoints are computed by backward and forward differences, respectively.

Returns:

centered_diff (array-like) – Same shape as array.

geomstats.numerics.finite_differences.forward_difference(array, delta=None, axis=-1)[source]#

Forward difference in a Euclidean space.

Points live in R^m, but are a k-dim embedding (e.g. a curve). Assumes points are in correspondence for cases higher than dim=1.

Parameters:
• array (array-like) – Values of a function.

• delta (float) – Spacing between points.

• axis (int) – Axis in which perform the difference. Must be given backwards.

Returns:

forward_diff (array-like) – Shape in the specified axis reduces by one.

geomstats.numerics.finite_differences.second_centered_difference(array, delta=None, axis=-1)[source]#

Second centered difference in a Euclidean space.

Points live in R^m, but are a k-dim embedding (e.g. a curve). Assumes points are in correspondence for cases higher than dim=1.

Parameters:
• array (array-like) – Values of a function.

• delta (float) – Spacing between points.

• axis (int) – Axis in which perform the difference. Must be given backwards.

Returns:

second_centered_diff (array-like) – Shape in the specified axis reduces by two (endpoints).

## geomstats.numerics.geodesic module#

Geodesic solvers implementation.

class geomstats.numerics.geodesic.ExpODESolver(space, integrator=None)[source]#

Bases: `ExpSolver`

Geodesic initial value problem solver.

Integrate geodesic equation.

Parameters:
• space (Manifold) – Equipped manifold.

• integrator (ODEIVPIntegrator) – Instance of ODEIVP integrator.

exp(tangent_vec, base_point)[source]#

Exponential map.

Parameters:
• tangent_vec (array-like, shape=[…, *space.shape]) – Tangent vector at the base point.

• base_point (array-like, shape=[…, *space.shape]) – Point on the manifold.

Returns:

end_point (array-like, shape=[…, *space.shape]) – Point on the manifold.

geodesic_ivp(tangent_vec, base_point)[source]#

Geodesic curve for initial value problem.

Parameters:
• tangent_vec (array-like, shape=[…, *space.shape]) – Tangent vector at the base point.

• base_point (array-like, shape=[…, *space.shape]) – Point on the manifold.

Returns:

path (callable) – Time parametrized geodesic curve. f(t).

property integrator#

An instance of ODEIVPIntegrator.

class geomstats.numerics.geodesic.ExpSolver(solves_ivp=False)[source]#

Bases: `ABC`

Abstract class for geodesic initial value problem solvers.

Parameters:

solves_ivp (bool) – Informs if solver is able to solve for geodesic at different t.

abstract exp(tangent_vec, base_point)[source]#

Exponential map.

Parameters:
• tangent_vec (array-like, shape=[…, dim]) – Tangent vector at the base point.

• base_point (array-like, shape=[…, dim]) – Point on the manifold.

Returns:

end_point (array-like, shape=[…, dim]) – Point on the manifold.

geodesic_ivp(tangent_vec, base_point)[source]#

Geodesic curve for initial value problem.

Parameters:
• tangent_vec (array-like, shape=[…, dim]) – Tangent vector at the base point.

• base_point (array-like, shape=[…, dim]) – Point on the manifold.

Returns:

path (callable) – Time parametrized geodesic curve. f(t).

class geomstats.numerics.geodesic.LogODESolver(space, n_nodes=10, integrator=None, initialization=None, use_jac=True)[source]#

Bases: `_LogBatchMixins`, `LogSolver`

Geodesic boundary value problem using an ODE solver.

Parameters:
• space (Manifold) – Equipped manifold.

• n_nodes (Number of mesh nodes.)

• integrator (ScipySolveBVP) – Instance of ScipySolveBVP.

• initialization (callable) – Function to provide initial solution. f( point, base_point). Defaults to linear initialization.

geodesic_bvp(point, base_point)[source]#

Geodesic curve for boundary value problem.

Parameters:
• end_point (array-like, shape=[…, dim]) – Point on the manifold.

• base_point (array-like, shape=[…, dim]) – Point on the manifold.

Returns:

path (callable) – Time parametrized geodesic curve. f(t). 0 <= t <= 1.

class geomstats.numerics.geodesic.LogShootingSolver(space, optimizer=None, initialization=None, flatten=True)[source]#

Bases: `object`

Geodesic boundary value problem solver using shooting.

Parameters:
• space (Manifold) – Equipped manifold.

• optimizer (ScipyMinimize) – Instance of ScipyMinimize.

• initialization (callable) – Function to provide initial solution. f(point, base_point). Defaults to linear initialization.

• flatten (bool) – If True, the optimization problem is solved together for all the batch points.

class geomstats.numerics.geodesic.LogSolver(solves_bvp=False)[source]#

Bases: `ABC`

Abstract class for geodesic boundary value problem solvers.

Parameters:

solves_bvp (bool) – Informs if solver is able to solve for geodesic at different t.

geodesic_bvp(point, base_point)[source]#

Geodesic curve for boundary value problem.

Parameters:
• end_point (array-like, shape=[…, dim]) – Point on the manifold.

• base_point (array-like, shape=[…, dim]) – Point on the manifold.

Returns:

path (callable) – Time parametrized geodesic curve. f(t).

abstract log(space, point, base_point)[source]#

Logarithm map.

Parameters:
• end_point (array-like, shape=[…, dim]) – Point on the manifold.

• base_point (array-like, shape=[…, dim]) – Point on the manifold.

Returns:

tangent_vec (array-like, shape=[…, dim]) – Tangent vector at the base point.

class geomstats.numerics.geodesic.PathStraightening(space, path_energy=None, n_nodes=100, optimizer=None, initialization=None)[source]#

Bases: `LogSolver`

Class to solve the geodesic boundary value problem with path-straightening.

Parameters:
• space (Manifold) – Equipped manifold.

• path_energy (callable) – Method to compute Riemannian path energy.

• n_nodes (int) – Number of midpoints.

• optimizer (ScipyMinimize) – An optimizer to solve path energy minimization problem.

• initialization (callable) – A method to get initial guess for optimization.

References

[HSKCB2022]

“Elastic shape analysis of surfaces with second-order Sobolev metrics: a comprehensive numerical framework”. arXiv:2204.04238 [cs.CV], 25 Sep 2022

discrete_geodesic_bvp(point, base_point)[source]#

Solve boundary value problem (BVP).

Given an initial point and an end point, solve the geodesic equation via minimizing the Riemannian path energy.

Parameters:
• point (array-like, shape=[…, *point_shape])

• base_point (array-like, shape=[…, *point_shape])

Returns:

discr_geod_path (array-like, shape=[…, n_times, *point_shape]) – Discrete geodesic.

geodesic_bvp(point, base_point)[source]#

Geodesic curve for boundary value problem.

Parameters:
• end_point (array-like, shape=[…, *point_shape]) – Point on the manifold.

• base_point (array-like, shape=[…, *point_shape]) – Point on the manifold.

Returns:

path (callable) – Time parametrized geodesic curve. f(t).

log(point, base_point)[source]#

Logarithm map.

Parameters:
• end_point (array-like, shape=[…, *point_shape]) – Point on the manifold.

• base_point (array-like, shape=[…, *point_shape]) – Point on the manifold.

Returns:

tangent_vec (array-like, shape=[…, *point_shape]) – Tangent vector at the base point.

## geomstats.numerics.interpolation module#

Interpolation machinery.

class geomstats.numerics.interpolation.Interpolator[source]#

Bases: `ABC`

Abstract class for interpolator.

abstract interpolate(t)[source]#

Interpolate data.

Parameters:

t (array-like, shape=[n_time]) – Interpolation time.

Returns:

point (array-like, shape=[…, n_time, *point_shape])

class geomstats.numerics.interpolation.LinearInterpolator1D(times, data, point_ndim=1)[source]#

Bases: `_LinearInterpolator1D`

A 1D linear interpolator.

Assumes interpolation occurs in the unit interval.

Parameters:
• times (array-like, [n_times]) – Times. Must be sorted.

• data (array-like, […, *point_shape])

• point_ndim (int) – Dimension of point.

class geomstats.numerics.interpolation.UniformUnitIntervalLinearInterpolator(data, point_ndim=1)[source]#

Bases: `_LinearInterpolator1D`

A 1D linear interpolator.

Assumes interpolation occurs in the unit interval and data is uniformly sampled.

Parameters:
• data (array-like, […, *point_shape])

• point_ndim (int) – Dimension of point.

## geomstats.numerics.ivp module#

Initial value problem solvers implementation.

class geomstats.numerics.ivp.GSIVPIntegrator(n_steps=10, step_type='euler', save_result=False)[source]#

In-house ODE integrator.

Parameters:
• n_steps (int) – Number of steps to perform.

• step_type (str) – Type of integration step. Possible values are euler, rk2, rk4.

• save_result (bool) – If True, result is stored after calling integrate or integrate_t.

integrate(force, initial_state, end_time=1.0)[source]#

Integrate force.

Parameters:
• force (callable) – Function to integrate: f(state, t).

• initial_state (array-like, shape=[…, n_vars, *point_shape]) – Initial state.

• end_time (float or None) – Integration end time.

Returns:

result (OdeResult)

property step_type#

Integrator step type.

class geomstats.numerics.ivp.ODEIVPIntegrator(save_result=False, tchosen=False)[source]#

Bases: `ABC`

Abstract class for ode ivp solvers.

Parameters:
• save_result (bool) – If True, result is stored after calling integrate or integrate_t.

• tchosen (bool) – Informs about ability to solve at chosen times. If False, then does not implement integrate_t.

abstract integrate(force, initial_state, end_time)[source]#

Integrate force.

Parameters:
• force (callable) – Function to integrate: f(state, t).

• initial_state (array-like, shape=[…, n_vars, *point_shape]) – Initial state.

• end_time (float or None) – Integration end time.

Returns:

result (OdeResult)

integrate_t(force, initial_state, t_eval)[source]#

Integrate force while choosing evaluating points.

Parameters:
• force (callable) – Function to integrate: f(state, t).

• initial_state (array-like, shape=[…, n_vars, *point_shape]) – Initial state.

• t_eval (array-like) – Times at which to store the computed solution.

Returns:

result (OdeResult)

class geomstats.numerics.ivp.OdeResult[source]#

Bases: `OptimizeResult`

Bunch object (follows scipy).

Its purposes is to homogenize output of different integrators.

get_last_y()[source]#

Get value for last y.

Allows to have y represented as an gs.array or list[gs.array] (latter for cases where they have different shapes).

Assumes last t is the same.

class geomstats.numerics.ivp.ScipySolveIVP(method='RK45', save_result=False, point_ndim=1, **options)[source]#

Wrapper for scipy.integrate.solve_ivp.

Check https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html for additional options.

Parameters:
• method (str) – Integration method.

• save_result (bool) – If True, result is stored after calling integrate or integrate_t.

• point_ndim = int – Dimension of array representing a point in the space.

integrate(force, initial_state, end_time=1.0)[source]#

Integrate force.

Parameters:
• force (callable) – Function to integrate: f(state, t).

• initial_state (array-like, shape=[…, n_vars, *point_shape]) – Initial state.

• end_time (float or None) – Integration end time.

Returns:

result (OdeResult)

integrate_t(force, initial_state, t_eval)[source]#

Integrate force at t_eval points.

Parameters:
• force (callable) – Function to integrate: f(state, t).

• initial_state (array-like, shape=[…, n_vars, *point_shape]) – Initial state.

• t_eval (array-like) – Times at which to store the computed solution.

Returns:

result (OdeResult)

## geomstats.numerics.optimizers module#

Optimizers implementations.

class geomstats.numerics.optimizers.ScipyMinimize(method='L-BFGS-B', jac=None, hess=None, bounds=None, constraints=(), tol=None, callback=None, options=None, save_result=False)[source]#

Bases: `object`

Wrapper for scipy.optimize.minimize.

Only jac differs from scipy: if autodiff, then gs.autodiff.value_and_grad is used to compute the jacobian.

minimize(fun, x0, fun_jac=None, fun_hess=None, hessp=None)[source]#

Minimize objective function.

Parameters:
• fun (callable) – The objective function to be minimized.

• x0 (array-like) – Initial guess.

• fun_jac (callable) – If not None, jac is ignored.

• fun_hess (callable) – If not None, hess is ignored.

• hessp (callable)

## geomstats.numerics.path module#

Discrete-path related machinery.

class geomstats.numerics.path.UniformlySampledDiscretePath(path, interpolator=None, **interpolator_kwargs)[source]#

Bases: `object`

A uniformly-sampled discrete path.

Parameters:
• path (array-like, […, *point_shape])

• interpolator (Interpolator1D)

class geomstats.numerics.path.UniformlySampledPathEnergy(space)[source]#

Bases: `object`

Riemannian path energy of a uniformly-sampled path.

Parameters:

space (Manifold) – Equipped manifold.

energy(path)[source]#

Compute Riemannian path energy.

Parameters:

path (array-like, shape=[…, n_times, *point_shape]) – Piecewise linear path.

Returns:

energy (array-like, shape=[…,]) – Path energy.

energy_per_time(path)[source]#

Compute Riemannian path enery per time.

Parameters:

path (array-like, shape=[…, n_times, *point_shape]) – Piecewise linear path.

Returns:

energy (array-like, shape=[…, n_times - 1,]) – Stepwise path energy.