ptd.splinelib

Fit smooth spline.

class potamides.splinelib.CostFn(*args, **kwargs)

Bases: Protocol

Protocol for a cost function.

Examples

>>> import jax.numpy as jnp
>>> isinstance(jnp.sum, CostFn)
True

potamides.splinelib.acceleration(spline: Interpolator1D, gamma: Real[Array, ''], /) → Real[Array, 'F']

Compute the acceleration vector \(\vec{a} = d^2\vec{x}/d\gamma^2\).

This is the second derivative of the spline position with respect to the curve parameter \(\gamma\). Equivalently, it’s the derivative of the tangent vector:

\[\vec{a}(\gamma) = \frac{d^2\vec{x}}{d\gamma^2} = \frac{d}{d\gamma} \left(\frac{d\vec{x}}{d\gamma}\right)\]

Parameters:

• spline (interpax.Interpolator1D) – A twice-differentiable 1D spline for \(\vec{x}(\gamma)\).

• gamma (float) – The scalar parameter value at which to compute the acceleration.

Returns:

The acceleration vector \(\frac{d^2\vec{x}}{d\gamma^2}\) of length F, where F is the spatial dimension of the spline.

Return type:

Array[float, (F,)]

Examples

>>> import jax
>>> import jax.numpy as jnp
>>> import interpax
>>> from potamides.splinelib import acceleration

>>> gamma = jnp.linspace(0, 2 * jnp.pi, 10_000)
>>> xy = 2 * jnp.stack([jnp.cos(gamma), jnp.sin(gamma)], axis=-1)
>>> spline = interpax.Interpolator1D(gamma, xy, method="cubic2")

>>> gamma = jnp.array([0, jnp.pi / 2, jnp.pi])
>>> acc = jax.vmap(acceleration, in_axes=(None, 0))(spline, gamma)
>>> print(acc.round(5))
[[-2.  0.]
 [ 0. -2.]
 [ 2.  0.]]

potamides.splinelib.arc_length(spline: Interpolator1D, gamma0: Real[Array, ''] | float | int = -1, gamma1: Real[Array, ''] | float | int = 1, *, method: Literal['p2p', 'quad', 'ode'] = 'p2p', method_kw: dict[str, Any] | None = None) → Real[Array, '']

Return the arc-length of the track.

\[s(\gamma_0, \gamma_1) = \int_{\gamma_0}^{\gamma_1} \left\| \frac{d\mathbf{x}(\gamma)}{d\gamma} \right\| \, d\gamma\]

Computing the arc-length requires computing an integral over the norm of the tangent vector. This can be done using many different methods. We provide three options, specified by the method parameter.

Parameters:

• gamma0 – The starting / ending gamma value between which to compute the arc-length. The default is [-1, 1], which is the full range of gamma for the track.

• gamma1 – The starting / ending gamma value between which to compute the arc-length. The default is [-1, 1], which is the full range of gamma for the track.

• method –

  The method to use for computing the arc-length. Options are “p2p”, “quad”, or “ode”. The default is “p2p”.

  • ”p2p”: point-to-point distance. This method computes the distance

    between each pair of points along the track and sums them up. Accuracy is limited by the 1e5 points used.

  • ”quad”: quadrature. This method uses fixed quadrature to compute

    the integral. It is the default method. It also uses 1e5 points.

  • ”ode”: ODE integration. This method uses ODE integration to

    compute the integral.

See also

None

This method auto-vectorizes to support arbitrarily shaped gamma inputs.

None

This method computes the distance between each pair of points along the track and sums them up. Accuracy is limited by the number of points used. This can be selected by setting method=”p2p”.

None

This method uses fixed quadrature to compute the integral. It is also limited in accuracy by the number of points used. This can be selected by setting method=”quad”.

None

This method uses ODE integration to compute the integral. This can be selected by setting method=”ode”.

Examples

>>> import jax
>>> import jax.numpy as jnp
>>> import interpax
>>> import potamides.splinelib as splib

>>> gamma = jnp.linspace(0, 2 * jnp.pi, 10_000)
>>> xy = 2 * jnp.stack([jnp.cos(gamma), jnp.sin(gamma)], axis=-1)
>>> spline = interpax.Interpolator1D(gamma, xy, method="cubic2")

>>> s = splib.arc_length(spline, 0, 2 * jnp.pi) / jnp.pi
>>> print(s.round(5))
4.0

