diff --git a/CHANGELOG.md b/CHANGELOG.md index 84ced0a..58514ee 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,40 @@ and the project adheres to [Semantic Versioning](https://semver.org/). ## [Unreleased] +### Changed + +- **`LinOp` endpoints are bound to `VectorSpace`** rather than `CoordinateSpace`. + A linear operator needs only `zeros`/`add`/`scale` from its endpoints; + coordinates are a representation, not part of the contract. Operators can now + be defined over any vector space, including one with no `shape`, `size`, + `flatten` or `stacked`. Concrete operators are unaffected — `SparseLinOp` and + `DiagonalLinOp` declare `LinOp[CoordinateSpace, CoordinateSpace]` themselves. + + `VectorSpace` and not `InnerProductSpace`: the two are siblings, so binding to + the latter would have *narrowed* the accepted set and dropped bare `TreeSpace` + and `StackedSpace`. + +- **Operations that need more than a vector space now say so.** Materialization + (`to_dense`, `to_matrix`, `fuse`, basis probes) requires coordinates, batching + (`vapply`, `rvapply`) requires the batched coordinate surface, and `rapply` + requires an inner product for its metric-adjoint contract. Each raises + `CapabilityError` naming the operation and the missing capability, where + previously an unsupported endpoint surfaced as an `AttributeError`. + +- `MatrixFreeLinOp`'s strict adjoint-consistency probe no longer requires + coordinates. It builds its probe element from the space's own `ones` when one + exists, so a non-coordinate space is checked too, and is skipped only when no + non-zero element can be built at all. + +- Operator endpoint compatibility no longer compares `shape`. Type identity, + context agreement and convert-equality already decide it, and + `CoordinateSpace` equality compares shape. + +### Added + +- `CapabilityError`, exported from `spacecore` and `spacecore.space`, raised when + a space lacks a capability the requested operation needs. + ## [0.4.3] — 2026-09-12 ### Removed diff --git a/spacecore/__init__.py b/spacecore/__init__.py index e65d33e..a8d76e6 100644 --- a/spacecore/__init__.py +++ b/spacecore/__init__.py @@ -110,6 +110,7 @@ Space, StarSpace, SpaceCheck, + CapabilityError, SpaceValidationError, SquareMatrixCheck, VectorSpace, @@ -243,6 +244,7 @@ "EuclideanElementwiseJordanSpace", "Space", "SpaceCheck", + "CapabilityError", "SpaceValidationError", "SquareMatrixCheck", "DenseArray", diff --git a/spacecore/_batching.py b/spacecore/_batching.py index bb7b8f8..73ccc4e 100644 --- a/spacecore/_batching.py +++ b/spacecore/_batching.py @@ -4,6 +4,7 @@ from typing import Any from .space.checks import _run_checks +from .space._capabilities import _CAP_BATCH, _CAP_INNER, require # Shared batched-evaluation helpers used by both LinOp and Functional. They live # here (rather than in functional/_base) so the linop and functional batched @@ -59,6 +60,8 @@ def _warn_vmap_fallback_once(obj: Any, method: str, batch_size: int) -> None: def _batched_inner(space: Any, xs: Any, ys: Any) -> Any: """Return ``space.inner(xs[i], ys[i])`` for a leading-axis batch.""" + require(space, _CAP_BATCH, "_batched_inner") + require(space, _CAP_INNER, "_batched_inner") xs_flat = space.flatten_batch(xs) ys_dual = ys if space.is_euclidean else space.riesz(ys) ys_flat = space.flatten_batch(ys_dual) diff --git a/spacecore/_errors.py b/spacecore/_errors.py index 80ea110..51d51ab 100644 --- a/spacecore/_errors.py +++ b/spacecore/_errors.py @@ -27,3 +27,7 @@ class ContextConflictError(ContextError): class UnknownBackendError(ContextError): """No backend is registered under the requested family name.""" + + +class CapabilityError(TypeError): + """A space lacks a capability the requested operation needs.""" diff --git a/spacecore/linop/_algebra.py b/spacecore/linop/_algebra.py index ce8f33c..cd7e991 100644 --- a/spacecore/linop/_algebra.py +++ b/spacecore/linop/_algebra.py @@ -1,13 +1,14 @@ from __future__ import annotations import warnings -from math import prod from typing import Any, Callable, Sequence, cast from ._base import LinOp, Domain, Codomain from ._metric import _requires_euclidean_or_riesz, metric_rapply, metric_rvapply from .._check_policy import CheckLevel, minimum_check_level from .._checks import checked_method +from ..space.base import InnerProductSpace +from ..space._capabilities import _CAP_BATCH, _CAP_COORDINATE, probe_element, require from .._repr import summarize_value from ..contextual import Context from ..kernels import core_kernels @@ -58,8 +59,6 @@ def _same_space_for_algebra(left: Any, right: Any) -> bool: return True if type(left) is not type(right): return False - if tuple(left.shape) != tuple(right.shape): - return False if not left.ctx.same_math(right.ctx): return False try: @@ -269,10 +268,14 @@ def rapply(self, y: Any) -> Any: @checked_method(in_space="domain", out_space="codomain", in_batched=True, out_batched=True) def vapply(self, xs: Any) -> Any: """Return ``scalar * op.vapply(xs)``.""" + require(self.domain, _CAP_BATCH, "ScaledLinOp.vapply") + require(self.codomain, _CAP_BATCH, "ScaledLinOp.vapply") return self._vapply_core(xs) def rvapply(self, ys: Any) -> Any: """Return ``conj(scalar) * op.rvapply(ys)``.""" + require(self.domain, _CAP_BATCH, "ScaledLinOp.rvapply") + require(self.codomain, _CAP_BATCH, "ScaledLinOp.rvapply") xs = self.op.rvapply(ys) return self.domain.scale_batch(_conjugate_scalar(self.scalar), xs) @@ -419,11 +422,15 @@ def rapply(self, y: Any) -> Any: @checked_method(in_space="domain", out_space="codomain", in_batched=True, out_batched=True) def vapply(self, xs: Any) -> Any: """Return ``sum_i ops[i].vapply(xs)``.""" + require(self.domain, _CAP_BATCH, "SumLinOp.vapply") + require(self.codomain, _CAP_BATCH, "SumLinOp.vapply") return self._vapply_core(xs) @checked_method(in_space="codomain", out_space="domain", in_batched=True, out_batched=True) def rvapply(self, ys: Any) -> Any: """Return ``sum_i ops[i].rvapply(ys)``.""" + require(self.domain, _CAP_BATCH, "SumLinOp.rvapply") + require(self.codomain, _CAP_BATCH, "SumLinOp.rvapply") add_batch = self.domain.add_batch acc = self.ops_tuple[0].rvapply(ys) for op in self.ops_tuple[1:]: @@ -578,10 +585,14 @@ def rapply(self, z: Any) -> Any: @checked_method(in_space="domain", out_space="codomain", in_batched=True, out_batched=True) def vapply(self, xs: Any) -> Any: """Return ``left.vapply(right.vapply(xs))``.""" + require(self.domain, _CAP_BATCH, "ComposedLinOp.vapply") + require(self.codomain, _CAP_BATCH, "ComposedLinOp.vapply") return self._vapply_core(xs) def rvapply(self, zs: Any) -> Any: """Return ``right.rvapply(left.rvapply(zs))``.""" + require(self.domain, _CAP_BATCH, "ComposedLinOp.rvapply") + require(self.codomain, _CAP_BATCH, "ComposedLinOp.rvapply") return self.right.rvapply(self.left.rvapply(zs)) def fuse(self, *, materialize: bool = False) -> LinOp: @@ -705,11 +716,15 @@ def rapply(self, y: Any) -> Any: @checked_method(in_space="domain", in_batched=True) def vapply(self, xs: Any) -> Any: """Return the batched zero element of the codomain.""" + require(self.domain, _CAP_BATCH, "ZeroLinOp.vapply") + require(self.codomain, _CAP_BATCH, "ZeroLinOp.vapply") return self._vapply_core(xs) @checked_method(in_space="codomain", in_batched=True) def rvapply(self, ys: Any) -> Any: """Return the batched zero element of the domain.""" + require(self.domain, _CAP_BATCH, "ZeroLinOp.rvapply") + require(self.codomain, _CAP_BATCH, "ZeroLinOp.rvapply") return _batched_zeros(self.domain, _leading_shape(self.codomain, ys)) def to_dense(self) -> Any: @@ -718,6 +733,8 @@ def to_dense(self) -> Any: The returned array has shape ``self.codomain.shape + self.domain.shape``. """ + require(self.domain, _CAP_COORDINATE, "ZeroLinOp.to_dense") + require(self.codomain, _CAP_COORDINATE, "ZeroLinOp.to_dense") return self.ops.zeros( tuple(self.codomain.shape) + tuple(self.domain.shape), dtype=self.dtype ) @@ -802,11 +819,15 @@ def rapply(self, x: Any) -> Any: @checked_method(in_space="domain", in_batched=True) def vapply(self, xs: Any) -> Any: """Return ``xs`` after batched domain validation.""" + require(self.domain, _CAP_BATCH, "IdentityLinOp.vapply") + require(self.codomain, _CAP_BATCH, "IdentityLinOp.vapply") return xs @checked_method(in_space="codomain", in_batched=True) def rvapply(self, xs: Any) -> Any: """Return ``xs`` after batched codomain validation.""" + require(self.domain, _CAP_BATCH, "IdentityLinOp.rvapply") + require(self.codomain, _CAP_BATCH, "IdentityLinOp.rvapply") return xs def to_dense(self) -> Any: @@ -815,6 +836,8 @@ def to_dense(self) -> Any: The returned array has shape ``self.codomain.shape + self.domain.shape``. """ + require(self.domain, _CAP_COORDINATE, "IdentityLinOp.to_dense") + require(self.codomain, _CAP_COORDINATE, "IdentityLinOp.to_dense") size = 1 for dim in self.domain.shape: size *= dim @@ -999,7 +1022,12 @@ def __init__( raise TypeError(f"rvapply must be callable, got {type(rvapply).__name__}.") super().__init__(dom, cod, ctx, check_level=check_level) - non_euclidean = not (self.domain.is_euclidean and self.codomain.is_euclidean) + has_inner = all( + isinstance(space, InnerProductSpace) for space in (self.domain, self.codomain) + ) + non_euclidean = has_inner and not ( + self.domain.is_euclidean and self.codomain.is_euclidean + ) if euclidean_adjoint is True: # The caller supplied the *coordinate* adjoint A^dagger; wrap it into # the metric adjoint R_X^-1 A^dagger R_Y once, here. @@ -1062,20 +1090,21 @@ def _wrapped_rvapply(ys: Any) -> Any: def _check_adjoint_consistency(self) -> None: """Probe the declared adjoint identity on deterministic space elements.""" - if not all(hasattr(space, "inner") for space in (self.domain, self.codomain)): + if not all( + isinstance(space, InnerProductSpace) + for space in (self.domain, self.codomain) + ): return - def probe(space: Any) -> Any: - if hasattr(space, "ones"): - return space.ones() - if hasattr(space, "unflatten") and hasattr(space, "shape"): - flat = self.ops.ones((prod(space.shape),), dtype=self.dtype) - return space.unflatten(flat) - raise TypeError(f"{type(space).__name__} cannot build a strict probe element.") + x = probe_element(self.domain, self.ops, self.dtype) + y = probe_element(self.codomain, self.ops, self.dtype) + if x is None or y is None: + # Neither coordinates nor ``ones``: the space cannot build a non-zero + # element, so there is nothing to probe with. Construction-time + # geometry warnings still apply. + return try: - x = probe(self.domain) - y = probe(self.codomain) ax = self.apply_fn(x) ahy = self.rapply_fn(y) self.codomain._check_member(ax) @@ -1204,6 +1233,8 @@ def rapply(self, y: Any) -> Any: Any Element of ``self.domain`` returned by ``rapply_fn``. """ + require(self.domain, InnerProductSpace, "MatrixFreeLinOp.rapply") + require(self.codomain, InnerProductSpace, "MatrixFreeLinOp.rapply") return self._rapply_core(y) @checked_method(in_space="domain", out_space="codomain", in_batched=True, out_batched=True) @@ -1222,6 +1253,8 @@ def vapply(self, xs: Any) -> Any: Batched element of ``self.codomain`` produced by ``vapply_fn`` or by the fallback batching implementation. """ + require(self.domain, _CAP_BATCH, "MatrixFreeLinOp.vapply") + require(self.codomain, _CAP_BATCH, "MatrixFreeLinOp.vapply") if self.vapply_fn is None: return super().vapply(xs) return self.vapply_fn(xs) @@ -1242,6 +1275,8 @@ def rvapply(self, ys: Any) -> Any: Batched element of ``self.domain`` produced by ``rvapply_fn`` or by the fallback batching implementation. """ + require(self.domain, _CAP_BATCH, "MatrixFreeLinOp.rvapply") + require(self.codomain, _CAP_BATCH, "MatrixFreeLinOp.rvapply") if self.rvapply_fn is None: return super().rvapply(ys) return self.rvapply_fn(ys) @@ -1379,10 +1414,14 @@ def rapply(self, x: Any) -> Any: def vapply(self, ys: Any) -> Any: """Return ``op.rvapply(ys)`` over a batch.""" + require(self.domain, _CAP_BATCH, "_AdjointViewLinOp.vapply") + require(self.codomain, _CAP_BATCH, "_AdjointViewLinOp.vapply") return self.op.rvapply(ys) def rvapply(self, xs: Any) -> Any: """Return ``op.vapply(xs)`` over a batch.""" + require(self.domain, _CAP_BATCH, "_AdjointViewLinOp.rvapply") + require(self.codomain, _CAP_BATCH, "_AdjointViewLinOp.rvapply") return self.op.vapply(xs) def fuse(self, *, materialize: bool = False) -> LinOp: diff --git a/spacecore/linop/_base.py b/spacecore/linop/_base.py index 5fa30cc..fdc30e9 100644 --- a/spacecore/linop/_base.py +++ b/spacecore/linop/_base.py @@ -10,14 +10,15 @@ from .._batching import _leading_batch_size, _warn_vmap_fallback_once from .._check_policy import CheckLevel from .._checks import checked_method +from ..space._capabilities import _CAP_BATCH, _CAP_COORDINATE, require from ..backend import PyTreeNode from .._repr import describe_space -from ..space import CoordinateSpace +from ..space import VectorSpace from ..contextual import Context from ..contextual import ContextBound -Domain = TypeVar("Domain", bound=CoordinateSpace) -Codomain = TypeVar("Codomain", bound=CoordinateSpace) +Domain = TypeVar("Domain", bound=VectorSpace) +Codomain = TypeVar("Codomain", bound=VectorSpace) class LinOp(PyTreeNode, ContextBound, Generic[Domain, Codomain]): @@ -183,12 +184,16 @@ def is_hermitian(self) -> bool | None: @checked_method(in_space="domain", in_batched=True) def vapply(self, xs: Any) -> Any: """Apply over a leading batch axis. Input must have shape ``(N,) + domain.shape``; use ``moveaxis`` for other layouts.""" + require(self.domain, _CAP_BATCH, "LinOp.vapply") + require(self.codomain, _CAP_BATCH, "LinOp.vapply") _warn_vmap_fallback_once(self, "vapply", _leading_batch_size(self.domain, xs)) return self.ops.vmap(self.apply, in_axes=0, out_axes=0)(xs) @checked_method(in_space="codomain", in_batched=True) def rvapply(self, ys: Any) -> Any: """Apply the adjoint over a leading batch axis. Input must have shape ``(N,) + codomain.shape``; use ``moveaxis`` for other layouts.""" + require(self.domain, _CAP_BATCH, "LinOp.rvapply") + require(self.codomain, _CAP_BATCH, "LinOp.rvapply") _warn_vmap_fallback_once(self, "rvapply", _leading_batch_size(self.codomain, ys)) return self.ops.vmap(self.rapply, in_axes=0, out_axes=0)(ys) @@ -327,6 +332,8 @@ def to_dense(self) -> Any: that already store a dense or sparse matrix should override this method for efficiency. """ + require(self.domain, _CAP_COORDINATE, "LinOp.to_dense") + require(self.codomain, _CAP_COORDINATE, "LinOp.to_dense") return self.ops.reshape( self.to_matrix(), tuple(self.codomain.shape) + tuple(self.domain.shape) ) @@ -346,6 +353,8 @@ def to_matrix(self) -> Any: for small/testing use; concrete storage-backed subclasses should override it when they can expose a matrix directly. """ + require(self.domain, _CAP_COORDINATE, "LinOp.to_matrix") + require(self.codomain, _CAP_COORDINATE, "LinOp.to_matrix") domain_size = prod(self.domain.shape) codomain_size = prod(self.codomain.shape) eye = self.ops.eye(domain_size, dtype=self.dtype) diff --git a/spacecore/linop/_metric.py b/spacecore/linop/_metric.py index 2fc8450..8a70122 100644 --- a/spacecore/linop/_metric.py +++ b/spacecore/linop/_metric.py @@ -2,7 +2,8 @@ import warnings -from ..space.base import InnerProduct, InnerProductSpace +from ..space.base import CoordinateSpace, InnerProduct, InnerProductSpace +from ..space._capabilities import require _METRIC_BATCH_FALLBACK_ERRORS = (TypeError, ValueError, NotImplementedError) @@ -50,6 +51,8 @@ def _requires_euclidean_or_riesz(dom, cod, opname: str) -> None: def _metric_is_hermitian_by_basis(op) -> bool | None: """Check self-adjointness by comparing forward and adjoint basis actions.""" + if not isinstance(op.domain, CoordinateSpace): + return None if op.domain != op.codomain: return False if op.domain.size > _METRIC_HERMITIAN_BASIS_CHECK_MAX_SIZE: @@ -95,6 +98,8 @@ def metric_rapply(domain, codomain, euclidean_rapply, y): Getting this wrong is silent: a coordinate transpose satisfies the identity on every Euclidean space, so only a non-Euclidean test can detect it. """ + require(domain, InnerProductSpace, "metric_rapply") + require(codomain, InnerProductSpace, "metric_rapply") if domain.is_euclidean and codomain.is_euclidean: return euclidean_rapply(y) return domain.riesz_inverse(euclidean_rapply(codomain.riesz(y))) diff --git a/spacecore/space/__init__.py b/spacecore/space/__init__.py index ae0f1df..9f60be9 100644 --- a/spacecore/space/__init__.py +++ b/spacecore/space/__init__.py @@ -1,5 +1,6 @@ """Vector space abstractions, concrete spaces, and validation checks.""" +from .._errors import CapabilityError from .checks import ( BackendCheck, DTypeCheck, @@ -59,6 +60,7 @@ "TreeSpectralDecomposition", "Space", "SpaceCheck", + "CapabilityError", "SpaceValidationError", "SquareMatrixCheck", "VectorSpace", diff --git a/spacecore/space/_capabilities.py b/spacecore/space/_capabilities.py new file mode 100644 index 0000000..dc9f627 --- /dev/null +++ b/spacecore/space/_capabilities.py @@ -0,0 +1,92 @@ +"""Class-based space capability queries and operation requirements.""" +from __future__ import annotations + +from typing import Any, Protocol, runtime_checkable + +from .._errors import CapabilityError +from .base import ( + CoordinateSpace, + EuclideanJordanAlgebraSpace, + InnerProductSpace, + JordanAlgebraSpace, + Space, + StarSpace, +) + +CapabilitySet = frozenset[type] +_CAP_COORDINATE = CoordinateSpace +_CAP_BATCH = CoordinateSpace # Batching uses the coordinate surface (ADR-006). +_CAP_INNER = InnerProductSpace +_CAP_STAR = StarSpace +_CAP_JORDAN = JordanAlgebraSpace +_CAP_EUCLIDEAN_JORDAN = EuclideanJordanAlgebraSpace + + +def _space_capabilities(space: Space) -> CapabilitySet: + """Return structural capabilities advertised by one leaf space.""" + capabilities: set[type] = set() + if isinstance(space, CoordinateSpace): + capabilities.add(_CAP_COORDINATE) + if isinstance(space, InnerProductSpace): + capabilities.add(_CAP_INNER) + if isinstance(space, StarSpace): + capabilities.add(_CAP_STAR) + if isinstance(space, JordanAlgebraSpace): + capabilities.add(_CAP_JORDAN) + if isinstance(space, EuclideanJordanAlgebraSpace): + capabilities.add(_CAP_EUCLIDEAN_JORDAN) + return frozenset(capabilities) + + +@runtime_checkable +class SupportsOnes(Protocol): + """A space that can build a deterministic all-ones element. + + A declared capability type rather than an ad-hoc ``hasattr`` probe (ADR-005). + It is structural because ``ones`` is optional on spaces that are otherwise + unrelated -- ``TreeSpace`` and ``StackedSpace`` define it, and so may a + user's non-coordinate space, which is the case that matters: without it, + such a space can build no non-zero element at all and cannot be probed. + """ + + def ones(self) -> Any: + """Return the all-ones element of this space.""" + + +def probe_element(space: Space, ops: Any, dtype: Any) -> Any | None: + """Return a deterministic non-zero element, or ``None`` if none can be built. + + Prefers the space's own ``ones``; falls back to unflattening a flat ones + vector for a coordinate space. A bare :class:`VectorSpace` offers only + ``zeros``, which is useless as a probe, so it yields ``None`` and the caller + skips whatever check it wanted the element for. + """ + if isinstance(space, SupportsOnes): + return space.ones() + if isinstance(space, CoordinateSpace): + size = 1 + for dim in space.shape: + size *= int(dim) + return space.unflatten(ops.ones((size,), dtype=dtype)) + return None + + +# Capabilities every member of a dispatch registry already has, so they carry no +# information for class selection. Both registration and lookup go through +# ``registry_key`` so that adding a capability to ``_space_capabilities`` cannot +# silently stop the registries from matching. +_NON_DISPATCH_CAPABILITIES = frozenset({_CAP_COORDINATE}) + + +def registry_key(capabilities: CapabilitySet) -> CapabilitySet: + """Return the capability-dispatch key for a capability set.""" + return capabilities - _NON_DISPATCH_CAPABILITIES + + +def require(space: Space, capability: type, operation: str) -> None: + """Raise CapabilityError if ``space`` does not provide ``capability``.""" + if not isinstance(space, capability): + raise CapabilityError( + f"{operation} requires {capability.__name__}; " + f"{type(space).__name__} lacks this capability." + ) diff --git a/spacecore/space/concrete/_stacked.py b/spacecore/space/concrete/_stacked.py index a4b7677..89931e2 100644 --- a/spacecore/space/concrete/_stacked.py +++ b/spacecore/space/concrete/_stacked.py @@ -10,7 +10,12 @@ Space, StarSpace, ) -from ._tree_space import TreeSpace, _space_capabilities +from ._tree_space import TreeSpace +from ._dense_vector import EuclideanElementwiseJordanSpace, _validate_euclidean_elementwise_jordan +from .._capabilities import ( + _CAP_INNER, _CAP_STAR, _CAP_JORDAN, + _CAP_EUCLIDEAN_JORDAN, _space_capabilities, registry_key, +) from ..._checks import checked_method from ..._check_policy import CheckLevel from ...contextual import resolve_context_priority @@ -22,11 +27,6 @@ _STACKED_FALLBACK_ERRORS = (TypeError, ValueError, AttributeError, IndexError) -_CAP_INNER = InnerProductSpace -_CAP_STAR = StarSpace -_CAP_JORDAN = JordanAlgebraSpace -_CAP_EUCLIDEAN_JORDAN = EuclideanJordanAlgebraSpace - _STACKED_REGISTRY: dict[frozenset[type], type[StackedSpace]] = {} @@ -54,12 +54,14 @@ def _validate_count(count: int, owner: str = "StackedSpace") -> int: def _stacked_capabilities(base: Space) -> frozenset[type]: """Return capabilities copied from the stacked base space.""" + if isinstance(base, EuclideanElementwiseJordanSpace): + _validate_euclidean_elementwise_jordan(base, base.geometry) return _space_capabilities(base) def _stacked_class_for(capabilities: frozenset[type]) -> type[StackedSpace]: """Return the deterministic concrete class for a stacked capability set.""" - return _STACKED_REGISTRY.get(capabilities, StackedSpace) + return _STACKED_REGISTRY.get(registry_key(capabilities), StackedSpace) def _require_base(base: Space, capability: type, owner: str) -> None: @@ -452,7 +454,7 @@ class _StackedEuclideanJordanStarSpace( _STACKED_REGISTRY.update( - { + {registry_key(_caps): _cls for _caps, _cls in { frozenset(): StackedSpace, frozenset({_CAP_INNER}): _StackedInnerProductSpace, frozenset({_CAP_STAR}): _StackedStarSpace, @@ -467,7 +469,7 @@ class _StackedEuclideanJordanStarSpace( frozenset( {_CAP_INNER, _CAP_STAR, _CAP_JORDAN, _CAP_EUCLIDEAN_JORDAN} ): _StackedEuclideanJordanStarSpace, - } + }.items()} ) diff --git a/spacecore/space/concrete/_tree_space.py b/spacecore/space/concrete/_tree_space.py index b04a3d4..4e970d4 100644 --- a/spacecore/space/concrete/_tree_space.py +++ b/spacecore/space/concrete/_tree_space.py @@ -26,11 +26,11 @@ _validate_euclidean_elementwise_jordan, ) -CapabilitySet = frozenset[type] -_CAP_INNER = InnerProductSpace -_CAP_STAR = StarSpace -_CAP_JORDAN = JordanAlgebraSpace -_CAP_EUCLIDEAN_JORDAN = EuclideanJordanAlgebraSpace +from .._capabilities import ( + CapabilitySet, _CAP_INNER, _CAP_STAR, _CAP_JORDAN, registry_key, + _CAP_EUCLIDEAN_JORDAN, _space_capabilities, +) + _TREE_REGISTRY: dict[CapabilitySet, type[TreeSpace]] = {} @@ -63,26 +63,13 @@ def _validate_leaf_spaces( return spaces -def _space_capabilities(space: Space) -> CapabilitySet: - """Return structural capabilities advertised by one leaf space.""" - capabilities: set[type] = set() - if isinstance(space, InnerProductSpace): - capabilities.add(_CAP_INNER) - if isinstance(space, StarSpace): - capabilities.add(_CAP_STAR) - if isinstance(space, JordanAlgebraSpace): - capabilities.add(_CAP_JORDAN) - if isinstance(space, EuclideanJordanAlgebraSpace): - if isinstance(space, EuclideanElementwiseJordanSpace): - _validate_euclidean_elementwise_jordan(space, space.geometry) - capabilities.add(_CAP_EUCLIDEAN_JORDAN) - return frozenset(capabilities) - - def _tree_capabilities(spaces: Sequence[Space]) -> CapabilitySet: """Return capabilities shared by every tree leaf.""" if not spaces: return frozenset() + for space in spaces: + if isinstance(space, EuclideanElementwiseJordanSpace): + _validate_euclidean_elementwise_jordan(space, space.geometry) shared = set(_space_capabilities(spaces[0])) for space in spaces[1:]: shared.intersection_update(_space_capabilities(space)) @@ -330,7 +317,7 @@ def __new__( spaces = _validate_leaf_spaces(leaf_spaces) resolved_ctx = resolve_context_priority(ctx, *spaces) converted = tuple(space.convert(resolved_ctx) for space in spaces) - cls = _TREE_REGISTRY.get(_tree_capabilities(converted), TreeSpace) + cls = _TREE_REGISTRY.get(registry_key(_tree_capabilities(converted)), TreeSpace) return super(TreeSpace, cls).__new__(cls) def __init__( @@ -964,7 +951,7 @@ class _TreeEuclideanJordanStarSpace( _TREE_REGISTRY.update( - { + {registry_key(_caps): _cls for _caps, _cls in { frozenset(): TreeSpace, frozenset({_CAP_INNER}): TreeInnerProductSpace, frozenset({_CAP_STAR}): _TreeStarSpace, @@ -979,7 +966,7 @@ class _TreeEuclideanJordanStarSpace( frozenset( {_CAP_INNER, _CAP_STAR, _CAP_JORDAN, _CAP_EUCLIDEAN_JORDAN} ): _TreeEuclideanJordanStarSpace, - } + }.items()} ) # The capability-dispatch subclasses above register themselves: each derives from diff --git a/tests/linop/test_noncoordinate_endpoints.py b/tests/linop/test_noncoordinate_endpoints.py new file mode 100644 index 0000000..12f68cb --- /dev/null +++ b/tests/linop/test_noncoordinate_endpoints.py @@ -0,0 +1,217 @@ +"""Linear operators need coordinates only for coordinate operations.""" +import numpy as np +import pytest + +from spacecore import Context, NumpyOps, DenseCoordinateSpace +from spacecore._batching import _batched_inner +from spacecore import CapabilityError +from spacecore.linop import IdentityLinOp, MatrixFreeLinOp, ZeroLinOp +from spacecore.linop._algebra import _same_space_for_algebra +from spacecore.linop._metric import _metric_is_hermitian_by_basis, metric_rapply +from spacecore.space.base import CoordinateSpace, InnerProduct, InnerProductSpace, VectorSpace +from spacecore.space._capabilities import ( + _CAP_BATCH, _NON_DISPATCH_CAPABILITIES, _space_capabilities, probe_element, + registry_key, require, +) + + +class ArraySpace(VectorSpace): + """Array-backed vectors with no public coordinate or geometry surface.""" + + def zeros(self): + return self.ctx.asarray([0., 0.]) + + def add(self, x, y): + return x + y + + def scale(self, a, x): + return a * x + + def _convert(self, new_ctx): + return type(self)(new_ctx, check_level=self.check_level) + + +class DiagonalGeometry(InnerProduct): + def __init__(self, weights): + self.weights = np.asarray(weights) + + def inner(self, ops, x, y): + return ops.vdot(x, self.weights * y) + + def riesz(self, ops, x): + return self.weights * x + + def riesz_inverse(self, ops, x): + return x / self.weights + + +class MetricArraySpace(ArraySpace, InnerProductSpace): + def ones(self): + """Deterministic non-zero probe element; no coordinate surface needed.""" + return self.ctx.asarray([1., 1.]) + + def __init__(self, ctx, weights=(2., 5.), check_level=None): + super().__init__(ctx, check_level=check_level) + self.geometry = DiagonalGeometry(weights) + + def _eq_algebra(self, other): + return super()._eq_algebra(other) and np.array_equal( + self.geometry.weights, other.geometry.weights + ) + + def _convert(self, new_ctx): + return type(self)(new_ctx, self.geometry.weights, self.check_level) + + +@pytest.fixture(params=['none', 'standard', 'strict']) +def example(request): + ctx = Context(NumpyOps(), dtype=np.float64) + dom = MetricArraySpace(ctx, (2., 5.), request.param) + cod = MetricArraySpace(ctx, (3., 7.), request.param) + matrix = np.array([[1., 2.], [-3., 4.]]) + op = MatrixFreeLinOp( + lambda x: matrix @ x, lambda y: matrix.T @ y, dom, cod, + euclidean_adjoint=True, check_level=request.param, + ) + return op, matrix + + +def test_noncoordinate_algebra_and_metric_adjoint(example): + op, matrix = example + x, y = np.array([1., -2.]), np.array([3., 4.]) + expected = (matrix.T @ (op.codomain.geometry.weights * y)) / op.domain.geometry.weights + assert isinstance(op.domain, VectorSpace) + assert isinstance(op.domain, InnerProductSpace) + assert not isinstance(op.domain, CoordinateSpace) + for name in ('shape', 'size', 'flatten', 'unflatten', 'stacked'): + assert not hasattr(op.domain, name) + np.testing.assert_allclose(op.apply(x), matrix @ x) + np.testing.assert_allclose(op.rapply(y), expected) + np.testing.assert_allclose(op.H.apply(y), expected) + np.testing.assert_allclose(op.H.rapply(x), matrix @ x) + np.testing.assert_allclose((op + op).apply(x), 2 * (matrix @ x)) + np.testing.assert_allclose((op + op).rapply(y), 2 * expected) + np.testing.assert_allclose((3 * op).apply(x), 3 * (matrix @ x)) + np.testing.assert_allclose((op * 3).rapply(y), 3 * expected) + np.testing.assert_allclose((op.H @ op).apply(x), op.rapply(op.apply(x))) + np.testing.assert_allclose((op.H @ op).rapply(x), op.rapply(op.apply(x))) + np.testing.assert_allclose(op.codomain.inner(op.apply(x), y), op.domain.inner(x, op.rapply(y))) + assert not np.allclose(expected, matrix.T @ y) + assert _metric_is_hermitian_by_basis(op.H @ op) is None + + +@pytest.mark.parametrize('method', ['to_dense', 'to_matrix', 'vapply', 'rvapply']) +@pytest.mark.parametrize('variant', ['original', 'scaled', 'sum', 'composed', 'adjoint', 'zero', 'identity']) +def test_coordinate_operations_reject_noncoordinate_endpoints(example, method, variant): + op, _ = example + op = { + 'original': op, 'scaled': 2 * op, 'sum': op + op, + 'composed': op.H @ op, 'adjoint': op.H, + 'zero': ZeroLinOp(op.domain, op.codomain, check_level=op.check_level), + 'identity': IdentityLinOp(op.domain, check_level=op.check_level), + }[variant] + args = (np.ones((3, 2)),) if method in ('vapply', 'rvapply') else () + with pytest.raises(CapabilityError, match='CoordinateSpace'): + getattr(op, method)(*args) + + +@pytest.mark.parametrize('method', ['vapply', 'rvapply']) +def test_explicit_batch_callbacks_require_coordinates(example, method): + op, _ = example + batch_op = MatrixFreeLinOp( + op.apply, op.rapply, op.domain, op.codomain, + vapply=lambda xs: xs, rvapply=lambda ys: ys, euclidean_adjoint=False, + check_level=op.check_level, + ) + with pytest.raises(CapabilityError, match='CoordinateSpace'): + getattr(batch_op, method)(np.ones((3, 2))) + + +@pytest.mark.parametrize('missing_endpoint', ['domain', 'codomain']) +def test_plain_vector_space_forward_works_but_adjoint_requires_inner(missing_endpoint): + ctx = Context(NumpyOps(), dtype=np.float64) + plain, metric = ArraySpace(ctx), MetricArraySpace(ctx) + dom, cod = (plain, metric) if missing_endpoint == 'domain' else (metric, plain) + op = MatrixFreeLinOp(lambda x: 2 * x, lambda y: 2 * y, dom, cod) + x = np.array([1., 2.]) + np.testing.assert_array_equal(op.apply(x), 2 * x) + with pytest.raises(CapabilityError, match='InnerProductSpace'): + op.rapply(x) + with pytest.raises(CapabilityError, match='InnerProductSpace'): + metric_rapply(dom, cod, lambda y: y, x) + + +def test_require_names_space_capability_and_operation(): + space = ArraySpace(Context(NumpyOps(), dtype=np.float64)) + with pytest.raises(CapabilityError, match='probe.*CoordinateSpace.*ArraySpace'): + require(space, _CAP_BATCH, 'probe') + require(space, VectorSpace, 'add') + with pytest.raises(CapabilityError, match='CoordinateSpace'): + _batched_inner(space, np.ones((2, 2)), np.ones((2, 2))) + assert _space_capabilities(space) == frozenset() + assert CoordinateSpace in _space_capabilities(DenseCoordinateSpace((2,), space.ctx)) + + +def test_algebra_compatibility_does_not_inspect_shape(): + ctx = Context(NumpyOps(), dtype=np.float64) + assert not _same_space_for_algebra(MetricArraySpace(ctx, (2., 5.)), MetricArraySpace(ctx, (3., 7.))) + + +def test_strict_probe_runs_on_noncoordinate_space_and_catches_a_wrong_adjoint(): + """A space with `ones` is probed even without coordinates (finding 1).""" + ctx = Context(NumpyOps(), dtype=np.float64) + dom, cod = MetricArraySpace(ctx, (2., 5.)), MetricArraySpace(ctx, (3., 7.)) + matrix = np.array([[1., 2.], [-3., 4.]]) + assert not isinstance(dom, CoordinateSpace) + + # The coordinate transpose is *not* the metric adjoint on this geometry, and + # euclidean_adjoint=False asserts that it is. Strict checks must catch that. + with pytest.raises(ValueError, match='adjoint consistency check failed'): + MatrixFreeLinOp( + lambda x: matrix @ x, lambda y: matrix.T @ y, dom, cod, + euclidean_adjoint=False, check_level='strict', + ) + + # The wrapped coordinate adjoint is correct, so the same probe passes. + MatrixFreeLinOp( + lambda x: matrix @ x, lambda y: matrix.T @ y, dom, cod, + euclidean_adjoint=True, check_level='strict', + ) + + +def test_strict_probe_skipped_when_no_element_can_be_built(): + """A bare VectorSpace offers only zeros, so the probe is skipped, not crashed.""" + ctx = Context(NumpyOps(), dtype=np.float64) + + class NoOnes(MetricArraySpace): + ones = None # withdraw the capability + + dom = NoOnes(ctx, (2., 5.)) + assert probe_element(dom, dom.ctx.ops, dom.dtype) is None + MatrixFreeLinOp( + lambda x: 2. * x, lambda y: 2. * y, dom, dom, + euclidean_adjoint=False, check_level='strict', + ) + + +def test_registry_key_strips_non_dispatch_capabilities(): + """Registry lookup survives a capability that does not select a class (finding 2).""" + ctx = Context(NumpyOps(), dtype=np.float64) + leaf = DenseCoordinateSpace((2,), ctx) + caps = _space_capabilities(leaf) + + assert CoordinateSpace in caps + assert CoordinateSpace not in registry_key(caps) + assert registry_key(caps) == caps - _NON_DISPATCH_CAPABILITIES + # Dispatch still reaches the inner-product specialization rather than the base. + assert isinstance(leaf.stacked(3), InnerProductSpace) + + +def test_capability_error_is_public(): + """Callers catch it, so it is exported rather than reached through a private module.""" + import spacecore + import spacecore.space + + assert spacecore.CapabilityError is CapabilityError + assert spacecore.space.CapabilityError is CapabilityError + assert "CapabilityError" in spacecore.__all__