Riemannian Manifold API

Manifolds and their geodesic operations, plus the epsilon-net construction that approximates a manifold by a quantum graph.

kmeanssa_ng.riemannian_manifold

Riemannian Manifold module for kmeanssa-ng.

BolzaSurface

Bases: RiemannianManifold

The Bolza surface as a quotient-aware Riemannian manifold.

Points are (..., 2) real (Re, Im) coordinates in the fundamental octagon of the Poincaré disk. The geodesic operations are intrinsic to the surface: :meth:log takes each target to its nearest copy under the Fuchsian group, and :meth:exp folds every step back into the fundamental domain. With :meth:random_uniform, :meth:embed, :meth:log, :meth:exp and :meth:norm all quotient-aware, the intrinsic :class:RepulsionNet builds a correct net, whereas an embedding-based (extrinsic) search would fail across the identified boundary.

Parameters:

Name	Type	Description	Default
ball		Group ball used by the quotient metric (defaults to the cached geometric ball, exact over the fundamental domain).	None

Source code in kmeanssa_ng/riemannian_manifold/bolza.py

class BolzaSurface(RiemannianManifold):
    """The Bolza surface as a quotient-aware Riemannian manifold.

    Points are ``(..., 2)`` real ``(Re, Im)`` coordinates in the fundamental
    octagon of the Poincaré disk. The geodesic operations are intrinsic to the
    *surface*: :meth:`log` takes each target to its nearest copy under the Fuchsian
    group, and :meth:`exp` folds every step back into the fundamental domain. With
    :meth:`random_uniform`, :meth:`embed`, :meth:`log`, :meth:`exp` and :meth:`norm`
    all quotient-aware, the intrinsic :class:`RepulsionNet` builds a correct net,
    whereas an embedding-based (extrinsic) search would fail across the identified
    boundary.

    Args:
        ball: Group ball used by the quotient metric (defaults to the cached
            geometric ball, exact over the fundamental domain).
    """

    def __init__(self, ball=None) -> None:
        super().__init__(_PoincareDisk())
        self.ball = default_ball() if ball is None else ball

    # -- quotient-aware geodesic operations on (..., 2) real arrays -------- #
    def exp(self, base_point: np.ndarray, tangent_vec: np.ndarray) -> np.ndarray:
        base = _to_complex(base_point)
        tangent = _to_complex(tangent_vec)
        moved = exp_map(base, tangent)
        return _to_real(fold_to_domain(moved))

    def log(self, base_point: np.ndarray, point: np.ndarray) -> np.ndarray:
        base = _to_complex(base_point)
        target = nearest_copy(base, _to_complex(point), self.ball)
        return _to_real(log_map(base, target))

    def norm(self, base_point: np.ndarray, tangent_vec: np.ndarray) -> np.ndarray:
        return tangent_norm(_to_complex(base_point), _to_complex(tangent_vec))

    def to_tangent(self, base_point: np.ndarray, ambient_vec: np.ndarray) -> np.ndarray:
        # The disk chart is 2-dimensional: every ambient vector is already tangent.
        return np.asarray(ambient_vec, dtype=float)

    def random_tangent(
        self, base_point: np.ndarray, rng: np.random.Generator
    ) -> np.ndarray:
        """Metric-isotropic tangent draw on the conformal disk chart.

        The hyperbolic metric is conformal with factor 2/(1 - |z|^2), so
        scaling an ambient Gaussian by (1 - |z|^2)/2 yields a vector whose
        components are standard normal in an orthonormal frame of the metric
        -- the same Brownian law at every point of the surface.
        """
        gaussian = rng.standard_normal(self.shape)
        z = _to_complex(np.asarray(base_point, dtype=float))
        return gaussian * (1.0 - np.abs(z) ** 2) / 2.0

    def embed(self, points: np.ndarray) -> np.ndarray:
        """Not available: the quotient has no faithful Euclidean embedding.

        ``embed`` exists so that an *extrinsic* neighbour search (a KD-tree on
        ambient coordinates) can accelerate strategies on isometrically embedded
        manifolds. The Bolza surface is a quotient: points glued across the octagon
        boundary are close on the surface yet far apart in any disk-coordinate
        image, so no such embedding is faithful. Returning the disk coordinates
        anyway would make extrinsic neighbour search silently wrong, so this raises
        instead. Use the intrinsic strategies (:class:`RepulsionNet`, which searches
        by geodesic distance), not the extrinsic ones.
        """
        raise NotImplementedError(
            "BolzaSurface has no faithful Euclidean embedding for neighbour search; "
            "use the intrinsic epsilon-net strategies (e.g. RepulsionNet), not the "
            "extrinsic ones (RepulsionNetExtrinsicSpeedup, build_epsilon_net_graph)."
        )

    def random_uniform(
        self, n: int, random_state: int | np.random.Generator | None = None
    ) -> np.ndarray:
        """``(n, 2)`` points uniform by hyperbolic area in the fundamental octagon."""
        rng = np.random.default_rng(random_state)
        return _to_real(sample_fundamental_domain(n, rng))

    # -- point-level metric (Space interface), quotient distance ---------- #
    def distance(self, point1, point2) -> float:
        d = quotient_distance(
            _to_complex(point1.coordinates), _to_complex(point2.coordinates), self.ball
        )
        return float(np.asarray(d).item())

embed(points)

Not available: the quotient has no faithful Euclidean embedding.

embed exists so that an extrinsic neighbour search (a KD-tree on ambient coordinates) can accelerate strategies on isometrically embedded manifolds. The Bolza surface is a quotient: points glued across the octagon boundary are close on the surface yet far apart in any disk-coordinate image, so no such embedding is faithful. Returning the disk coordinates anyway would make extrinsic neighbour search silently wrong, so this raises instead. Use the intrinsic strategies (:class:RepulsionNet, which searches by geodesic distance), not the extrinsic ones.

Source code in kmeanssa_ng/riemannian_manifold/bolza.py

def embed(self, points: np.ndarray) -> np.ndarray:
    """Not available: the quotient has no faithful Euclidean embedding.

    ``embed`` exists so that an *extrinsic* neighbour search (a KD-tree on
    ambient coordinates) can accelerate strategies on isometrically embedded
    manifolds. The Bolza surface is a quotient: points glued across the octagon
    boundary are close on the surface yet far apart in any disk-coordinate
    image, so no such embedding is faithful. Returning the disk coordinates
    anyway would make extrinsic neighbour search silently wrong, so this raises
    instead. Use the intrinsic strategies (:class:`RepulsionNet`, which searches
    by geodesic distance), not the extrinsic ones.
    """
    raise NotImplementedError(
        "BolzaSurface has no faithful Euclidean embedding for neighbour search; "
        "use the intrinsic epsilon-net strategies (e.g. RepulsionNet), not the "
        "extrinsic ones (RepulsionNetExtrinsicSpeedup, build_epsilon_net_graph)."
    )

random_tangent(base_point, rng)

Metric-isotropic tangent draw on the conformal disk chart.

The hyperbolic metric is conformal with factor 2/(1 - |z|^2), so scaling an ambient Gaussian by (1 - |z|^2)/2 yields a vector whose components are standard normal in an orthonormal frame of the metric -- the same Brownian law at every point of the surface.

Source code in kmeanssa_ng/riemannian_manifold/bolza.py

def random_tangent(
    self, base_point: np.ndarray, rng: np.random.Generator
) -> np.ndarray:
    """Metric-isotropic tangent draw on the conformal disk chart.

    The hyperbolic metric is conformal with factor 2/(1 - |z|^2), so
    scaling an ambient Gaussian by (1 - |z|^2)/2 yields a vector whose
    components are standard normal in an orthonormal frame of the metric
    -- the same Brownian law at every point of the surface.
    """
    gaussian = rng.standard_normal(self.shape)
    z = _to_complex(np.asarray(base_point, dtype=float))
    return gaussian * (1.0 - np.abs(z) ** 2) / 2.0

random_uniform(n, random_state=None)

(n, 2) points uniform by hyperbolic area in the fundamental octagon.

Source code in kmeanssa_ng/riemannian_manifold/bolza.py

def random_uniform(
    self, n: int, random_state: int | np.random.Generator | None = None
) -> np.ndarray:
    """``(n, 2)`` points uniform by hyperbolic area in the fundamental octagon."""
    rng = np.random.default_rng(random_state)
    return _to_real(sample_fundamental_domain(n, rng))

EpsilonNetStrategy

Bases: ABC

Abstract base class for epsilon-net placement strategies.

Source code in kmeanssa_ng/riemannian_manifold/epsilon_net.py

class EpsilonNetStrategy(ABC):
    """Abstract base class for epsilon-net placement strategies."""

    def __init__(self, random_state: int | np.random.Generator | None = None):
        self.random_state = random_state

    @abstractmethod
    def build(self, manifold: RiemannianManifold, n: int) -> np.ndarray:
        """Return ``n`` points on ``manifold`` forming a quasi-uniform net."""
        raise NotImplementedError

build(manifold, n) abstractmethod

Return n points on manifold forming a quasi-uniform net.

Source code in kmeanssa_ng/riemannian_manifold/epsilon_net.py

@abstractmethod
def build(self, manifold: RiemannianManifold, n: int) -> np.ndarray:
    """Return ``n`` points on ``manifold`` forming a quasi-uniform net."""
    raise NotImplementedError

FibonacciNet

Bases: EpsilonNetStrategy

Deterministic Fibonacci lattice on the 2-sphere.

A near-optimal, low-discrepancy net available only for S^2; random_state is ignored since the construction is deterministic.

Source code in kmeanssa_ng/riemannian_manifold/epsilon_net.py

class FibonacciNet(EpsilonNetStrategy):
    """Deterministic Fibonacci lattice on the 2-sphere.

    A near-optimal, low-discrepancy net available only for S^2; ``random_state``
    is ignored since the construction is deterministic.
    """

    def build(self, manifold: RiemannianManifold, n: int) -> np.ndarray:
        if not (manifold.is_sphere and manifold.dim == 2):
            raise ValueError("FibonacciNet is only defined for the 2-sphere S^2.")
        i = np.arange(n)
        golden = (1 + 5**0.5) / 2
        z = 1 - 2 * (i + 0.5) / n
        r = np.sqrt(1 - z * z)
        theta = 2 * np.pi * i / golden
        return np.c_[r * np.cos(theta), r * np.sin(theta), z]