>>> s = splib.arc_length(spline, 0, 2 * jnp.pi, method="quad") / jnp.pi
>>> print(s.round(4))
4.0

>>> s = splib.arc_length(spline, 0, 2 * jnp.pi, method="ode") / jnp.pi
>>> print(s.round(5))
4.0

potamides.splinelib.arc_length_odeint(spline: Interpolator1D, gamma0: Real[Array, ''] | float | int = -1, gamma1: Real[Array, ''] | float | int = 1, *, rtol: float = 1.4e-08, atol: float = 1.4e-08, mxstep: float = inf, hmax: float = inf) → Real[Array, '']

Compute the arc-length using ODE integration.

Parameters:

• spline – The spline interpolator.

• gamma0 – The starting / ending gamma value between which to compute the arc-length. The default is [-1, 1], which is the full range of gamma for the track.

• gamma1 – The starting / ending gamma value between which to compute the arc-length. The default is [-1, 1], which is the full range of gamma for the track.

• rtol – The relative and absolute tolerances for the ODE solver. The default is 1.4e-8.

• atol – The relative and absolute tolerances for the ODE solver. The default is 1.4e-8.

• mxstep – The maximum number of steps and maximum step size for the ODE solver. The default is inf.

• hmax – The maximum number of steps and maximum step size for the ODE solver. The default is inf.

See also

None

This method auto-vectorizes to support arbitrarily shaped gamma inputs. It also allows for other methods of computing the arc-length.

None

This method computes the distance between each pair of points along the track and sums them up. Accuracy is limited by the number of points used.

None

This method uses fixed quadrature to compute the integral. It is also limited in accuracy by the number of points used.

Examples

>>> import jax
>>> import jax.numpy as jnp
>>> import interpax
>>> import potamides.splinelib as splib

>>> gamma = jnp.linspace(0, 2 * jnp.pi, 10_000)
>>> xy = 2 * jnp.stack([jnp.cos(gamma), jnp.sin(gamma)], axis=-1)
>>> spline = interpax.Interpolator1D(gamma, xy, method="cubic2")

>>> s = splib.arc_length_odeint(spline, 0, 2 * jnp.pi) / jnp.pi
>>> print(s.round(5))
4.0

potamides.splinelib.arc_length_p2p(spline: Interpolator1D, gamma0: Real[Array, ''] | float | int = -1, gamma1: Real[Array, ''] | float | int = 1, *, num: int = 100000) → Real[Array, '']

Compute the arc-length using a point-to-point distance approximation.

Parameters:

• spline – The spline interpolator.

• gamma0 – The starting / ending gamma value between which to compute the arc-length. The default is [-1, 1], which is the recommended range of gamma for the track.

• gamma1 – The starting / ending gamma value between which to compute the arc-length. The default is [-1, 1], which is the recommended range of gamma for the track.

• num – The integer number of points to use for the quadrature. The default is 100,000.

See also

None

This method auto-vectorizes to support arbitrarily shaped gamma inputs. It also allows for other methods of computing the arc-length.

None

This method uses fixed quadrature to compute the integral. It is also limited in accuracy by the number of points used.

None

This method uses ODE integration to compute the integral.

Examples

>>> import jax
>>> import jax.numpy as jnp
>>> import interpax
>>> import potamides.splinelib as splib

>>> gamma = jnp.linspace(0, 2 * jnp.pi, 10_000)
>>> xy = 2 * jnp.stack([jnp.cos(gamma), jnp.sin(gamma)], axis=-1)
>>> spline = interpax.Interpolator1D(gamma, xy, method="cubic2")

>>> s = splib.arc_length_p2p(spline, 0, 2 * jnp.pi) / jnp.pi
>>> print(s.round(5))
4.0

potamides.splinelib.arc_length_quadrature(spline: Interpolator1D, gamma0: Real[Array, ''] | float | int = -1, gamma1: Real[Array, ''] | float | int = 1, *, num: int = 100000) → Real[Array, '']

Compute the arc-length using fixed quadrature.

Parameters:

• spline – The spline interpolator.

• gamma0 – The starting / ending gamma value between which to compute the arc-length. The default is [-1, 1], which is the full range of gamma for the track.

• gamma1 – The starting / ending gamma value between which to compute the arc-length. The default is [-1, 1], which is the full range of gamma for the track.

