✨ feat(manifolds): add chord_distance, the straight line through the ambient

docs/api/manifolds.md

@@ -49,6 +49,7 @@ ang = cxm.angle_between(cxc.cart3d, uvec, vvec, at=at)
 - `angle_between`: return the metric angle between two tangent-vector CDicts
 - `norm`: compute the Riemannian norm $\|v\|_g = \sqrt{g_p(v,v)}$ of a tangent vector in a chart. Requires a **positive-definite** metric; raises `NotImplementedError` for an indefinite one (e.g. Minkowski), where the square root would be `nan`
 - `geodesic_distance`: length of the shortest path between two points _along the manifold_. Symmetric and chart-invariant, computed from the manifold's geometry: the straight line in flat space, the great circle on a sphere. A manifold with no closed-form geodesic raises `NotImplementedError` rather than approximating, as does Minkowski, whose indefinite metric admits no distance
+- `chord_distance`: length of the straight line between two points _through the ambient space_ they are embedded in -- the tunnel rather than the surface path. A different measurement from `geodesic_distance`, not an approximation to it: on a sphere of radius $R$ separated by a central angle $\theta$, the geodesic is $R\theta$ and the chord is $2R\sin(\theta/2)$. Needs an embedding, so a manifold that is its own ambient (anything Euclidean) raises `NotImplementedError` and points at `geodesic_distance`
 - `interval`: signed squared interval $\Delta s^2 = \Delta x^\top G\,\Delta x$ of the _coordinate difference_, with the metric taken at the first point. Defined for **every** metric, including indefinite ones, and the causal invariant when Lorentzian. It is not the squared `geodesic_distance` except where the metric is constant along the path, i.e. on a flat manifold in Cartesian coordinates -- flatness alone is not enough, since a curvilinear chart on flat space has a varying metric. The verbs that read its sign need a timelike direction and live in the `coordinax.manifolds.lorentzian` sub-namespace below
 - `pt_embed`: embed intrinsic coordinates into ambient coordinates
 - `pt_project`: project ambient coordinates back to intrinsic chart coordinates

packages/coordinaxs.api/src/coordinaxs/api/manifolds.py

@@ -11,6 +11,7 @@
     "pt_map",
     "norm",
     "geodesic_distance",
+    "chord_distance",
     "interval",
     "causal_character",
     "proper_time",
@@ -95,6 +96,17 @@ def geodesic_distance(*args: Any, **kwargs: Any) -> Any:
     raise NotImplementedError  # pragma: no cover
 
 
+@plum.dispatch.abstract
+def chord_distance(*args: Any, **kwargs: Any) -> Any:
+    """Straight-line distance between two points through their ambient space.
+
+    The chord, as opposed to `geodesic_distance`'s path along the manifold.
+    Defined wherever the manifold carries an embedding; a manifold that is its
+    own ambient space has no distinct chord and is refused.
+    """
+    raise NotImplementedError  # pragma: no cover
+
+
 @plum.dispatch.abstract
 def angle_between(
     chart: Any, uvec: Any, vvec: Any, /, *args: Any, **kwargs: Any

src/coordinax/__init__.py

@@ -36,6 +36,7 @@
     "EmbeddedManifold",
     "CustomAtlas",
     "CustomManifold",
+    "chord_distance",
     "geodesic_distance",
     # frames -- frames
     "noframe",
@@ -132,6 +133,7 @@
     EuclideanManifold,
     FlatMetric,
     Rn,
+    chord_distance,
     embedded_twosphere,
     geodesic_distance,
 )

src/coordinax/_src/manifolds/__init__.py

@@ -3,6 +3,7 @@
 __all__: tuple[str, ...] = ()
 
 from .angle_between import *
+from .chord_distance import *
 from .geodesic_distance import *
 from .guess import *
 from .interval import *

src/coordinax/_src/manifolds/chord_distance.py