KarcherFrechetMean

Bases: LloydUpdateStrategy

Update strategy that computes the new center as the Fréchet mean (Karcher mean) of the points in the cluster.

The mean is computed by the intrinsic Karcher iteration mean <- exp_mean(average of log_mean(x_i)), driven entirely by the space's own exp/log maps. This works on every space of the package -- including quotient surfaces like Bolza, whose log picks the nearest copy under the group action, which a chart-level estimator (e.g. geomstats' FrechetMean on the underlying manifold) would silently ignore.

The iteration is deterministic, fast and locally exact, but it needs the space to provide exp/log and it descends to the nearest critical point. For a stochastic, globally-minded alternative that works on any metric space of the package (graphs included), see SimulatedAnnealingFrechetMean.

Parameters:

Name	Type	Description	Default
max_iter	int	Maximum number of Karcher iterations.	64
tol	float	Stop when the Riemannian norm of the mean update step falls below this threshold.	1e-09

Source code in kmeanssa_ng/riemannian_manifold/update.py

class KarcherFrechetMean(LloydUpdateStrategy):
    """Update strategy that computes the new center as the Fréchet mean
    (Karcher mean) of the points in the cluster.

    The mean is computed by the intrinsic Karcher iteration
    ``mean <- exp_mean(average of log_mean(x_i))``, driven entirely by the
    space's own ``exp``/``log`` maps. This works on every space of the
    package -- including quotient surfaces like Bolza, whose ``log`` picks
    the nearest copy under the group action, which a chart-level estimator
    (e.g. geomstats' ``FrechetMean`` on the underlying manifold) would
    silently ignore.

    The iteration is deterministic, fast and locally exact, but it needs the
    space to provide ``exp``/``log`` and it descends to the nearest critical
    point. For a stochastic, globally-minded alternative that works on any
    metric space of the package (graphs included), see
    ``SimulatedAnnealingFrechetMean``.

    Args:
        max_iter: Maximum number of Karcher iterations.
        tol: Stop when the Riemannian norm of the mean update step falls
            below this threshold.
    """

    def __init__(self, max_iter: int = 64, tol: float = 1e-9):
        self.max_iter = max_iter
        self.tol = tol

    def update(
        self, points: list[RiemannianPoint], space: "RiemannianManifold"
    ) -> "RiemannianCenter":
        """Compute the new center for a given cluster of points.

        Args:
            points: A list of points belonging to a single cluster.
            space: The Riemannian manifold in which the points and center exist.

        Returns:
            The new center for the cluster.
        """
        if not points:
            return None

        # Extract coordinates from RiemannianPoint objects
        points_coords = []
        for p in points:
            if not isinstance(p, space.get_point_type()):
                raise TypeError("All points must be RiemannianPoint instances.")
            points_coords.append(p.coordinates)

        coords = np.asarray(points_coords, dtype=float)

        # Karcher iteration: follow the mean of the log directions until the
        # step vanishes (exp/log are batched over leading axes).
        mean = coords[0].copy()
        for _ in range(self.max_iter):
            logs = np.asarray(space.log(mean, coords))
            step = logs.mean(axis=0)
            if float(np.asarray(space.norm(mean, step))) < self.tol:
                break
            mean = np.asarray(space.exp(mean, step))

        # Create a RiemannianCenter from the mean coordinates
        return space.center_from_point(space.get_point_type()(space, mean))

update(points, space)

Compute the new center for a given cluster of points.

Parameters:

Name	Type	Description	Default
points	list[RiemannianPoint]	A list of points belonging to a single cluster.	required
space	'RiemannianManifold'	The Riemannian manifold in which the points and center exist.	required

Returns:

Type	Description
'RiemannianCenter'	The new center for the cluster.

Source code in kmeanssa_ng/riemannian_manifold/update.py

def update(
    self, points: list[RiemannianPoint], space: "RiemannianManifold"
) -> "RiemannianCenter":
    """Compute the new center for a given cluster of points.

    Args:
        points: A list of points belonging to a single cluster.
        space: The Riemannian manifold in which the points and center exist.

    Returns:
        The new center for the cluster.
    """
    if not points:
        return None

    # Extract coordinates from RiemannianPoint objects
    points_coords = []
    for p in points:
        if not isinstance(p, space.get_point_type()):
            raise TypeError("All points must be RiemannianPoint instances.")
        points_coords.append(p.coordinates)

    coords = np.asarray(points_coords, dtype=float)

    # Karcher iteration: follow the mean of the log directions until the
    # step vanishes (exp/log are batched over leading axes).
    mean = coords[0].copy()
    for _ in range(self.max_iter):
        logs = np.asarray(space.log(mean, coords))
        step = logs.mean(axis=0)
        if float(np.asarray(space.norm(mean, step))) < self.tol:
            break
        mean = np.asarray(space.exp(mean, step))

    # Create a RiemannianCenter from the mean coordinates
    return space.center_from_point(space.get_point_type()(space, mean))

RepulsionNet

Bases: EpsilonNetStrategy

Uniform init relaxed by truncated Riesz repulsion (Riemannian gradient flow).

Fully intrinsic: each point is pushed away from its k nearest neighbours -- found by geodesic distance -- along geodesics, with a step that is a shrinking fraction of the local spacing, until the configuration freezes into a near-regular net. Being driven only by the manifold's intrinsic log, exp and norm, it is correct on any compact Riemannian manifold, including quotient spaces (e.g. the Bolza surface) where an ambient embedding would misrepresent neighbours across identified boundaries.

Neighbour search is O(n^2) in the point count. When the manifold is isometrically embedded and ambient (chordal) proximity matches geodesic proximity -- as for S^2 in R^3 -- :class:RepulsionNetExtrinsicSpeedup replaces it by an O(n log n) KD-tree search on the embedding.

Parameters:

Name	Type	Description	Default
k	int	Number of nearest neighbours each point repels.	12
n_iter	int	Number of relaxation steps.	600
riesz_s	float	Exponent of the repulsion kernel (force ~ 1 / d^(s+1)).	1.0
cooling	float	Larger values cool more slowly (step ~ 1 / (1 + t / cooling)).	300.0
random_state	int | Generator | None	Seed or Generator for the uniform initialisation.	None

Source code in kmeanssa_ng/riemannian_manifold/epsilon_net.py