• num – The number of points to use for the quadrature. The default is 100,000.

See also

None

This method auto-vectorizes to support arbitrarily shaped gamma inputs. It also allows for other methods of computing the arc-length.

None

This method computes the distance between each pair of points along the track and sums them up. Accuracy is limited by the number of points used.

None

This method uses ODE integration to compute the integral.

Examples

>>> import jax
>>> import jax.numpy as jnp
>>> import interpax
>>> import potamides.splinelib as splib

>>> gamma = jnp.linspace(0, 2 * jnp.pi, 10_000)
>>> xy = 2 * jnp.stack([jnp.cos(gamma), jnp.sin(gamma)], axis=-1)
>>> spline = interpax.Interpolator1D(gamma, xy, method="cubic2")

>>> s = splib.arc_length_quadrature(spline, 0, 2 * jnp.pi, num=500_000) / jnp.pi
>>> print(s.round(4))  # harder to make accurate
4.0

potamides.splinelib.concavity_change_cost_fn(knots: Real[Array, 'N 2'], gamma: Real[Array, 'N'], /, data_gamma: Real[Array, 'data'], scale: float = 100.0, num_points: int = 1000) → Real[Array, '']

Cost function to penalize changes in signed curvature for 2D curves.

The integrand of the cost function is the derivative of the arctangent of the signed curvature multiplied by a large number \(\lambda\).

\[\left( \frac{d}{ds} \atan\left(\lambda \kappa_{\text{signed}}(s)\right) \right)^2\]

where \(\kappa_{\text{signed}}(s)\) is the signed curvature at \(s\) and \(\lambda\) is a large number that controls the width of the smoothing. The \(\atan\) function differentiably mimics the non-differentiable \(\text{sign}\) function. The cost is the integral over the arc-length.

Parameters:

• knots – Output values of spline at gamma – e.g. x, y values.

• gamma – The gamma values at which the spline is anchored. There are N of these, one per knots. These are fixed while the knots are optimized. These should be from the same distribution as data_gamma.

• data_gamma – gamma of the target data.

• scale – Inverse Smoothing width.

• num_points – The number of points to use to compute the cost function. This should be large enough to capture the curvature of the spline. Default is 1_000.

potamides.splinelib.curvature(spline: Interpolator1D, gamma: Real[Array, ''], /) → Real[Array, 'F']

Return the curvature vector at a given position along the stream.

This method computes the curvature by taking the ratio of the gamma derivative of the unit tangent vector to the derivative of the arc-length with respect to gamma. In other words, if

\[\frac{d\hat{T}}{d\gamma} = \frac{ds}{d\gamma} \frac{d\hat{T}}{ds}\]

and since the curvature vector is defined as

\[\frac{d\hat{T}}{ds} = \kappa \hat{N}\]

where \(\kappa\) is the curvature and \(\hat{N}\) the unit normal vector, then dividing \(\frac{d\hat{T}}{d\gamma}\) by \(\frac{ds}{d\gamma}\) yields

\[\kappa \hat{N} = \frac{d\hat{T}/d\gamma}{ds/d\gamma}\]

Here, \(\frac{d\hat{T}}{d\gamma}\) (computed by dThat_dgamma) describes how the direction of the tangent changes with respect to the affine parameter \(\gamma\), and \(\frac{ds}{d\gamma}\) (obtained from state_speed) represents the state speed (i.e. the rate of change of arc-length with respect to \(\gamma\)).

This formulation assumes that \(\gamma\) is chosen to be proportional to the arc-length of the track.

Parameters:

• spline – The spline interpolator.

• gamma – The gamma value at which to evaluate the curvature.

See also

None

This method auto-vectorizes to support arbitrarily shaped gamma inputs.

Examples

>>> import jax
>>> import jax.numpy as jnp
>>> import interpax
>>> import potamides.splinelib as splib

>>> gamma = jnp.linspace(0, 2 * jnp.pi, 10_000)
>>> xy = 2 * jnp.stack([jnp.cos(gamma), jnp.sin(gamma)], axis=-1)
>>> spline = interpax.Interpolator1D(gamma, xy, method="cubic2")

>>> gamma = jnp.array([0, jnp.pi / 2, jnp.pi])
>>> kappa_vec = jax.vmap(splib.curvature, (None, 0))(spline, gamma)
>>> print(kappa_vec.round(5))
[[-0.5  0. ]
 [ 0.  -0.5]
 [ 0.5  0. ]]