@@ -0,0 +1,236 @@
+r"""Dispatch implementations for `coordinaxs.api.manifolds.chord_distance`.
+
+The straight-line distance between two points *through the ambient space* they
+are embedded in -- the tunnel rather than the surface path.
+
+This is a different measurement from
+`~coordinax.manifolds.geodesic_distance`, not an approximation to it. Both are
+exact and symmetric; they answer different questions. On a sphere of radius
+$R$ separated by a central angle $\theta$:
+
+.. math::
+
+    d_{\mathrm{geodesic}} = R\,\theta,
+    \qquad
+    d_{\mathrm{chord}} = 2R \sin(\theta/2).
+
+They agree to first order and diverge as the points separate -- 0.927 against
+0.964 at $\theta = 1$ on the unit sphere, and $2R$ against $\pi R$ at
+antipodes. Which one is wanted depends on the question: a great-circle flight
+path is the geodesic, a line of sight through the body is the chord.
+
+A chord is only defined relative to an embedding, so this is implemented for
+manifolds that carry one. A manifold that is its own ambient space -- anything
+Euclidean -- has no distinct chord, and is directed to `geodesic_distance`
+rather than being given the same number under a second name.
+"""
+
+__all__: tuple[str, ...] = ()
+
+from typing import Any
+
+import plum
+
+import coordinaxs.api.charts as cxcapi
+import coordinaxs.api.manifolds as cxmapi
+from coordinax._src.base import AbstractChart, AbstractManifold
+from coordinax._src.custom_types import OptUSys
+from coordinax._src.embedded.chart import EmbeddedChart
+from coordinax._src.embedded.manifold import EmbeddedManifold
+from coordinax._src.euclidean.manifold import EuclideanManifold
+from coordinax._src.spherical.chart import sph2
+from coordinax._src.spherical.embed import TwoSphereIn3D
+from coordinax._src.spherical.manifold import HyperSphericalManifold
+from coordinaxs.api.custom_types import CDict
+
+
+def _ambient_distance(
+    embedded: EmbeddedChart[Any, Any],
+    chart: AbstractChart,
+    intrinsic_chart: AbstractChart,
+    a: CDict,
+    b: CDict,
+    usys: OptUSys,
+    /,
+) -> Any:
+    """Embed both points, then measure the straight line in the ambient space.
+
+    The ambient manifold is flat, so its `geodesic_distance` *is* the straight
+    line -- there is no separate implementation of the chord itself.
+
+    No fast path here, deliberately. Profiling puts ~99.5% of an eager call in
+    the ambient `geodesic_distance`, and ~0.5% in the embedding steps this
+    could skip; within the former, the cost is a single `pt_map` between charts
+    (~5ms eagerly, ~15us under `jit`). Hand-rolling the sphere's embedding
+    formula here to dodge one `pt_map` would duplicate `TwoSphereIn3D` for a
+    fraction of a percent.
+    """
+    ambient_chart = embedded.ambient
+
+    def embed(p: CDict) -> CDict:
+        intrinsic: CDict = cxcapi.pt_map(p, chart, intrinsic_chart, usys=usys)  # ty: ignore[invalid-assignment]
+        out: CDict = cxmapi.pt_embed(intrinsic, embedded, usys=usys)  # ty: ignore[invalid-assignment]
+        return out
+
+    return cxmapi.geodesic_distance(
+        ambient_chart.M, ambient_chart, embed(a), embed(b), usys=usys
+    )
+
+
+@plum.dispatch
+def chord_distance(
+    chart: AbstractChart, a: CDict, b: CDict, /, *, usys: OptUSys = None
+) -> Any:
+    """Return the ambient straight-line distance, on the chart's manifold.
+
+    >>> import jax.numpy as jnp
+    >>> import unxt as u
+    >>> import coordinax.charts as cxc
+    >>> import coordinax.manifolds as cxm
+
+    A quarter turn along the equator of the unit sphere: the arc is ``pi / 2``,
+    the chord through the interior is ``sqrt(2)``.
+
+    >>> a = {"theta": u.Angle(jnp.pi / 2, "rad"), "phi": u.Angle(0.0, "rad")}
+    >>> b = {"theta": u.Angle(jnp.pi / 2, "rad"), "phi": u.Angle(jnp.pi / 2, "rad")}
+    >>> round(float(cxm.chord_distance(cxc.sph2, a, b)), 6)
+    1.414214
+    >>> round(float(cxm.geodesic_distance(cxc.sph2, a, b).ustrip("rad")), 6)
+    1.570796
+
+    Antipodes are one diameter apart through the middle, half the great-circle
+    distance around the outside:
+
+    >>> n = {"theta": u.Angle(0.0, "rad"), "phi": u.Angle(0.0, "rad")}
+    >>> s = {"theta": u.Angle(jnp.pi, "rad"), "phi": u.Angle(0.0, "rad")}
+    >>> round(float(cxm.chord_distance(cxc.sph2, n, s)), 6)
+    2.0
+
+    """
+    return cxmapi.chord_distance(chart.M, chart, a, b, usys=usys)
+
+
+@plum.dispatch
+def chord_distance(
+    M: HyperSphericalManifold,
+    chart: AbstractChart,
+    a: CDict,
+    b: CDict,
+    /,
+    *,
+    usys: OptUSys = None,
+) -> Any:
+    """Return the chord of the unit hypersphere, through its canonical embedding.
+
+    >>> import jax.numpy as jnp
+    >>> import unxt as u
+    >>> import coordinax.charts as cxc
+    >>> import coordinax.manifolds as cxm
+
+    Any chart on the sphere gives the same answer:
+
+    >>> a = {"lon": u.Angle(0.0, "rad"), "lat": u.Angle(0.0, "rad")}
+    >>> b = {"lon": u.Angle(jnp.pi / 2, "rad"), "lat": u.Angle(0.0, "rad")}
+    >>> round(float(cxm.chord_distance(cxc.lonlat_sph2, a, b)), 6)
+    1.414214
+
+    """
+    if M.ndim != 2:
+        msg = (
+            f"chord_distance is only implemented for the two-sphere; {M} is "
+            f"{M.ndim}-dimensional."
+        )
+        raise NotImplementedError(msg)
+    unit_sphere = EmbeddedChart(TwoSphereIn3D(radius=1.0))
+    return _ambient_distance(unit_sphere, chart, sph2, a, b, usys)
+
+
+@plum.dispatch
+def chord_distance(
+    M: EmbeddedManifold,
+    chart: AbstractChart,
+    a: CDict,
+    b: CDict,
+    /,
+    *,
+    usys: OptUSys = None,
+) -> Any:
+    """Return the chord through the manifold's own ambient space.
+
+    >>> import jax.numpy as jnp
+    >>> import unxt as u
+    >>> import coordinax.charts as cxc
+    >>> import coordinax.manifolds as cxm
+
+    A sphere of radius 2 m: antipodes are one diameter apart.
+
+    >>> M = cxm.EmbeddedManifold(
+    ...     intrinsic=cxm.S2, ambient=cxm.R3,
+    ...     embed_map=cxm.TwoSphereIn3D(radius=u.Q(2.0, "m")),
+    ... )
+    >>> n = {"theta": u.Angle(0.0, "rad"), "phi": u.Angle(0.0, "rad")}
+    >>> s = {"theta": u.Angle(jnp.pi, "rad"), "phi": u.Angle(0.0, "rad")}
+    >>> cxm.chord_distance(M, cxc.sph2, n, s).round(6)
+    Distance(4., 'm')
+
+    """
+    embedded = EmbeddedChart(M.embed_map)
+    return _ambient_distance(embedded, chart, M.embed_map.intrinsic, a, b, usys)
+
+
+@plum.dispatch
+def chord_distance(
+    M: EuclideanManifold,
+    chart: AbstractChart,
+    a: CDict,
+    b: CDict,
+    /,
+    *,
+    usys: OptUSys = None,
+) -> Any:
+    """Refuse: flat space is its own ambient, so the chord is the geodesic.
+
+    Returning the same number under a second name invites the reader to think
+    two things were measured.
+
+    >>> import unxt as u
+    >>> import coordinax.charts as cxc
+    >>> import coordinax.manifolds as cxm
+
+    >>> a = {"x": u.Q(3.0, "m"), "y": u.Q(0.0, "m"), "z": u.Q(0.0, "m")}
+    >>> b = {"x": u.Q(0.0, "m"), "y": u.Q(4.0, "m"), "z": u.Q(0.0, "m")}
+    >>> try: cxm.chord_distance(cxc.cart3d, a, b)
+    ... except NotImplementedError as e: print(e)
+    chord_distance is a measurement through an ambient space, and Rn(3) is its
+    own ambient -- its chord is the straight line, which is what
+    `geodesic_distance` already returns.
+
+    """
+    del chart, a, b, usys
+    msg = (
+        "chord_distance is a measurement through an ambient space, and "
+        f"{M} is its own ambient -- its chord is the straight line, which is "
+        "what `geodesic_distance` already returns."
+    )
+    raise NotImplementedError(msg)
+
+
+@plum.dispatch
+def chord_distance(
+    M: AbstractManifold,
+    chart: AbstractChart,
+    a: CDict,
+    b: CDict,
+    /,
+    *,
+    usys: OptUSys = None,
+) -> Any:
+    """Refuse: without an embedding there is no ambient space to cut through."""
+    del chart, a, b, usys
+    msg = (
+        f"no chord distance is implemented for {M}: a chord is measured "
+        "through an ambient space, and this manifold carries no embedding. "
+        "Wrap it in an `EmbeddedManifold`, or use `geodesic_distance` for the "
+        "distance along the manifold."
+    )
+    raise NotImplementedError(msg)