class RepulsionNet(EpsilonNetStrategy):
    """Uniform init relaxed by truncated Riesz repulsion (Riemannian gradient flow).

    Fully *intrinsic*: each point is pushed away from its ``k`` nearest neighbours
    -- found by geodesic distance -- along geodesics, with a step that is a
    shrinking fraction of the local spacing, until the configuration freezes into
    a near-regular net. Being driven only by the manifold's intrinsic ``log``,
    ``exp`` and ``norm``, it is correct on any compact Riemannian manifold,
    including quotient spaces (e.g. the Bolza surface) where an ambient embedding
    would misrepresent neighbours across identified boundaries.

    Neighbour search is O(n^2) in the point count. When the manifold is
    isometrically embedded and ambient (chordal) proximity matches geodesic
    proximity -- as for S^2 in R^3 -- :class:`RepulsionNetExtrinsicSpeedup`
    replaces it by an O(n log n) KD-tree search on the embedding.

    Args:
        k: Number of nearest neighbours each point repels.
        n_iter: Number of relaxation steps.
        riesz_s: Exponent of the repulsion kernel (force ~ 1 / d^(s+1)).
        cooling: Larger values cool more slowly (step ~ 1 / (1 + t / cooling)).
        random_state: Seed or Generator for the uniform initialisation.
    """

    # Cap on the transient (chunk * n, dim) buffer of the neighbour search.
    _NEIGHBOUR_BUFFER = 2_000_000

    def __init__(
        self,
        k: int = 12,
        n_iter: int = 600,
        riesz_s: float = 1.0,
        cooling: float = 300.0,
        random_state: int | np.random.Generator | None = None,
    ):
        super().__init__(random_state)
        self.k = k
        self.n_iter = n_iter
        self.riesz_s = riesz_s
        self.cooling = cooling

    def _neighbours(self, manifold: RiemannianManifold, X: np.ndarray) -> np.ndarray:
        """Indices of each point's ``k`` nearest neighbours by *geodesic* distance.

        Distances come from the manifold's own ``log``/``norm`` (so the metric --
        quotient-aware or not -- is whatever the space defines), evaluated in row
        chunks to bound memory. O(n^2); self is excluded. Returns ``(n, k)``.
        """
        n = X.shape[0]
        idx = np.empty((n, self.k), dtype=int)
        chunk = max(1, self._NEIGHBOUR_BUFFER // n)
        for lo in range(0, n, chunk):
            hi = min(lo + chunk, n)
            m = hi - lo
            base = np.repeat(X[lo:hi], n, axis=0)
            tiled = np.tile(X, (m,) + (1,) * (X.ndim - 1))
            d = manifold.norm(base, manifold.log(base, tiled)).reshape(m, n)
            d[np.arange(m), np.arange(lo, hi)] = np.inf  # exclude self
            idx[lo:hi] = np.argpartition(d, self.k - 1, axis=1)[:, : self.k]
        return idx

    def build(self, manifold: RiemannianManifold, n: int) -> np.ndarray:
        # NOTE: on non-compact domains (e.g. raw hyperbolic space) unbounded
        # repulsion pushes points to the boundary and leaves interior voids.
        # Supporting those needs a bounded region with a soft confining potential;
        # this is orthogonal to the intrinsic/extrinsic neighbour choice. Compact
        # spaces (spheres, the Bolza surface) need no confinement.
        rng = np.random.default_rng(self.random_state)
        X = manifold.random_uniform(n, rng)
        dim = X.shape[1]
        s = self.riesz_s
        for t in range(self.n_iter):
            frac = 0.25 / (1 + t / self.cooling)
            idx = self._neighbours(manifold, X)
            base = np.repeat(X, self.k, axis=0)
            neighbours = X[idx].reshape(-1, dim)
            tangent = manifold.log(base, neighbours)  # base -> neighbour
            dist = np.maximum(manifold.norm(base, tangent), 1e-6)
            unit = tangent / dist[:, None]
            weight = 1.0 / dist ** (s + 1)  # Riesz repulsion
            force = -(weight[:, None] * unit).reshape(n, self.k, dim).sum(axis=1)
            local = np.median(dist.reshape(n, self.k), axis=1, keepdims=True)
            f_norm = np.maximum(manifold.norm(X, force), 1e-12)[:, None]
            # Move at most a shrinking fraction of the local spacing (stable).
            step = np.minimum(
                frac * local, frac * local * f_norm / (f_norm.mean() + 1e-12)
            )
            X = manifold.exp(X, step * (force / f_norm))
        return X

RepulsionNetExtrinsicSpeedup

Bases: RepulsionNet

:class:RepulsionNet with an extrinsic neighbour search (optimisation).

Overrides only the neighbour step, replacing the O(n^2) geodesic search by an O(n log n) Euclidean k-NN (KD-tree) on manifold.embed(X). The repulsion dynamics are otherwise identical.

Valid ONLY when embed is an isometric embedding whose ambient (chordal) neighbour ordering matches the geodesic one -- true for a compact manifold embedded in ambient space (e.g. S^2 in R^3), FALSE for a quotient space: there, points glued across an identified boundary are geodesic neighbours yet lie far apart in the embedding, so the KD-tree returns the wrong neighbours. Use the intrinsic :class:RepulsionNet in that case.

Source code in kmeanssa_ng/riemannian_manifold/epsilon_net.py

class RepulsionNetExtrinsicSpeedup(RepulsionNet):
    """:class:`RepulsionNet` with an *extrinsic* neighbour search (optimisation).

    Overrides only the neighbour step, replacing the O(n^2) geodesic search by an
    O(n log n) Euclidean k-NN (KD-tree) on ``manifold.embed(X)``. The repulsion
    dynamics are otherwise identical.

    Valid ONLY when ``embed`` is an isometric embedding whose ambient (chordal)
    neighbour ordering matches the geodesic one -- true for a compact manifold
    embedded in ambient space (e.g. S^2 in R^3), FALSE for a quotient space: there,
    points glued across an identified boundary are geodesic neighbours yet lie far
    apart in the embedding, so the KD-tree returns the wrong neighbours. Use the
    intrinsic :class:`RepulsionNet` in that case.
    """

    def _neighbours(self, manifold: RiemannianManifold, X: np.ndarray) -> np.ndarray:
        embedded = manifold.embed(X)
        _, idx = (
            NearestNeighbors(n_neighbors=self.k + 1).fit(embedded).kneighbors(embedded)
        )
        return idx[:, 1:]  # drop self

RiemannianCenter

Bases: RiemannianPoint, Center

A movable cluster center on a Riemannian manifold.

Centers can perform: - Brownian motion: Random walk using Geomstats BrownianMotion - Drift: Directed movement toward target points along geodesics

Attributes:

Name	Type	Description
space	RiemannianManifold	The Riemannian manifold this center belongs to.
coordinates	RiemannianManifold	The coordinates of the center on the manifold.

Example

from geomstats.geometry.hypersphere import Hypersphere
manifold = Hypersphere(dim=2)
space = RiemannianManifold(manifold)
point = RiemannianPoint(space, coordinates=np.array([1.0, 0.0, 0.0]))
center = RiemannianCenter(point)
center.brownian_motion(0.1)
center.drift(target_point, 0.5)

Source code in kmeanssa_ng/riemannian_manifold/center.py

class RiemannianCenter(RiemannianPoint, AbstractCenter):
    """A movable cluster center on a Riemannian manifold.

    Centers can perform:
    - Brownian motion: Random walk using Geomstats BrownianMotion
    - Drift: Directed movement toward target points along geodesics

    Attributes:
        space: The Riemannian manifold this center belongs to.
        coordinates: The coordinates of the center on the manifold.

    Example:
        ```python
        from geomstats.geometry.hypersphere import Hypersphere
        manifold = Hypersphere(dim=2)
        space = RiemannianManifold(manifold)
        point = RiemannianPoint(space, coordinates=np.array([1.0, 0.0, 0.0]))
        center = RiemannianCenter(point)
        center.brownian_motion(0.1)
        center.drift(target_point, 0.5)
        ```
    """

    def __init__(
        self,
        point: RiemannianPoint,
        rng: Generator | None = None,
    ) -> None:
        """Initialize a center from a point.

        Args:
            point: The initial point location.
            rng: Random number generator. If None, creates a new default_rng().
        """
        super().__init__(point.space, point.coordinates)
        self._rng = rng if rng is not None else default_rng()

    def brownian_motion(self, time_to_travel: float) -> None:
        """Perform Brownian motion on the Riemannian manifold.

        Implements a simple Brownian motion by:
        1. Generating a random tangent vector
        2. Scaling it by sqrt(time_to_travel)
        3. Moving along the geodesic in that direction

        Args:
            time_to_travel: Time parameter (distance ~ sqrt(time)).

        Raises:
            ValueError: If time_to_travel is negative or not numeric.
        """
        # Validate time_to_travel
        try:
            time_float = float(time_to_travel)
        except (TypeError, ValueError) as e:
            raise ValueError(
                f"time_to_travel must be a number, got {type(time_to_travel).__name__}"
            ) from e

        if time_float < 0:
            raise ValueError(f"time_to_travel must be non-negative, got {time_float}")

        if time_float == 0:
            return  # No movement

        # The Brownian increment is sqrt(time) * V, where V is a standard
        # Gaussian tangent vector (iid N(0,1) components in a metric-orthonormal
        # frame). random_tangent already returns that full Gaussian vector --
        # random in both direction *and* length -- so it is the whole
        # increment; multiplying it by an extra scalar N(0,1) would double-count
        # the noise (a product of two normals: kurtosis 9 instead of 3, wrong
        # radial law at finite step_size) and diverge from the 1-D graph step,
        # which draws a single Gaussian. It comes from the space (spaces without
        # a metric-orthonormal frame refuse) and from this center's own
        # generator, not geomstats' global RNG, so runs stay reproducible.
        tangent_vec = self.space.random_tangent(self.coordinates, self._rng)

        # Move along the exponential map (self.space dispatches to a closed form
        # on known manifolds, e.g. the sphere, and to geomstats otherwise).
        self.coordinates = self.space.exp(
            self.coordinates, np.sqrt(time_float) * tangent_vec
        )

    def drift(self, target_point: RiemannianPoint, prop_to_travel: float) -> None:
        """Move toward a target point along the geodesic.

        Moves a proportion of the geodesic distance to the target point.

        Args:
            target_point: The point to move toward.
            prop_to_travel: Proportion of distance to travel (0 to 1).

        Raises:
            ValueError: If target_point is None, prop_to_travel is not numeric,
                or prop_to_travel is not in [0, 1].

        Note:
            On manifolds where geodesics are not unique (e.g., antipodal points
            on a sphere), the drift may not move the center. This is a known
            limitation of Geomstats' geodesic computation. In practice, this
            rarely occurs due to the measure-zero probability of exact antipodal
            configurations, and the Brownian motion component of the simulated
            annealing algorithm provides thermal agitation to escape such
            degenerate configurations.
        """
        # Validate target_point
        if target_point is None:
            raise ValueError("target_point cannot be None")

        # Validate prop_to_travel
        try:
            prop_float = float(prop_to_travel)
        except (TypeError, ValueError) as e:
            raise ValueError(
                f"prop_to_travel must be a number, got {type(prop_to_travel).__name__}"
            ) from e

        if prop_float < 0 or prop_float > 1:
            raise ValueError(f"prop_to_travel must be in [0, 1], got {prop_float}")

        if prop_float == 0:
            return  # No movement

        # Move a proportion prop_float along the geodesic toward the target:
        # exp_p(prop * log_p(target)). This equals geomstats' geodesic evaluated at
        # prop, without building a geodesic callable per call, and dispatches to a
        # closed form on known manifolds (e.g. the sphere) via self.space.
        log_vec = self.space.log(self.coordinates, target_point.coordinates)
        new_coords = self.space.exp(self.coordinates, prop_float * log_vec)

        self.coordinates = new_coords

    def clone(self) -> RiemannianCenter:
        """Create an independent copy of this center.

        The cloned center shares the same manifold space but has
        independent coordinates. This is much faster than deepcopy.

        Returns:
            A new RiemannianCenter with the same location but independent state.

        Example:
            ```python
            original = RiemannianCenter(...)
            copy = original.clone()
            original.brownian_motion(0.1)  # Doesn't affect copy
            ```
        """
        # Create new center bypassing validation for speed
        new_center = object.__new__(RiemannianCenter)
        new_center._space = self._space
        new_center.coordinates = self.coordinates.copy()
        new_center._rng = self._rng
        return new_center

    def __repr__(self) -> str:
        """Detailed string representation."""
        return f"RiemannianCenter(coordinates={self.coordinates})"

    def __str__(self) -> str:
        """User-friendly string representation."""
        manifold_name = (
            self._space.manifold.__class__.__name__
            if hasattr(self._space, "manifold")
            else "Unknown"
        )
        return f"Center on {manifold_name} at {self.coordinates}"

__init__(point, rng=None)

Initialize a center from a point.

Parameters:

Name	Type	Description	Default
point	RiemannianPoint	The initial point location.	required
rng	Generator | None	Random number generator. If None, creates a new default_rng().	None

Source code in kmeanssa_ng/riemannian_manifold/center.py

def __init__(
    self,
    point: RiemannianPoint,
    rng: Generator | None = None,
) -> None:
    """Initialize a center from a point.

    Args:
        point: The initial point location.
        rng: Random number generator. If None, creates a new default_rng().
    """
    super().__init__(point.space, point.coordinates)
    self._rng = rng if rng is not None else default_rng()

__repr__()

Detailed string representation.

Source code in kmeanssa_ng/riemannian_manifold/center.py

def __repr__(self) -> str:
    """Detailed string representation."""
    return f"RiemannianCenter(coordinates={self.coordinates})"

__str__()

User-friendly string representation.

Source code in kmeanssa_ng/riemannian_manifold/center.py

def __str__(self) -> str:
    """User-friendly string representation."""
    manifold_name = (
        self._space.manifold.__class__.__name__
        if hasattr(self._space, "manifold")
        else "Unknown"
    )
    return f"Center on {manifold_name} at {self.coordinates}"

brownian_motion(time_to_travel)

Perform Brownian motion on the Riemannian manifold.

Implements a simple Brownian motion by: 1. Generating a random tangent vector 2. Scaling it by sqrt(time_to_travel) 3. Moving along the geodesic in that direction

Parameters:

Name	Type	Description	Default
time_to_travel	float	Time parameter (distance ~ sqrt(time)).	required

Raises:

Type	Description
ValueError	If time_to_travel is negative or not numeric.

Source code in kmeanssa_ng/riemannian_manifold/center.py

def brownian_motion(self, time_to_travel: float) -> None:
    """Perform Brownian motion on the Riemannian manifold.

    Implements a simple Brownian motion by:
    1. Generating a random tangent vector
    2. Scaling it by sqrt(time_to_travel)
    3. Moving along the geodesic in that direction

    Args:
        time_to_travel: Time parameter (distance ~ sqrt(time)).

    Raises:
        ValueError: If time_to_travel is negative or not numeric.
    """
    # Validate time_to_travel
    try:
        time_float = float(time_to_travel)
    except (TypeError, ValueError) as e:
        raise ValueError(
            f"time_to_travel must be a number, got {type(time_to_travel).__name__}"
        ) from e

    if time_float < 0:
        raise ValueError(f"time_to_travel must be non-negative, got {time_float}")

    if time_float == 0:
        return  # No movement

    # The Brownian increment is sqrt(time) * V, where V is a standard
    # Gaussian tangent vector (iid N(0,1) components in a metric-orthonormal
    # frame). random_tangent already returns that full Gaussian vector --
    # random in both direction *and* length -- so it is the whole
    # increment; multiplying it by an extra scalar N(0,1) would double-count
    # the noise (a product of two normals: kurtosis 9 instead of 3, wrong
    # radial law at finite step_size) and diverge from the 1-D graph step,
    # which draws a single Gaussian. It comes from the space (spaces without
    # a metric-orthonormal frame refuse) and from this center's own
    # generator, not geomstats' global RNG, so runs stay reproducible.
    tangent_vec = self.space.random_tangent(self.coordinates, self._rng)

    # Move along the exponential map (self.space dispatches to a closed form
    # on known manifolds, e.g. the sphere, and to geomstats otherwise).
    self.coordinates = self.space.exp(
        self.coordinates, np.sqrt(time_float) * tangent_vec
    )

clone()

Create an independent copy of this center.

The cloned center shares the same manifold space but has independent coordinates. This is much faster than deepcopy.

Returns:

Type	Description
RiemannianCenter	A new RiemannianCenter with the same location but independent state.

Example

original = RiemannianCenter(...)
copy = original.clone()
original.brownian_motion(0.1)  # Doesn't affect copy

Source code in kmeanssa_ng/riemannian_manifold/center.py

def clone(self) -> RiemannianCenter:
    """Create an independent copy of this center.

    The cloned center shares the same manifold space but has
    independent coordinates. This is much faster than deepcopy.

    Returns:
        A new RiemannianCenter with the same location but independent state.

    Example:
        ```python
        original = RiemannianCenter(...)
        copy = original.clone()
        original.brownian_motion(0.1)  # Doesn't affect copy
        ```
    """
    # Create new center bypassing validation for speed
    new_center = object.__new__(RiemannianCenter)
    new_center._space = self._space
    new_center.coordinates = self.coordinates.copy()
    new_center._rng = self._rng
    return new_center

drift(target_point, prop_to_travel)

Move toward a target point along the geodesic.

Moves a proportion of the geodesic distance to the target point.

Parameters:

Name	Type	Description	Default
target_point	RiemannianPoint	The point to move toward.	required
prop_to_travel	float	Proportion of distance to travel (0 to 1).	required

Raises:

Type	Description
ValueError	If target_point is None, prop_to_travel is not numeric, or prop_to_travel is not in [0, 1].

Note

On manifolds where geodesics are not unique (e.g., antipodal points on a sphere), the drift may not move the center. This is a known limitation of Geomstats' geodesic computation. In practice, this rarely occurs due to the measure-zero probability of exact antipodal configurations, and the Brownian motion component of the simulated annealing algorithm provides thermal agitation to escape such degenerate configurations.

Source code in kmeanssa_ng/riemannian_manifold/center.py

def drift(self, target_point: RiemannianPoint, prop_to_travel: float) -> None:
    """Move toward a target point along the geodesic.

    Moves a proportion of the geodesic distance to the target point.

    Args:
        target_point: The point to move toward.
        prop_to_travel: Proportion of distance to travel (0 to 1).

    Raises:
        ValueError: If target_point is None, prop_to_travel is not numeric,
            or prop_to_travel is not in [0, 1].

    Note:
        On manifolds where geodesics are not unique (e.g., antipodal points
        on a sphere), the drift may not move the center. This is a known
        limitation of Geomstats' geodesic computation. In practice, this
        rarely occurs due to the measure-zero probability of exact antipodal
        configurations, and the Brownian motion component of the simulated
        annealing algorithm provides thermal agitation to escape such
        degenerate configurations.
    """
    # Validate target_point
    if target_point is None:
        raise ValueError("target_point cannot be None")

    # Validate prop_to_travel
    try:
        prop_float = float(prop_to_travel)
    except (TypeError, ValueError) as e:
        raise ValueError(
            f"prop_to_travel must be a number, got {type(prop_to_travel).__name__}"
        ) from e

    if prop_float < 0 or prop_float > 1:
        raise ValueError(f"prop_to_travel must be in [0, 1], got {prop_float}")

    if prop_float == 0:
        return  # No movement

    # Move a proportion prop_float along the geodesic toward the target:
    # exp_p(prop * log_p(target)). This equals geomstats' geodesic evaluated at
    # prop, without building a geodesic callable per call, and dispatches to a
    # closed form on known manifolds (e.g. the sphere) via self.space.
    log_vec = self.space.log(self.coordinates, target_point.coordinates)
    new_coords = self.space.exp(self.coordinates, prop_float * log_vec)

    self.coordinates = new_coords

RiemannianManifold

Bases: Space

A Riemannian manifold space using geomstats.

This class wraps a geomstats manifold object and implements the Space interface for k-means clustering on Riemannian manifolds.

Attributes:

Name	Type	Description
manifold		The geomstats manifold object.

Note

On manifolds with non-unique geodesics (e.g., antipodal points on spheres), the drift operation may exhibit degenerate behavior where centers do not move toward their targets. This is a known limitation of geodesic computation. The Brownian motion in the simulated annealing algorithm provides thermal agitation to escape such configurations.

Example

from geomstats.geometry.hypersphere import Hypersphere
sphere = Hypersphere(dim=2)
space = RiemannianManifold(sphere)
points = space.sample_points(100, strategy=UniformManifoldSampling())
centers = [space.center_from_point(p) for p in points[:3]]
energy = space.calculate_energy(centers, observations=points)

Source code in kmeanssa_ng/riemannian_manifold/space.py

class RiemannianManifold(Space):
    """A Riemannian manifold space using geomstats.

    This class wraps a geomstats manifold object and implements the Space
    interface for k-means clustering on Riemannian manifolds.

    Attributes:
        manifold: The geomstats manifold object.

    Note:
        On manifolds with non-unique geodesics (e.g., antipodal points on spheres),
        the drift operation may exhibit degenerate behavior where centers do not
        move toward their targets. This is a known limitation of geodesic
        computation. The Brownian motion in the simulated annealing algorithm
        provides thermal agitation to escape such configurations.

    Example:
        ```python
        from geomstats.geometry.hypersphere import Hypersphere
        sphere = Hypersphere(dim=2)
        space = RiemannianManifold(sphere)
        points = space.sample_points(100, strategy=UniformManifoldSampling())
        centers = [space.center_from_point(p) for p in points[:3]]
        energy = space.calculate_energy(centers, observations=points)
        ```
    """

    def __init__(self, manifold) -> None:
        """Initialize a Riemannian manifold space.

        Args:
            manifold: A geomstats manifold object (e.g., Hypersphere, Hyperboloid).
        """
        self.manifold = manifold

    def distance(self, point1: RiemannianPoint, point2: RiemannianPoint) -> float:
        """Compute the geodesic distance between two points.

        Uses the manifold's Riemannian metric to compute the distance.

        Args:
            point1: First point.
            point2: Second point.

        Returns:
            The geodesic distance between point1 and point2.
        """
        dist = self.manifold.metric.dist(point1.coordinates, point2.coordinates)
        # Handle both scalar and 0-d array results
        return float(np.asarray(dist).item())

    def center_from_point(self, point: RiemannianPoint) -> RiemannianCenter:
        """Create a RiemannianCenter object from a RiemannianPoint object."""
        return RiemannianCenter(point)

    # ------------------------------------------------------------------
    # Geodesic operations on ambient coordinate arrays.
    # These let geodesic strategies (e.g. epsilon-net construction) drive the
    # manifold without touching geomstats or the RiemannianPoint wrapper.
    # ------------------------------------------------------------------
    @property
    def dim(self) -> int:
        """Intrinsic dimension of the manifold."""
        return self.manifold.dim

    @property
    def is_sphere(self) -> bool:
        """Whether the underlying manifold is a hypersphere."""
        return isinstance(self.manifold, Hypersphere)

    def exp(self, base_point: np.ndarray, tangent_vec: np.ndarray) -> np.ndarray:
        """Riemannian exponential: retract ``tangent_vec`` from ``base_point``.

        Batched over leading axes; arrays are ambient (extrinsic) coordinates.
        """
        return np.asarray(self.manifold.metric.exp(tangent_vec, base_point))

    def log(self, base_point: np.ndarray, point: np.ndarray) -> np.ndarray:
        """Riemannian logarithm: tangent at ``base_point`` pointing to ``point``.

        Inverse of :meth:`exp`; its Riemannian norm is the geodesic distance.
        """
        return np.asarray(self.manifold.metric.log(point, base_point))

    def norm(self, base_point: np.ndarray, tangent_vec: np.ndarray) -> np.ndarray:
        """Riemannian norm of ``tangent_vec`` in the tangent space at ``base_point``."""
        return np.asarray(self.manifold.metric.norm(tangent_vec, base_point))

    def to_tangent(self, base_point: np.ndarray, ambient_vec: np.ndarray) -> np.ndarray:
        """Project an ambient vector onto the tangent space at ``base_point``."""
        return np.asarray(self.manifold.to_tangent(ambient_vec, base_point))

    def random_tangent(
        self, base_point: np.ndarray, rng: np.random.Generator
    ) -> np.ndarray:
        """Tangent vector with iid standard-normal components in an
        orthonormal basis of the metric at ``base_point``.

        This is the Brownian direction of the annealing: drawing it in a
        metric-orthonormal frame is what makes the exploration isotropic and
        its law identical at every point. On a hypersphere the ambient
        Euclidean frame is orthonormal for the induced metric, so projecting
        an ambient Gaussian is exact. For a general metric it is **not**
        (on a conformal chart the ambient draw is up to several times longer
        near the boundary than at the center), so spaces without a known
        orthonormal frame refuse rather than silently bias the dynamics.

        Raises:
            NotImplementedError: If no metric-isotropic draw is implemented
                for this manifold. Override this method with a draw in an
                orthonormal frame of the metric.
        """
        if self.is_sphere:
            ambient = rng.standard_normal(self.shape)
            return self.to_tangent(base_point, ambient)
        raise NotImplementedError(
            "random_tangent (the metric-isotropic Brownian direction) is not "
            f"implemented for {type(self.manifold).__name__}: an ambient "
            "Gaussian is not isotropic for a general Riemannian metric and "
            "would silently bias the annealing. Override random_tangent with "
            "a draw in an orthonormal frame of the metric."
        )

    @property
    def shape(self) -> tuple:
        """Ambient shape of a single point on the manifold."""
        return self.manifold.shape

    def embed(self, points: np.ndarray) -> np.ndarray:
        """Ambient coordinates of ``points`` (for neighbour search).

        Points are already extrinsic here, so this is the identity; it exists so
        geodesic strategies can stay agnostic to the coordinate representation.
        """
        return np.asarray(points)

    def random_uniform(
        self, n: int, random_state: int | np.random.Generator | None = None
    ) -> np.ndarray:
        """Sample ``n`` points uniformly, reproducibly from ``random_state``.

        Returns an ``(n, dim + 1)`` array of ambient coordinates.
        """
        rng = np.random.default_rng(random_state)
        if self.is_sphere:
            # Gaussian in the ambient space, normalised -> uniform on the sphere.
            x = rng.standard_normal((n, self.dim + 1))
            return x / np.linalg.norm(x, axis=1, keepdims=True)
        # TODO: area-uniform sampling for other manifolds. Non-compact domains
        # (e.g. hyperbolic) also need a bounded region + soft confinement; this
        # is handled at the epsilon-net strategy level.
        raise NotImplementedError(
            "random_uniform is currently implemented for hyperspheres only, not "
            f"{type(self.manifold).__name__}."
        )

    def distances_from_centers(
        self, centers: list[RiemannianCenter], target: RiemannianPoint
    ) -> np.ndarray:
        """Compute distances from multiple centers to a single target point.

        Args:
            centers: List of k centers to compute distances from.
            target: The target point.

        Returns:
            Array of shape (k,) with distances from each center to target.

        Example:
            ```python
            centers = space.sample_centers(5)
            target = space.sample_points(1)[0]
            distances = space.distances_from_centers(centers, target)
            closest_idx = np.argmin(distances)
            closest_center = centers[closest_idx]
            ```
        """
        distances = np.empty(len(centers))
        for i, center in enumerate(centers):
            distances[i] = self.distance(center, target)
        return distances

    def calculate_energy(
        self,
        centers: list[RiemannianCenter],
        how: str = "empirical",
        observations: list[RiemannianPoint] | None = None,
    ) -> float:
        """Calculate the k-means energy for the given centers.

        The energy is the mean squared distance from each observation to its
        nearest center. The observations must be passed explicitly: they
        belong to the algorithm evaluating the energy, not to the manifold,
        so several algorithms can share one space without interfering.

        Args:
            centers: List of cluster centers.
            how: Only ``"empirical"`` is supported on a manifold: there is no
                uniform distribution over the points of a continuous
                manifold, and no node measure to carry. Any other mode is an
                error rather than a silent reinterpretation.
            observations: The points defining the empirical measure
                (``RiemannianPoint`` instances or coordinate arrays);
                required.

        Returns:
            The k-means energy (mean squared distance to nearest center).

        Raises:
            ValueError: If the mode is not "empirical", observations are
                missing, or the centers list is empty.
        """
        if how == "obs":
            raise ValueError(
                "energy mode 'obs' is retired: on a manifold use "
                "how='empirical' with the explicit 'observations' list"
            )
        if how != "empirical":
            raise ValueError(
                f"a Riemannian manifold only supports how='empirical', got "
                f"{how!r}: there is no uniform or node-carried measure on a "
                "continuous manifold"
            )
        if not observations:
            raise ValueError(
                "energy mode 'empirical' requires the explicit "
                "'observations' list (there is no reference measure on the "
                "space itself)"
            )

        if len(centers) == 0:
            raise ValueError("Centers list cannot be empty")

        total_energy = 0.0

        # For each observation, find squared distance to nearest center
        for obs in observations:
            if isinstance(obs, RiemannianPoint):
                obs_point = obs
            else:  # It's a numpy array (coordinates)
                obs_point = RiemannianPoint(self, obs)

            min_dist_sq = min(
                self.distance(center, obs_point) ** 2 for center in centers
            )
            total_energy += min_dist_sq

        return total_energy / len(observations)

    def get_point_type(self) -> type[RiemannianPoint]:
        """Return the type of points in this space."""
        return RiemannianPoint

dim property

Intrinsic dimension of the manifold.

is_sphere property

Whether the underlying manifold is a hypersphere.

shape property

Ambient shape of a single point on the manifold.

__init__(manifold)

Initialize a Riemannian manifold space.

Parameters:

Name	Type	Description	Default
manifold		A geomstats manifold object (e.g., Hypersphere, Hyperboloid).	required

Source code in kmeanssa_ng/riemannian_manifold/space.py

def __init__(self, manifold) -> None:
    """Initialize a Riemannian manifold space.

    Args:
        manifold: A geomstats manifold object (e.g., Hypersphere, Hyperboloid).
    """
    self.manifold = manifold

calculate_energy(centers, how='empirical', observations=None)

Calculate the k-means energy for the given centers.

The energy is the mean squared distance from each observation to its nearest center. The observations must be passed explicitly: they belong to the algorithm evaluating the energy, not to the manifold, so several algorithms can share one space without interfering.

Parameters:

Name	Type	Description	Default
centers	list[RiemannianCenter]	List of cluster centers.	required
how	str	Only "empirical" is supported on a manifold: there is no uniform distribution over the points of a continuous manifold, and no node measure to carry. Any other mode is an error rather than a silent reinterpretation.	'empirical'
observations	list[RiemannianPoint] | None	The points defining the empirical measure (RiemannianPoint instances or coordinate arrays); required.	None

Returns:

Type	Description
float	The k-means energy (mean squared distance to nearest center).

Raises:

Type	Description
ValueError	If the mode is not "empirical", observations are missing, or the centers list is empty.

Source code in kmeanssa_ng/riemannian_manifold/space.py

def calculate_energy(
    self,
    centers: list[RiemannianCenter],
    how: str = "empirical",
    observations: list[RiemannianPoint] | None = None,
) -> float:
    """Calculate the k-means energy for the given centers.

    The energy is the mean squared distance from each observation to its
    nearest center. The observations must be passed explicitly: they
    belong to the algorithm evaluating the energy, not to the manifold,
    so several algorithms can share one space without interfering.

    Args:
        centers: List of cluster centers.
        how: Only ``"empirical"`` is supported on a manifold: there is no
            uniform distribution over the points of a continuous
            manifold, and no node measure to carry. Any other mode is an
            error rather than a silent reinterpretation.
        observations: The points defining the empirical measure
            (``RiemannianPoint`` instances or coordinate arrays);
            required.

    Returns:
        The k-means energy (mean squared distance to nearest center).

    Raises:
        ValueError: If the mode is not "empirical", observations are
            missing, or the centers list is empty.
    """
    if how == "obs":
        raise ValueError(
            "energy mode 'obs' is retired: on a manifold use "
            "how='empirical' with the explicit 'observations' list"
        )
    if how != "empirical":
        raise ValueError(
            f"a Riemannian manifold only supports how='empirical', got "
            f"{how!r}: there is no uniform or node-carried measure on a "
            "continuous manifold"
        )
    if not observations:
        raise ValueError(
            "energy mode 'empirical' requires the explicit "
            "'observations' list (there is no reference measure on the "
            "space itself)"
        )

    if len(centers) == 0:
        raise ValueError("Centers list cannot be empty")

    total_energy = 0.0

    # For each observation, find squared distance to nearest center
    for obs in observations:
        if isinstance(obs, RiemannianPoint):
            obs_point = obs
        else:  # It's a numpy array (coordinates)
            obs_point = RiemannianPoint(self, obs)

        min_dist_sq = min(
            self.distance(center, obs_point) ** 2 for center in centers
        )
        total_energy += min_dist_sq

    return total_energy / len(observations)

center_from_point(point)

Create a RiemannianCenter object from a RiemannianPoint object.

Source code in kmeanssa_ng/riemannian_manifold/space.py

def center_from_point(self, point: RiemannianPoint) -> RiemannianCenter:
    """Create a RiemannianCenter object from a RiemannianPoint object."""
    return RiemannianCenter(point)

distance(point1, point2)

Compute the geodesic distance between two points.

Uses the manifold's Riemannian metric to compute the distance.

Parameters:

Name	Type	Description	Default
point1	RiemannianPoint	First point.	required
point2	RiemannianPoint	Second point.	required

Returns:

Type	Description
float	The geodesic distance between point1 and point2.

Source code in kmeanssa_ng/riemannian_manifold/space.py

def distance(self, point1: RiemannianPoint, point2: RiemannianPoint) -> float:
    """Compute the geodesic distance between two points.

    Uses the manifold's Riemannian metric to compute the distance.

    Args:
        point1: First point.
        point2: Second point.

    Returns:
        The geodesic distance between point1 and point2.
    """
    dist = self.manifold.metric.dist(point1.coordinates, point2.coordinates)
    # Handle both scalar and 0-d array results
    return float(np.asarray(dist).item())

distances_from_centers(centers, target)

Compute distances from multiple centers to a single target point.

Parameters:

Name	Type	Description	Default
centers	list[RiemannianCenter]	List of k centers to compute distances from.	required
target	RiemannianPoint	The target point.	required

Returns:

Type	Description
ndarray	Array of shape (k,) with distances from each center to target.

Example

centers = space.sample_centers(5)
target = space.sample_points(1)[0]
distances = space.distances_from_centers(centers, target)
closest_idx = np.argmin(distances)
closest_center = centers[closest_idx]

Source code in kmeanssa_ng/riemannian_manifold/space.py

def distances_from_centers(
    self, centers: list[RiemannianCenter], target: RiemannianPoint
) -> np.ndarray:
    """Compute distances from multiple centers to a single target point.

    Args:
        centers: List of k centers to compute distances from.
        target: The target point.

    Returns:
        Array of shape (k,) with distances from each center to target.

    Example:
        ```python
        centers = space.sample_centers(5)
        target = space.sample_points(1)[0]
        distances = space.distances_from_centers(centers, target)
        closest_idx = np.argmin(distances)
        closest_center = centers[closest_idx]
        ```
    """
    distances = np.empty(len(centers))
    for i, center in enumerate(centers):
        distances[i] = self.distance(center, target)
    return distances

embed(points)

Ambient coordinates of points (for neighbour search).

Points are already extrinsic here, so this is the identity; it exists so geodesic strategies can stay agnostic to the coordinate representation.

Source code in kmeanssa_ng/riemannian_manifold/space.py

def embed(self, points: np.ndarray) -> np.ndarray:
    """Ambient coordinates of ``points`` (for neighbour search).

    Points are already extrinsic here, so this is the identity; it exists so
    geodesic strategies can stay agnostic to the coordinate representation.
    """
    return np.asarray(points)

exp(base_point, tangent_vec)

Riemannian exponential: retract tangent_vec from base_point.

Batched over leading axes; arrays are ambient (extrinsic) coordinates.

Source code in kmeanssa_ng/riemannian_manifold/space.py

def exp(self, base_point: np.ndarray, tangent_vec: np.ndarray) -> np.ndarray:
    """Riemannian exponential: retract ``tangent_vec`` from ``base_point``.

    Batched over leading axes; arrays are ambient (extrinsic) coordinates.
    """
    return np.asarray(self.manifold.metric.exp(tangent_vec, base_point))

get_point_type()

Return the type of points in this space.

Source code in kmeanssa_ng/riemannian_manifold/space.py

def get_point_type(self) -> type[RiemannianPoint]:
    """Return the type of points in this space."""
    return RiemannianPoint

log(base_point, point)

Riemannian logarithm: tangent at base_point pointing to point.

Inverse of :meth:exp; its Riemannian norm is the geodesic distance.

Source code in kmeanssa_ng/riemannian_manifold/space.py

def log(self, base_point: np.ndarray, point: np.ndarray) -> np.ndarray:
    """Riemannian logarithm: tangent at ``base_point`` pointing to ``point``.

    Inverse of :meth:`exp`; its Riemannian norm is the geodesic distance.
    """
    return np.asarray(self.manifold.metric.log(point, base_point))

norm(base_point, tangent_vec)

Riemannian norm of tangent_vec in the tangent space at base_point.

Source code in kmeanssa_ng/riemannian_manifold/space.py

def norm(self, base_point: np.ndarray, tangent_vec: np.ndarray) -> np.ndarray:
    """Riemannian norm of ``tangent_vec`` in the tangent space at ``base_point``."""
    return np.asarray(self.manifold.metric.norm(tangent_vec, base_point))

random_tangent(base_point, rng)

Tangent vector with iid standard-normal components in an orthonormal basis of the metric at base_point.

This is the Brownian direction of the annealing: drawing it in a metric-orthonormal frame is what makes the exploration isotropic and its law identical at every point. On a hypersphere the ambient Euclidean frame is orthonormal for the induced metric, so projecting an ambient Gaussian is exact. For a general metric it is not (on a conformal chart the ambient draw is up to several times longer near the boundary than at the center), so spaces without a known orthonormal frame refuse rather than silently bias the dynamics.

Raises:

Type	Description
NotImplementedError	If no metric-isotropic draw is implemented for this manifold. Override this method with a draw in an orthonormal frame of the metric.

Source code in kmeanssa_ng/riemannian_manifold/space.py

def random_tangent(
    self, base_point: np.ndarray, rng: np.random.Generator
) -> np.ndarray:
    """Tangent vector with iid standard-normal components in an
    orthonormal basis of the metric at ``base_point``.

    This is the Brownian direction of the annealing: drawing it in a
    metric-orthonormal frame is what makes the exploration isotropic and
    its law identical at every point. On a hypersphere the ambient
    Euclidean frame is orthonormal for the induced metric, so projecting
    an ambient Gaussian is exact. For a general metric it is **not**
    (on a conformal chart the ambient draw is up to several times longer
    near the boundary than at the center), so spaces without a known
    orthonormal frame refuse rather than silently bias the dynamics.

    Raises:
        NotImplementedError: If no metric-isotropic draw is implemented
            for this manifold. Override this method with a draw in an
            orthonormal frame of the metric.
    """
    if self.is_sphere:
        ambient = rng.standard_normal(self.shape)
        return self.to_tangent(base_point, ambient)
    raise NotImplementedError(
        "random_tangent (the metric-isotropic Brownian direction) is not "
        f"implemented for {type(self.manifold).__name__}: an ambient "
        "Gaussian is not isotropic for a general Riemannian metric and "
        "would silently bias the annealing. Override random_tangent with "
        "a draw in an orthonormal frame of the metric."
    )

random_uniform(n, random_state=None)

Sample n points uniformly, reproducibly from random_state.

Returns an (n, dim + 1) array of ambient coordinates.

Source code in kmeanssa_ng/riemannian_manifold/space.py

def random_uniform(
    self, n: int, random_state: int | np.random.Generator | None = None
) -> np.ndarray:
    """Sample ``n`` points uniformly, reproducibly from ``random_state``.

    Returns an ``(n, dim + 1)`` array of ambient coordinates.
    """
    rng = np.random.default_rng(random_state)
    if self.is_sphere:
        # Gaussian in the ambient space, normalised -> uniform on the sphere.
        x = rng.standard_normal((n, self.dim + 1))
        return x / np.linalg.norm(x, axis=1, keepdims=True)
    # TODO: area-uniform sampling for other manifolds. Non-compact domains
    # (e.g. hyperbolic) also need a bounded region + soft confinement; this
    # is handled at the epsilon-net strategy level.
    raise NotImplementedError(
        "random_uniform is currently implemented for hyperspheres only, not "
        f"{type(self.manifold).__name__}."
    )

to_tangent(base_point, ambient_vec)

Project an ambient vector onto the tangent space at base_point.

Source code in kmeanssa_ng/riemannian_manifold/space.py

def to_tangent(self, base_point: np.ndarray, ambient_vec: np.ndarray) -> np.ndarray:
    """Project an ambient vector onto the tangent space at ``base_point``."""
    return np.asarray(self.manifold.to_tangent(ambient_vec, base_point))

RiemannianPoint

Bases: Point

A point on a Riemannian manifold.

A Riemannian point is represented by its coordinates on the manifold. The coordinates are validated to ensure they belong to the manifold.

Attributes:

Name	Type	Description
space	RiemannianManifold	The Riemannian manifold space this point belongs to.
coordinates	RiemannianManifold	The coordinates of the point on the manifold.

Example

from geomstats.geometry.hypersphere import Hypersphere
manifold = Hypersphere(dim=2)
space = RiemannianManifold(manifold)
point = RiemannianPoint(space, coordinates=np.array([1.0, 0.0, 0.0]))

Source code in kmeanssa_ng/riemannian_manifold/point.py

class RiemannianPoint(AbstractPoint):
    """A point on a Riemannian manifold.

    A Riemannian point is represented by its coordinates on the manifold.
    The coordinates are validated to ensure they belong to the manifold.

    Attributes:
        space: The Riemannian manifold space this point belongs to.
        coordinates: The coordinates of the point on the manifold.

    Example:
        ```python
        from geomstats.geometry.hypersphere import Hypersphere
        manifold = Hypersphere(dim=2)
        space = RiemannianManifold(manifold)
        point = RiemannianPoint(space, coordinates=np.array([1.0, 0.0, 0.0]))
        ```
    """

    def __init__(
        self,
        space: RiemannianManifold,
        coordinates: np.ndarray,
    ) -> None:
        """Initialize a point on a Riemannian manifold.

        Args:
            space: The Riemannian manifold space containing this point.
            coordinates: The coordinates of the point on the manifold.

        Raises:
            ValueError: If space is None, coordinates are not a numpy array,
                or coordinates don't belong to the manifold.
        """
        if space is None:
            raise ValueError("space cannot be None")

        if not isinstance(coordinates, np.ndarray):
            raise ValueError(
                f"coordinates must be a numpy array, got {type(coordinates).__name__}"
            )

        self._space = space
        self._validate_and_set_coordinates(coordinates)

    def _validate_and_set_coordinates(self, coordinates: np.ndarray) -> None:
        """Validate and set the coordinates on the manifold.

        Args:
            coordinates: Coordinates on the manifold.

        Raises:
            ValueError: If the shape does not match the manifold, the
                coordinates are not finite, or they do not belong to the
                manifold.
        """
        manifold_shape = self._space.manifold.shape
        coords_shape = coordinates.shape if coordinates.ndim > 0 else ()

        if coords_shape != manifold_shape:
            raise ValueError(
                f"Coordinates shape {coords_shape} does not match the manifold "
                f"shape {manifold_shape}"
            )

        if not np.all(np.isfinite(coordinates)):
            raise ValueError(f"Coordinates must be finite, got {coordinates}")

        # Membership check: an off-manifold point would not crash, it would
        # silently produce plausible-looking geodesic distances.
        belongs = getattr(self._space.manifold, "belongs", None)
        if belongs is not None:
            try:
                on_manifold = bool(np.all(belongs(coordinates)))
            except NotImplementedError:
                on_manifold = True  # backend cannot decide; accept
            if not on_manifold:
                raise ValueError(
                    f"Coordinates {coordinates} do not belong to the manifold "
                    f"{type(self._space.manifold).__name__}"
                )

        self.coordinates = coordinates

    @property
    def space(self) -> RiemannianManifold:
        """The Riemannian manifold space this point belongs to."""
        return self._space

    def __str__(self) -> str:
        """String representation of the point."""
        manifold_name = (
            self._space.manifold.__class__.__name__
            if hasattr(self._space, "manifold")
            else "Unknown"
        )
        return f"RiemannianPoint on {manifold_name} at {self.coordinates}"

    def __repr__(self) -> str:
        """Detailed string representation."""
        return f"RiemannianPoint(coordinates={self.coordinates})"

space property

The Riemannian manifold space this point belongs to.

__init__(space, coordinates)

Initialize a point on a Riemannian manifold.

Parameters:

Name	Type	Description	Default
space	RiemannianManifold	The Riemannian manifold space containing this point.	required
coordinates	ndarray	The coordinates of the point on the manifold.	required

Raises:

Type	Description
ValueError	If space is None, coordinates are not a numpy array, or coordinates don't belong to the manifold.

Source code in kmeanssa_ng/riemannian_manifold/point.py

def __init__(
    self,
    space: RiemannianManifold,
    coordinates: np.ndarray,
) -> None:
    """Initialize a point on a Riemannian manifold.

    Args:
        space: The Riemannian manifold space containing this point.
        coordinates: The coordinates of the point on the manifold.

    Raises:
        ValueError: If space is None, coordinates are not a numpy array,
            or coordinates don't belong to the manifold.
    """
    if space is None:
        raise ValueError("space cannot be None")

    if not isinstance(coordinates, np.ndarray):
        raise ValueError(
            f"coordinates must be a numpy array, got {type(coordinates).__name__}"
        )

    self._space = space
    self._validate_and_set_coordinates(coordinates)

__repr__()

Detailed string representation.

Source code in kmeanssa_ng/riemannian_manifold/point.py

def __repr__(self) -> str:
    """Detailed string representation."""
    return f"RiemannianPoint(coordinates={self.coordinates})"

__str__()

String representation of the point.

Source code in kmeanssa_ng/riemannian_manifold/point.py

def __str__(self) -> str:
    """String representation of the point."""
    manifold_name = (
        self._space.manifold.__class__.__name__
        if hasattr(self._space, "manifold")
        else "Unknown"
    )
    return f"RiemannianPoint on {manifold_name} at {self.coordinates}"

Sphere

Bases: RiemannianManifold

Hypersphere with closed-form geodesic operations.

The unit sphere's exponential, logarithm, distance and tangent projection have simple closed forms. Overriding them here avoids geomstats' generic per-call overhead (~7x on a single point), which dominates the manifold annealing loop; the results match geomstats to machine precision. Every override is vectorised over leading axes, matching the base class.

Source code in kmeanssa_ng/riemannian_manifold/space.py

class Sphere(RiemannianManifold):
    """Hypersphere with closed-form geodesic operations.

    The unit sphere's exponential, logarithm, distance and tangent projection
    have simple closed forms. Overriding them here avoids geomstats' generic
    per-call overhead (~7x on a single point), which dominates the manifold
    annealing loop; the results match geomstats to machine precision. Every
    override is vectorised over leading axes, matching the base class.
    """

    def distance(self, point1: RiemannianPoint, point2: RiemannianPoint) -> float:
        inner = np.clip(np.dot(point1.coordinates, point2.coordinates), -1.0, 1.0)
        return float(np.arccos(inner))

    def exp(self, base_point: np.ndarray, tangent_vec: np.ndarray) -> np.ndarray:
        base_point = np.asarray(base_point, dtype=float)
        tangent_vec = np.asarray(tangent_vec, dtype=float)
        # Project onto the tangent space first (as geomstats' exp does), so any
        # radial component does not push the result off the sphere.
        radial = np.sum(tangent_vec * base_point, axis=-1, keepdims=True)
        tangent_vec = tangent_vec - radial * base_point
        norm = np.linalg.norm(tangent_vec, axis=-1, keepdims=True)
        # cos|v| * p + sin|v| * v/|v|, with the |v| -> 0 limit equal to p.
        direction = np.divide(
            tangent_vec, norm, out=np.zeros_like(tangent_vec), where=norm > 1e-12
        )
        return np.cos(norm) * base_point + np.sin(norm) * direction

    def log(self, base_point: np.ndarray, point: np.ndarray) -> np.ndarray:
        base_point = np.asarray(base_point, dtype=float)
        point = np.asarray(point, dtype=float)
        inner = np.clip(np.sum(base_point * point, axis=-1, keepdims=True), -1.0, 1.0)
        proj = point - inner * base_point  # tangent component of `point` at base
        proj_norm = np.linalg.norm(proj, axis=-1, keepdims=True)
        direction = np.divide(
            proj, proj_norm, out=np.zeros_like(proj), where=proj_norm > 1e-12
        )
        return np.arccos(inner) * direction  # scaled by the geodesic distance

    def norm(self, base_point: np.ndarray, tangent_vec: np.ndarray) -> np.ndarray:
        # A tangent vector's Riemannian norm on the unit sphere is its ambient norm.
        return np.linalg.norm(np.asarray(tangent_vec, dtype=float), axis=-1)

    def to_tangent(self, base_point: np.ndarray, ambient_vec: np.ndarray) -> np.ndarray:
        base_point = np.asarray(base_point, dtype=float)
        ambient_vec = np.asarray(ambient_vec, dtype=float)
        inner = np.sum(ambient_vec * base_point, axis=-1, keepdims=True)
        return ambient_vec - inner * base_point

UniformNet

Bases: EpsilonNetStrategy

Plain uniform sampling. A baseline: fast, but irregular for finite n.

Source code in kmeanssa_ng/riemannian_manifold/epsilon_net.py

class UniformNet(EpsilonNetStrategy):
    """Plain uniform sampling. A baseline: fast, but irregular for finite n."""

    def build(self, manifold: RiemannianManifold, n: int) -> np.ndarray:
        return manifold.random_uniform(n, self.random_state)

approximate_geodesic_space(manifold, n, *, net=None, ell=None, random_state=None, intrinsic=False)

Approximate manifold by a quantum graph on an n-point epsilon-net.

Parameters:

Name	Type	Description	Default
manifold	RiemannianManifold	The geodesic space to approximate.	required
n	int	Number of net points.	required
net	EpsilonNetStrategy | None	Placement strategy (defaults to :class:RepulsionNet).	None
ell	float | None	Connection radius (defaults to sqrt of the covering radius).	None
random_state	int | Generator | None	Seed for the net and the covering-radius estimate.	None
intrinsic	bool	Build the net and graph without an ambient embedding (required for quotient spaces such as the Bolza surface).	False

Source code in kmeanssa_ng/riemannian_manifold/graph.py

def approximate_geodesic_space(
    manifold: RiemannianManifold,
    n: int,
    *,
    net: EpsilonNetStrategy | None = None,
    ell: float | None = None,
    random_state: int | np.random.Generator | None = None,
    intrinsic: bool = False,
) -> QuantumGraph:
    """Approximate ``manifold`` by a quantum graph on an ``n``-point epsilon-net.

    Args:
        manifold: The geodesic space to approximate.
        n: Number of net points.
        net: Placement strategy (defaults to :class:`RepulsionNet`).
        ell: Connection radius (defaults to sqrt of the covering radius).
        random_state: Seed for the net and the covering-radius estimate.
        intrinsic: Build the net and graph without an ambient embedding (required
            for quotient spaces such as the Bolza surface).
    """
    if net is None:
        net = RepulsionNet(random_state=random_state)
    points = net.build(manifold, n)
    return build_epsilon_net_graph(
        manifold, points, ell, random_state=random_state, intrinsic=intrinsic
    )

build_epsilon_net_graph(manifold, points, ell=None, *, precompute=True, covering_test=10000, random_state=None, intrinsic=False)

Connect an epsilon-net into a QuantumGraph (edges within ell).

Parameters:

Name	Type	Description	Default
manifold	RiemannianManifold	The manifold the points live on.	required
points	ndarray	An (n, d) epsilon-net (see :mod:epsilon_net).	required
ell	float | None	Connection radius l(epsilon). Defaults to sqrt(covering radius).	None
precompute	bool	Precompute pairwise shortest-path distances on the graph.	True
covering_test	int	Sample size used to estimate the covering radius.	10000
random_state	int | Generator | None	Seed for the covering-radius estimate.	None
intrinsic	bool	Build edges from exhaustive geodesic distances instead of an ambient nearest-neighbour prefilter. Required for spaces without a faithful embedding (e.g. quotients like the Bolza surface); O(n^2).	False

Returns:

Type	Description
QuantumGraph	A QuantumGraph with unit node/edge weights and geodesic edge lengths.

Raises:

Type	Description
ValueError	If the resulting graph is disconnected (net too sparse for ell; use more points or a larger ell).

Source code in kmeanssa_ng/riemannian_manifold/graph.py

def build_epsilon_net_graph(
    manifold: RiemannianManifold,
    points: np.ndarray,
    ell: float | None = None,
    *,
    precompute: bool = True,
    covering_test: int = 10000,
    random_state: int | np.random.Generator | None = None,
    intrinsic: bool = False,
) -> QuantumGraph:
    """Connect an epsilon-net into a QuantumGraph (edges within ``ell``).

    Args:
        manifold: The manifold the points live on.
        points: An ``(n, d)`` epsilon-net (see :mod:`epsilon_net`).
        ell: Connection radius l(epsilon). Defaults to sqrt(covering radius).
        precompute: Precompute pairwise shortest-path distances on the graph.
        covering_test: Sample size used to estimate the covering radius.
        random_state: Seed for the covering-radius estimate.
        intrinsic: Build edges from exhaustive geodesic distances instead of an
            ambient nearest-neighbour prefilter. Required for spaces without a
            faithful embedding (e.g. quotients like the Bolza surface); O(n^2).

    Returns:
        A QuantumGraph with unit node/edge weights and geodesic edge lengths.

    Raises:
        ValueError: If the resulting graph is disconnected (net too sparse for
            ``ell``; use more points or a larger ``ell``).
    """
    n = len(points)
    if ell is None:
        eps = estimate_covering_radius(
            manifold, points, covering_test, random_state, intrinsic=intrinsic
        )
        ell = float(np.sqrt(eps))

    if intrinsic:
        # No embedding: filter on the exact all-pairs geodesic (quotient-aware).
        dmat = _geodesic_matrix(manifold, points, points)
        ii, jj = np.triu_indices(n, k=1)
        lengths = dmat[ii, jj]
        keep = lengths <= ell
        ii, jj, lengths = ii[keep], jj[keep], lengths[keep]
    else:
        # Ambient radius ell captures every pair within geodesic distance ell
        # (ambient distance <= geodesic distance); the exact geodesic filter follows.
        embedded = manifold.embed(points)
        adjacency = (
            NearestNeighbors(radius=ell)
            .fit(embedded)
            .radius_neighbors_graph(embedded, radius=ell, mode="connectivity")
            .tocoo()
        )
        upper = adjacency.row < adjacency.col
        ii, jj = adjacency.row[upper], adjacency.col[upper]
        lengths = _geodesic(manifold, points[ii], points[jj])
        keep = lengths <= ell
        ii, jj, lengths = ii[keep], jj[keep], lengths[keep]

    graph = nx.Graph()
    graph.add_nodes_from(range(n))
    graph.add_weighted_edges_from(
        zip(ii.tolist(), jj.tolist(), lengths.tolist()), weight="length"
    )
    if not nx.is_connected(graph):
        raise ValueError(
            "The epsilon-net graph is disconnected; use more points or a larger ell."
        )

    qg = QuantumGraph(graph, precompute=False)
    nx.set_node_attributes(qg, 1.0, "weight")
    nx.set_edge_attributes(qg, 1.0, "weight")
    if precompute:
        qg.precomputing()
    return qg

create_bolza_surface()

Create the Bolza surface, a compact genus-2 hyperbolic space.

The Bolza surface is the quotient of the Poincaré disk by the Fuchsian group gluing opposite sides of a regular hyperbolic octagon (interior angle pi/4). It has constant curvature -1 and is the most symmetric closed genus-2 surface -- a negatively curved counterpart to the sphere. Unlike the other factories it carries no geomstats backend: its geodesic operations are closed-form and quotient-aware, so it plugs into the intrinsic epsilon-net strategies.

Returns:

Name	Type	Description
A	BolzaSurface	class:BolzaSurface. Points are (2,) real (Re, Im) coordinates
	BolzaSurface	in the fundamental octagon.

Example

surface = create_bolza_surface()
net = surface.random_uniform(500, random_state=0)

Source code in kmeanssa_ng/riemannian_manifold/generators.py

def create_bolza_surface() -> BolzaSurface:
    """Create the Bolza surface, a compact genus-2 hyperbolic space.

    The Bolza surface is the quotient of the Poincaré disk by the Fuchsian group
    gluing opposite sides of a regular hyperbolic octagon (interior angle pi/4).
    It has constant curvature -1 and is the most symmetric closed genus-2 surface
    -- a negatively curved counterpart to the sphere. Unlike the other factories
    it carries no geomstats backend: its geodesic operations are closed-form and
    quotient-aware, so it plugs into the intrinsic epsilon-net strategies.

    Returns:
        A :class:`BolzaSurface`. Points are ``(2,)`` real ``(Re, Im)`` coordinates
        in the fundamental octagon.

    Example:
        ```python
        surface = create_bolza_surface()
        net = surface.random_uniform(500, random_state=0)
        ```
    """
    return BolzaSurface()

create_hyperbolic_space(dim=2, **kwargs)

Create a hyperbolic space wrapped in a RiemannianManifold.

Creates the hyperboloid model of hyperbolic space H^dim.

Parameters:

Name	Type	Description	Default
dim	int	Dimension of the hyperbolic space (default: 2).	2
**kwargs		Additional arguments passed to Hyperboloid constructor.	{}

Returns:

Type	Description
RiemannianManifold	A RiemannianManifold wrapping the hyperbolic space.

Raises:

Type	Description
ValueError	If dim is not a positive integer (raised by Geomstats).
TypeError	If dim is not a valid type (raised by Geomstats).

Example

# Create 2D hyperbolic space (Poincaré disk)
hyperbolic = create_hyperbolic_space(dim=2)
points = hyperbolic.sample_points(100)

Note

This uses the hyperboloid model of hyperbolic geometry.

Source code in kmeanssa_ng/riemannian_manifold/generators.py

def create_hyperbolic_space(dim: int = 2, **kwargs) -> RiemannianManifold:
    """Create a hyperbolic space wrapped in a RiemannianManifold.

    Creates the hyperboloid model of hyperbolic space H^dim.

    Args:
        dim: Dimension of the hyperbolic space (default: 2).
        **kwargs: Additional arguments passed to Hyperboloid constructor.

    Returns:
        A RiemannianManifold wrapping the hyperbolic space.

    Raises:
        ValueError: If dim is not a positive integer (raised by Geomstats).
        TypeError: If dim is not a valid type (raised by Geomstats).

    Example:
        ```python
        # Create 2D hyperbolic space (Poincaré disk)
        hyperbolic = create_hyperbolic_space(dim=2)
        points = hyperbolic.sample_points(100)
        ```

    Note:
        This uses the hyperboloid model of hyperbolic geometry.
    """
    from geomstats.geometry.hyperboloid import Hyperboloid

    # Geomstats handles validation
    hyperbolic_manifold = Hyperboloid(dim=dim, **kwargs)
    return RiemannianManifold(hyperbolic_manifold)

create_sphere(dim=2, **kwargs)

Create a hypersphere space wrapped in a RiemannianManifold.

Creates a hypersphere S^dim embedded in R^(dim+1). The sphere is equipped with the standard round metric inherited from the Euclidean ambient space.

Parameters:

Name	Type	Description	Default
dim	int	Dimension of the sphere (default: 2 for the standard 2-sphere S^2). - dim=1: Circle S^1 in R^2 - dim=2: Standard sphere S^2 in R^3 - dim=3: 3-sphere S^3 in R^4 etc.	2
**kwargs		Additional arguments passed to Hypersphere constructor.	{}

Returns:

Type	Description
RiemannianManifold	A RiemannianManifold wrapping the hypersphere.

Raises:

Type	Description
ValueError	If dim is not a positive integer (raised by Geomstats).
TypeError	If dim is not a valid type (raised by Geomstats).

Example

# Create a 2-sphere (surface of a ball in 3D)
sphere = create_sphere(dim=2)
points = sphere.sample_points(100)

# Create a circle
circle = create_sphere(dim=1)

Note

The sphere S^dim is the set of unit vectors in R^(dim+1): S^dim = {x ∈ R^(dim+1) : ||x|| = 1}

Source code in kmeanssa_ng/riemannian_manifold/generators.py

def create_sphere(dim: int = 2, **kwargs) -> RiemannianManifold:
    """Create a hypersphere space wrapped in a RiemannianManifold.

    Creates a hypersphere S^dim embedded in R^(dim+1). The sphere is equipped
    with the standard round metric inherited from the Euclidean ambient space.

    Args:
        dim: Dimension of the sphere (default: 2 for the standard 2-sphere S^2).
            - dim=1: Circle S^1 in R^2
            - dim=2: Standard sphere S^2 in R^3
            - dim=3: 3-sphere S^3 in R^4
            etc.
        **kwargs: Additional arguments passed to Hypersphere constructor.

    Returns:
        A RiemannianManifold wrapping the hypersphere.

    Raises:
        ValueError: If dim is not a positive integer (raised by Geomstats).
        TypeError: If dim is not a valid type (raised by Geomstats).

    Example:
        ```python
        # Create a 2-sphere (surface of a ball in 3D)
        sphere = create_sphere(dim=2)
        points = sphere.sample_points(100)

        # Create a circle
        circle = create_sphere(dim=1)
        ```

    Note:
        The sphere S^dim is the set of unit vectors in R^(dim+1):
        S^dim = {x ∈ R^(dim+1) : ||x|| = 1}
    """
    # Geomstats handles validation
    #  Note: intrinsic coordinates cause issues with belongs() and BrownianMotion
    # For now, use extrinsic (default) coordinates
    sphere_manifold = Hypersphere(dim=dim, **kwargs)
    return Sphere(sphere_manifold)

estimate_covering_radius(manifold, points, n_test=10000, random_state=None, *, intrinsic=False)

Estimate the covering radius max_z min_i d(z, x_i) by dense sampling.

Nearest net points are found in ambient coordinates (exact when the ambient distance is monotone in the geodesic one, e.g. on the sphere), then the geodesic distance to that nearest point is measured exactly. With intrinsic=True the nearest net point is instead found by exhaustive geodesic distance -- for spaces without a faithful embedding (e.g. quotients), at O(n_test * n) cost.

Source code in kmeanssa_ng/riemannian_manifold/graph.py

def estimate_covering_radius(
    manifold: RiemannianManifold,
    points: np.ndarray,
    n_test: int = 10000,
    random_state: int | np.random.Generator | None = None,
    *,
    intrinsic: bool = False,
) -> float:
    """Estimate the covering radius max_z min_i d(z, x_i) by dense sampling.

    Nearest net points are found in ambient coordinates (exact when the ambient
    distance is monotone in the geodesic one, e.g. on the sphere), then the
    geodesic distance to that nearest point is measured exactly. With
    ``intrinsic=True`` the nearest net point is instead found by exhaustive
    geodesic distance -- for spaces without a faithful embedding (e.g. quotients),
    at O(n_test * n) cost.
    """
    rng = np.random.default_rng(random_state)
    test = manifold.random_uniform(n_test, rng)
    if intrinsic:
        return float(_geodesic_matrix(manifold, test, points).min(axis=1).max())
    _, idx = (
        NearestNeighbors(n_neighbors=1)
        .fit(manifold.embed(points))
        .kneighbors(manifold.embed(test))
    )
    nearest = points[idx[:, 0]]
    return float(_geodesic(manifold, test, nearest).max())

:::