potamides.splinelib.data_distance_cost_fn(knots: Real[Array, 'N 2'], gamma: Real[Array, 'N'], /, data_gamma: Real[Array, 'data'], data_y: Real[Array, 'data'], *, sigmas: Real[Array, 'N'] | float = 1.0) → Real[Array, '']

Cost function to minimize that compares data to spline fit.

\[\text{cost} = \sum_i \left( \frac{y_i - f(\gamma_i)}{\sigma_i} \right)^2\]

where \(y_i\) is the target data, \(f(\gamma_i)\) is the spline evaluated at \(\gamma_i\), and \(\sigma_i\) is the uncertainty on \(y_i\).

Parameters:

• knots – Output values of spline at gamma – e.g. x, y values.

• gamma – The gamma values at which the spline is anchored. There are N of these, one per knots. These are fixed while the knots are optimized.

• data_gamma – gamma of the target data.

• data_y – The target data. This is the data that the spline is trying to fit.

• sigmas – The uncertainty on each datum in data_y.

potamides.splinelib.default_cost_fn(knots: Real[Array, 'N 2'], gamma: Real[Array, 'N'], data_gamma: Real[Array, 'data'], data_y: Real[Array, 'data'], /, *, sigmas: float = 1.0, data_weight: float = 1000.0, concavity_weight: float = 0.0, concavity_scale: float = 100.0) → Real[Array, '']

Cost function to minimize that compares data to spline fit.

Parameters:

• knots – Output values of spline at gamma – e.g. x or y values. This is the parameter to be optimized to minimize the cost function.

• gamma – The gamma values at which the spline is anchored. There are N of these, one per knots. These are fixed while the knots are optimized. These should be from the same distribution as data_gamma.

• data_gamma – gamma of the target data.

• data_y – The target data. This is the data that the spline is trying to fit.

• sigmas – The uncertainty on each datum in data_y.

• concavity_weight – The weight of the curvature penalization term. This should be tuned based on the desired smoothness of the spline. Default is 0.0.

• concavity_scale – The scale of the smoothing for the curvature penalization term. This should be tuned based on the desired smoothness of the spline. Default is 1e2.

potamides.splinelib.interpax_PPoly_from_scipy_UnivariateSpline(scipy_spl: UnivariateSpline, /) → PPoly

Convert a scipy.interpolate.UnivariateSpline to an interpax.PPoly.

Notes

This conversion relies on scipy’s internal _eval_args attribute, which may not be available for all scipy spline types. The function works by: 1. Converting the scipy UnivariateSpline to a scipy PPoly using from_spline 2. Extracting coefficients and breakpoints from the scipy PPoly 3. Constructing an equivalent interpax PPoly instance.

Examples

>>> import numpy as np
>>> from scipy.interpolate import UnivariateSpline
>>> import interpax

>>> x = np.linspace(0, 2*np.pi, 10)
>>> y = np.sin(x)

>>> scipy_spline = UnivariateSpline(x, y, s=0)
>>> interpax_ppoly = interpax_PPoly_from_scipy_UnivariateSpline(scipy_spline)

>>> x_test = np.linspace(0, 2*np.pi, 50)
>>> scipy_vals = scipy_spline(x_test)
>>> interpax_vals = interpax_ppoly(x_test)
>>> np.allclose(scipy_vals, interpax_vals)
True

potamides.splinelib.kappa(spline: Interpolator1D, gamma: Real[Array, ''], /) → Real[Array, '']

Return the curvature magnitude at a given position along the stream.

Parameters:

• spline – The spline interpolator.

• gamma – The gamma value at which to evaluate the curvature.

See also

None

This method auto-vectorizes to support arbitrarily shaped gamma inputs.

Examples

>>> import jax
>>> import jax.numpy as jnp
>>> import interpax
>>> import potamides.splinelib as splib

>>> gamma = jnp.linspace(0, 2 * jnp.pi, 10_000)
>>> xy = 2 * jnp.stack([jnp.cos(gamma), jnp.sin(gamma)], axis=-1)
>>> spline = interpax.Interpolator1D(gamma, xy, method="cubic2")

>>> gamma = jnp.array([0, jnp.pi / 2, jnp.pi])
>>> kappa_val = jax.vmap(splib.kappa, (None, 0))(spline, gamma)
>>> print(kappa_val.round(5))  # circles have constant curvature
[0.5 0.5 0.5]