src/coordinax/_src/manifolds/geodesic_distance.py

@@ -40,7 +40,11 @@
     AbstractMetricField,
     check_metric_is_charts,
 )
-from coordinax._src.charts.d3 import cart3d
+from coordinax._src.charts.d0 import Cart0D
+from coordinax._src.charts.d1 import Cart1D
+from coordinax._src.charts.d2 import Cart2D
+from coordinax._src.charts.d3 import Cart3D, cart3d
+from coordinax._src.charts.dn import CartND
 from coordinax._src.custom_types import OptUSys
 from coordinax._src.embedded.chart import EmbeddedChart
 from coordinax._src.embedded.manifold import EmbeddedManifold
@@ -417,3 +421,57 @@ def geodesic_distance(
     """
     check_metric_is_charts(metric, chart, "geodesic_distance")
     return cxmapi.geodesic_distance(chart.M, chart, a, b, usys=usys)
+
+
+# ===================================================================
+# Fast paths.
+#
+# The packed overloads unpack their operands into component dicts so that a
+# curvilinear chart can be mapped to Cartesian first. In a Cartesian chart
+# there is nothing to map, and plum can skip the unpacking entirely by
+# dispatching on the chart type: 33x for arrays, 3.3x for quantities.
+#
+# There is deliberately no matching `CDict` fast path. The two `pt_map` calls
+# it would skip cost 4.8us of a ~3000us call -- the time is in `norm` and in
+# wrapping the result as a `Distance` -- so it measured as no gain at all, and
+# would have been dispatch surface bought with nothing.
+
+CartesianChart = Cart0D | Cart1D | Cart2D | Cart3D | CartND
+
+
+@plum.dispatch.multi(
+    (CartesianChart, u.AbstractQuantity, u.AbstractQuantity),
+    (CartesianChart, Array, Array),
+)
+def geodesic_distance(
+    chart: CartesianChart, a: Any, b: Any, /, *, usys: OptUSys = None
+) -> Any:
+    """Packed operands in a Cartesian chart: measure without unpacking.
+
+    Registered as two signatures rather than one with a ``Quantity | Array``
+    union: each generic overload it overrides is typed for one of them, and a
+    union would be narrower in the chart but wider in the operands, so plum
+    would rank neither more specific and refuse the call as ambiguous.
+
+    The generic packed overloads split the operands into component dicts so a
+    curvilinear chart can be mapped to Cartesian. Here there is nothing to map,
+    and the components are already on the trailing axis in the right order, so
+    the norm is taken directly.
+
+    >>> import jax.numpy as jnp
+    >>> import unxt as u
+    >>> import coordinax.charts as cxc
+    >>> import coordinax.manifolds as cxm
+
+    >>> a = u.Q(jnp.asarray([3.0, 0.0, 0.0]), "m")
+    >>> b = u.Q(jnp.asarray([0.0, 4.0, 0.0]), "m")
+    >>> cxm.geodesic_distance(cxc.cart3d, a, b).round(2)
+    Distance(5., 'm')
+
+    >>> float(cxm.geodesic_distance(cxc.cart3d, jnp.asarray([3.0, 0.0, 0.0]),
+    ...                             jnp.asarray([0.0, 4.0, 0.0])))
+    5.0
+
+    """
+    del chart, usys
+    return _as_distance(jnp.linalg.norm(b - a, axis=-1))

src/coordinax/manifolds/__init__.py

@@ -11,6 +11,7 @@
     "metric_matrix",
     "metric_representation",
     "norm",
+    "chord_distance",
     "geodesic_distance",
     "interval",
     # Sub-namespaces
@@ -140,6 +141,7 @@
     from coordinaxs.api.charts import pt_map
     from coordinaxs.api.manifolds import (
         angle_between,
+        chord_distance,
         geodesic_distance,
         guess_manifold,
         interval,