potamides.splinelib.make_gamma_from_data(data: Real[Array, 'data 2'], /) → Real[Array, 'data-1']

Return \(\gamma\), the normalized arc-length of the data.

\[\gamma = 2\frac{s}{L} - 1 , \in [-1, 1]\]

where \(s\) is the arc-length at \(\gamma\) and \(L\) is the total arc-length.

Gamma is constructed approximately using a point-to-point approximation (the function point_to_point_arclength).

Notes

This is guaranteed to be monotonically non-decreasing since the point-to-point arc-length is always non-negative. However, this is not guaranteed to be monotonically increasing since adjacent data points can have 0 distance. See make_increasing_gamma_from_data for a function that trims the data such that gamma is monotonically increasing.

Examples

>>> import jax.numpy as jnp

>>> data = jnp.array([[0, 0], [1, 0], [1, 2], [0, 2]])
>>> make_gamma_from_data(data)
Array([-1.        ,  0.33333333,  1.        ], dtype=float64)

potamides.splinelib.make_increasing_gamma_from_data(data: Real[Array, 'data 2'], /) → tuple[Real[Array, 'data-1'], Real[Array, 'data-1 2']]

Return the trimmed data and gamma, the normalized arc-length.

\[\gamma = 2\frac{s}{L} - 1 , \in [-1, 1]\]

where \(s\) is the arc-length at \(\gamma\) and \(L\) is the total arc-length.

Gamma is constructed approximately using a point-to-point (P2P) approximation (the function point_to_point_arclength). Using the P2P arc-length is not guaranteed to be monotonically increasing since adjacent data points can have 0 distance. This function then trims the data such that gamma is monotonically increasing, keeping the first point of any plateau.

Returns:

• gamma (Array[real, (N-1,)]) – Monotonically increasing normalized arc-length.

• data_trimmed (Array[real, (N-1, 2)]) – The data, with points where gamma is non-increasing trimmed out.

Examples

>>> import jax.numpy as jnp

>>> data = jnp.array([[0, 0], [1, 0], [1, 2], [1, 2], [0, 2]])
>>> gamma, data2 = make_increasing_gamma_from_data(data)
>>> gamma, data2
(Array([-1.        ,  0.33333333,  1.        ], dtype=float64),
 Array([[0.5, 0. ],
       [1. , 1. ],
       [0.5, 2. ]], dtype=float64))

Note that the second point [1, 2] was removed since it was a repeat, resulting in a “plateau” in gamma. Then the point-to-point mean was returned as the new data.

potamides.splinelib.new_gamma_knots_from_spline(spline: Interpolator1D, /, *, nknots: int) → tuple[Real[Array, '{nknots}'], Real[Array, '{nknots} 2']]

Define new gamma (and knots) from an existing spline.

When the knots of a spline are changed the arc-length of the spline changes as well. It is often useful to define a new gamma that is the normalized arc-length of the spline. This function takes a spline and returns a new gamma (and corresponding knots) that is the normalized arc-length of the spline so that a new spline can be created with the new gamma (and knots).

Parameters:

• spline – The spline to use to define the new gamma.

• nknots – The number of knots to use in the new spline.

Returns:

• gamma_new – The new gamma values. One is at -1 and one is at 1. The rest are evenly spaced in between.

• points_new – The new points of the spline at the new gamma values.

potamides.splinelib.optimize_spline_knots(cost_fn: ~potamides._src.splinelib.opt.CostFn, /, init_knots: ~jaxtyping.Real[Array, 'N 2'], init_gamma: ~jaxtyping.Real[Array, 'N'], cost_args: tuple[~typing.Any, ...], *, cost_kwargs: ~xmmutablemap._core.ImmutableMap[str, ~typing.Any] | tuple[tuple[str, ~typing.Any], ...] | None = None, optimizer: ~optax._src.base.GradientTransformation = (<function chain.<locals>.init_fn>, <function chain.<locals>.update_fn>), nsteps: int = 10000, fixed_mask: tuple[bool, ...] | None = None) → Real[Array, 'N 2']

Optimize spline knots to fit data.

Warning

If you use this function to change the locations of the knots then this changes the arc-length of the spline. This can be problematic if gamma is the normalized arc-length of the data. If you change the knots then you should also change gamma accordingly. The easiest way to do this is to:

1. evaluate the optimized spline on a dense array of old gamma values

2. call make_gamma_from_data on the new data to define a new gamma,

3. create a new spline with the new gamma. This spline will have the same shape as the optimized spline but with the new gamma values.

Parameters:

• cost_fn – The cost function.

• init_knots – starting outputs of splines at init_gamma.

• init_gamma – anchor points for spline. median gamma in chunk.

• cost_args – Additional positional arguments to pass to the cost function. For example, cost_fn can be default_cost_fn which takes data_gamma and data_y as arguments.

• cost_kwargs – Additional keyword arguments to pass to the cost function. E.g. data_distance_cost_fn can take ‘sigmas’ and concavity_change_cost_fn can take concavity_weight and concavity_scale. default_cost_fn can take both. JAX treats these as static.

• fixed_mask – A mask that indicates which knots are fixed. If None then all knots are free to be optimized. If a mask is provided then the knots at the indices where the mask is True are fixed and the knots at the indices where the mask is False are free to be optimized. The fixed knots are not optimized and are not included in the cost function. This is useful if you want to optimize only a subset of the knots while keeping the others fixed. The mask should be the same length as init_knots. The fixed knots are reconstructed in the final output.

• optimizer – The optimizer to use. Defaults to Adam with a learning rate of 1e-3.

• nsteps – The number of optimization steps to take. Defaults to 10_000.

potamides.splinelib.point_to_point_arclength(data: Real[Array, 'data 2'], /) → Real[Array, 'data-1']

Return a P2P approximation of the arc-length.

The data should be sorted, otherwise this doesn’t make a lot of sense.

Examples

>>> import jax.numpy as jnp

>>> data = jnp.array([[0, 0], [1, 0], [1, 2], [0, 2]])
>>> point_to_point_arclength(data)
Array([1., 3., 4.], dtype=float64)

potamides.splinelib.point_to_point_distance(data: Real[Array, 'data 2'], /) → Real[Array, 'data-1']

Return the distance between points in data.

The data should be sorted, otherwise this doesn’t make a lot of sense.

Examples

>>> import jax.numpy as jnp

>>> data = jnp.array([[0, 0], [1, 0], [1, 2], [0, 2]])
>>> point_to_point_distance(data)
Array([1., 2., 1.], dtype=float64)

potamides.splinelib.position(spline: Interpolator1D, gamma: Real[Array, 'N'], /) → Real[Array, 'N F']

Compute \(\vec{x}(\gamma)\) for spline \(\vec{x}\) at gamma.

This is the Cartesian position vector at the given parameter values gamma. The output is an array with shape (N, F), where N is the number of input gamma values and F is the number of dimensions of the spline.

Parameters:

• spline – The spline interpolator.

• gamma – The gamma values at which to evaluate the spline. This can be a scalar or an array of shape (N,).

Returns:

The position vector \(\vec{x}(\gamma)\) at the specified positions. The shape is (N, F), where N is the number of input gamma values and F is the number of dimensions of the spline.

Return type:

Array[real, (N, F)]

See also

None

This method auto-vectorizes to support arbitrarily shaped gamma inputs.

Examples

>>> import jax
>>> import jax.numpy as jnp
>>> import interpax
>>> import potamides.splinelib as splib

>>> gamma = jnp.linspace(0, 2 * jnp.pi, 10_000)
>>> xy = jnp.stack([jnp.cos(gamma), jnp.sin(gamma)], axis=-1)
>>> spline = interpax.Interpolator1D(gamma, xy, method="cubic2")

>>> gamma = jnp.array([0, jnp.pi / 2, jnp.pi])
>>> pos = jax.vmap(splib.position, (None, 0))(spline, gamma)
>>> print(pos.round(4))
[[ 1.  0.]
 [ 0.  1.]
 [-1.  0.]]

potamides.splinelib.principle_unit_normal(spline: Interpolator1D, gamma: Real[Array, ''], /) → Real[Array, 'F']

Return the unit normal vector \(\hat{N}(\gamma)\) along the spline.

The unit normal vector is defined as the projection of the acceleration vector onto the plane orthogonal to the unit tangent vector, divided by its norm:

\[\hat{N}(\gamma) = \frac{d\hat{T}/d\gamma}{|d\hat{T}/d\gamma|}\]

where \(\hat{T}(\gamma)\) is the unit tangent vector at \(\gamma\) and \(\vec{a}(\gamma)\) is the acceleration vector at \(\gamma\). This function is scalar. To compute the unit normal vector at multiple positions, use jax.vmap.

Parameters:

• spline – The spline interpolator.

• gamma – The scalar gamma value at which to evaluate the unit normal vector. To evaluate the unit normal vector at multiple positions, use jax.vmap.

Returns:

The unit normal vector at the specified position. The shape is (F,), where F is the number of dimensions of the spline.

Return type:

Array[real, (F,)]

See also

None

This method auto-vectorizes to support arbitrarily shaped gamma inputs.

Examples

>>> import jax
>>> import jax.numpy as jnp
>>> import interpax
>>> import potamides.splinelib as splib

>>> gamma = jnp.linspace(0, 2 * jnp.pi, 10_000)
>>> xy = 2 * jnp.stack([jnp.cos(gamma), jnp.sin(gamma)], axis=-1)
>>> spline = interpax.Interpolator1D(gamma, xy, method="cubic2")

>>> gamma = jnp.array([0, jnp.pi / 2, jnp.pi])
>>> N_hat = jax.vmap(splib.principle_unit_normal, in_axes=(None, 0))(spline, gamma)
>>> print(N_hat.round(5))
[[-1.  0.]
 [ 0. -1.]
 [ 1.  0.]]

potamides.splinelib.reduce_point_density(gamma: ~jaxtyping.Real[Array, 'N'], data: ~jaxtyping.Real[Array, 'N 2'], *, num_splits: int, reduce_fn: ~potamides._src.splinelib.opt.ReduceFn = <PjitFunction of <function median>>) → tuple[Real[Array, '{num_splits + 2}'], Real[Array, '{num_splits + 2} 2']]

Split and reduce gamma, data into num_splits blocks, keeping ends.

A dataset representing the points along a stream’s track can have problematic small changes in curvature. If we reduce the number of points that represents the curve then it necessarily forces a greater degree of smoothness. Combining this with optimize_spline_knots can produce a spline curve that better represents the smooth stream track.

Parameters:

• gamma – The gamma values at which the spline is anchored.

• data – The data points of the spline.

• num_splits – The number of splits to make in the data. The spline will be reduced to num_splits + 2 points.

• reduce_fn – The function to use to reduce the data within each chunk to a single point. Defaults to jnp.median.

Examples

>>> import jax.numpy as jnp

>>> gamma = jnp.array([-1, 0, 0.5, 1])
>>> data = jnp.array([[0, 0], [1, 0], [1, 2], [0, 2]])

>>> gamma2, data2 = reduce_point_density(gamma, data, num_splits=1)
>>> gamma2
Array([-1.  ,  0.25,  1.  ], dtype=float64)
>>> data2
Array([[0. , 0. ],
       [0.5, 1. ],
       [0. , 2. ]], dtype=float64)

potamides.splinelib.speed(spline: Interpolator1D, gamma: Real[Array, ''], /) → Real[Array, 'F']

Return the speed in gamma of the track at a given position.

This is the norm of the tangent vector at the given position.

\[v(\gamma) = \| \frac{d\vec{x}(\gamma)}{d\gamma} \| = \|\vec{T}(\gamma) \|\]

An important note is that this is also the differential arc-length!

\[s = \int_{\gamma_0}^{\gamma} \|\frac{\vec{x}}{d\gamma'}\| d\gamma'.\]

Thus, the arc-length element is:

\[\frac{ds}{d\gamma} = \|\frac{\vec{x}}{d\gamma'}\|\]

If we are working in 2D in the flat-sky approximation for extragalactic streams, then it is recommended for \(\gamma\) to be proportional to the arc-length with \(\gamma \in [-1, 1] = \frac{2s}{L} - 1\), we have

\[\frac{ds}{d\gamma} = \frac{L}{2}\]

where \(L\) is the total arc-length of the stream.

Parameters:

• spline – The spline interpolator.

• gamma – The scalar gamma value at which to evaluate the spline. To evaluate the speed at multiple gammas, use jax.vmap.

Returns:

The speed at the specified position. The shape is (F,), where F is the number of dimensions of the spline.

Return type:

Array[real, (F,)]

See also

None

This method auto-vectorizes to support arbitrarily shaped gamma inputs.

Examples

>>> import jax
>>> import jax.numpy as jnp
>>> import interpax
>>> import potamides.splinelib as splib

>>> gamma = jnp.linspace(0, 2 * jnp.pi, 10_000)
>>> xy = 2 * jnp.stack([jnp.cos(gamma), jnp.sin(gamma)], axis=-1)
>>> spline = interpax.Interpolator1D(gamma, xy, method="cubic2")

>>> gamma = jnp.array([0, jnp.pi / 2, jnp.pi])
>>> speed = jax.vmap(splib.speed, (None, 0))(spline, gamma)
>>> print(speed)  # see 2 in xy
[2. 2. 2.]

potamides.splinelib.spherical_position(spline: Interpolator1D, gamma: Real[Array, ''], /) → Real[Array, 'F']

Compute \(|\vec{f}(gamma)|\) for spline \(\vec{f}\) at gamma.

This is the spherical coordinate at the given parameter values gamma. The output is an array with shape (N, F), where N is the number of input gamma values and F is the number of dimensions of the spline. The 0th feature is the radius, and the remaining features are the angular coordinates.

The radius is defined as the Euclidean norm of the position vector:

\[r(\gamma) = \left\| \vec{x}(\gamma) \right\|\]

The angular coordinates are computed recursively using:

\[\phi_i(\gamma) = \arctan2( R_{i+1}, x_i ), \quad \text{for } i = 0, \dots, F-2\]

where \(R_{i+1} = \sqrt{\sum_{j=i+1}^{F-1} x_j^2}\) is the partial radius from the i-th coordinate to the last coordinate.

The last angular coordinate is special-cased as it only depends on the last two coordinates:

\[\phi_{F-1}(\gamma) = \arctan2\left(x_F, x_{F-1}\right)\]

For more details, see https://en.wikipedia.org/wiki/N-sphere.

Parameters:

• spline – The spline interpolator.

• gamma – The scalar gamma value at which to evaluate the spline.

Returns:

The spherical coordinates at the specified position. The shape is (F,), where F is the number of dimensions of the spline. The 0th feature is the radius, and the remaining features are the angular coordinates. To evaluate the spherical coordinates at multiple positions, use jax.vmap.

Return type:

Array[real, (F,)]

See also

None

This method auto-vectorizes to support arbitrarily shaped gamma inputs.

Examples

>>> import jax
>>> import jax.numpy as jnp
>>> import interpax
>>> import potamides.splinelib as splib

>>> gamma = jnp.linspace(0, 2 * jnp.pi, 10_000)
>>> xy = 2 * jnp.stack([jnp.cos(gamma), jnp.sin(gamma)], axis=-1)
>>> spline = interpax.Interpolator1D(gamma, xy, method="cubic2")

>>> gamma = jnp.array([0, jnp.pi / 2, jnp.pi])
>>> r = jax.vmap(splib.spherical_position, (None, 0))(spline, gamma)
>>> print(r.round(4))
[[2.     0.    ]
 [2.     1.5708]
 [2.     3.1416]]

potamides.splinelib.tangent(spline: Interpolator1D, gamma: Real[Array, ''], /) → Real[Array, 'F']

Compute the tangent vector at a given position along the stream.

The tangent vector is defined as:

\[T(\gamma) = \frac{d\vec{x}}{d\gamma}\]

This function is scalar. To compute the unit tangent vector at multiple positions, use jax.vmap.

Parameters:

• spline – The spline interpolator.

• gamma – The scalar gamma value at which to evaluate the spline. To evaluate the tangent vector at multiple positions, use jax.vmap.

Returns:

The tangent vector at the specified position.

Return type:

Array[real, (F,)]

See also

None

This method auto-vectorizes to support arbitrarily shaped gamma inputs.

Examples

>>> import jax
>>> import jax.numpy as jnp
>>> import interpax
>>> import potamides.splinelib as splib

>>> gamma = jnp.linspace(0, 2 * jnp.pi, 10_000)
>>> xy = 2 * jnp.stack([jnp.cos(gamma), jnp.sin(gamma)], axis=-1)
>>> spline = interpax.Interpolator1D(gamma, xy, method="cubic2")

>>> gamma = jnp.array([0, jnp.pi / 2, jnp.pi])
>>> tangents = jax.vmap(splib.tangent, (None, 0))(spline, gamma)
>>> print(tangents.round(2))
[[ 0.  2.]
 [-2.  0.]
 [ 0. -2.]]