diff --git a/CHANGELOG.md b/CHANGELOG.md index 87158b947c..5e8d517699 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,61 @@ # Changelog +## Unreleased + +### Billing impact + +One change here raises prices, in one narrow case. **`full` and `full_like` +with a non-scalar `fill_value` no longer carry an inferred symmetry tag**, so +operations on their results are priced as the dense arrays they are — measured +at `+50%` on `fnp.sin` and `+55%` on `fnp.einsum("ij,ij->")` for a `(3, 3)` +result built from a length-3 fill. + +The rise reaches only calls that were being given a tag the data did not +support: a non-scalar fill into a shape with two or more equal-length axes. +A non-scalar fill into any other shape was never tagged and is unchanged +(a `(3, 4)` result still costs 384 on `fnp.sin`, a `(5,)` result still 160), +and a scalar `fill_value` keeps its tag and its price at every shape. Values +are unchanged throughout. + +**It ships as a new version opening a new phase, so no submission is +re-evaluated against it: nothing already scored is repriced.** + +### Symmetry + +- **A symmetry tag now matches the buffer it is attached to.** `as_symmetric` + and the public `SymmetricTensor` constructor validate within a tolerance, + but the cost model reads the resulting tag as exact. They now copy one + representative per orbit before attaching it, so the two agree. Data that is + already exactly invariant is passed through unchanged, so the common case + keeps its zero-copy view semantics — including its use as an `out=` + destination with a caller-chosen layout. Charges are unchanged. + +- **`symmetrize` gains `mode=`.** The default, `"reynolds-projection"`, is the + existing behaviour and is unchanged in both result and price. The new + `"canonical-copy"` keeps each orbit's lexicographically first entry instead + of averaging the orbit, bills `numel(data)` — the rate every other + materializing copy pays — and preserves the input dtype rather than + promoting to float64. It reads only the group's generators, so it stays + available for groups too large to enumerate. `fnp.random.symmetric` takes + the same argument. + +- **Reynolds symmetrization refuses an unenumerable group before charging for + it.** `symmetrize` and `fnp.random.symmetric` previously computed the cost + of a projection from a closed-form group order, charged it, and only then + hit the enumeration limit. They now raise `ValueError` first, naming + `mode="canonical-copy"`, and the refused call costs nothing. + +- **`full` and `full_like` no longer infer symmetry from shape when + `fill_value` is not a scalar.** The inferred group describes a constant + fill; a broadcast array fill writes distinct values into positions it would + otherwise report as redundant. + +- **Symmetry attachment inside the package routes through + `wrap_with_derived_symmetry`.** Trust is anchored to + `wrap_with_trusted_symmetry` alone, so `wrap_with_symmetry` — which nothing + internal calls — is validated and charged like any other caller-supplied + claim. + ## v0.12.0 (2026-08-21) ### Billing impact diff --git a/docs/reference/cost-model.md b/docs/reference/cost-model.md index c24c079825..fee260eb5b 100644 --- a/docs/reference/cost-model.md +++ b/docs/reference/cost-model.md @@ -859,7 +859,7 @@ each backed by a CI-enforced test you can open and read: | **No substitution arbitrage** | a bit-identical alias cannot bill cheaper than its canonical (e.g. `acos` *is* `arccos` — the 16× ufunc-alias fix); equivalent contractions (`dot`/`inner`/`matmul`/`einsum`) share one cost engine | `test_ufunc_alias_parity.py`, `test_random_weight_aliasing.py`; the shared einsum engine ([§Contraction](#contraction-einsum-family)) | | **No unpriceable or mispriceable dtype** | a non-numeric dtype (`dtype.kind` outside the allowlist `"biufc"` — object, string, bytes, structured/void, datetime64, timedelta64) is refused before any charge — as an operand, an explicit `dtype=`, or an `out=` destination — because it either has unbounded per-element cost (object) or a real per-element cost no flat rate captures (the rest); a dtype that is still zero-itemsize once NumPy materialises it is the one exception, since it carries no data either way — `'U0'`/`'S0'` are not, because NumPy promotes them to `'U1'`/`'S1'` on allocation | `tests/test_object_dtype_ban.py` | | **No cheap in-op path** | top-k `svd(k=)` cannot yield a *full* decomposition below full price (the `min(4mnk, economy)` cap + `k ≥ min → full` guard); invalid `k` (`< 1` or `> min(m, n)`) is rejected before any billing | `test_svd_topk_cost.py` (cap / guard / monotonicity); `test_linalg.py` (invalid-`k` `ValueError`) | -| **No free symmetry tag over a bare top-level constructor call** | a *bare, top-level* `SymmetricTensor(data, symmetry=…)` call cannot mint a tag over data that doesn't actually satisfy the claimed group: it now validates the claim through the same check `as_symmetric(data, symmetry=…)` already used, raising `SymmetryError` on a mismatch and charging `k·(7n − 1)` (`k` = non-identity generators, `n = data.size`) for a genuine one — and symmetry a flopscope op derives *internally* (e.g. `exp` propagating the tag of an already-validated operand, or a slice/transpose view) is exempt from re-paying, since it never reaches an unvalidated, caller-supplied claim. **Two narrower, in-process-only routes remain open and are NOT closed by this**: calling `flopscope._symmetry_utils.wrap_with_symmetry` directly (it only checks that the group's axes fit `ndim`, never buffer contents — trusted by this constructor for its ~30 legitimate internal uses in `_array_ops.py`), and constructing inside a participant callback a counted host op invokes (`fnp.apply_along_axis`, `fnp.piecewise`, …), which inherits the host op's trust. Neither is remotely reachable — `wrap_with_symmetry` is outside the server's op REGISTRY and unexported, and callback ops raise `RemoteCallbackError` on the server backend — so both are pinned as `strict=True` `xfail` regression tests rather than fixed, to avoid repricing the ~34 internal call sites that legitimately rely on the same trust | `tests/test_symmetric_tensor_new_validation.py` (incl. the two pinned `xfail`s); `tests/test_symmetric_cost.py` (`k·(7n−1)` rate) | +| **A symmetry tag matches the buffer it is attached to** | symmetry validation uses a tolerance (`np.allclose`, `atol=1e-6`, `rtol=1e-5`) while the cost model reads the tag as exact — it prices every position in an orbit after the first as a redundant degree of freedom and does not read the buffer again. The two untrusted ingress points, `as_symmetric(data, symmetry=…)` and a *bare, top-level* `SymmetricTensor(data, symmetry=…)`, therefore validate the claim (raising `SymmetryError` on a mismatch, charging `k·(7n − 1)` for a genuine one, `k` = non-identity generators, `n = data.size`) and then canonicalize: each orbit takes the value at its lexicographically smallest index, so values that survived only within the tolerance do not reach the tag. Data already exactly invariant is passed through unchanged, so the charge and the zero-copy view semantics are unaffected; the caller's array is never modified. Symmetry a flopscope op derives *internally* (e.g. `exp` propagating an already-validated operand's tag, or a slice/transpose view) is exempt, since it never carries a fresh caller-supplied claim. **Scope:** this is a property of those two ingress points, not a package-wide invariant — a Reynolds projection sums each orbit in a fixed element order, so its own output is typically invariant only to about an ulp, as is a symmetric matmul. That residue is rounding rather than caller-placed information, so it is sound for accounting; code needing bit-exactness should use an ingress point or `_canonical_symmetry.is_exactly_invariant`. In-process code can attach a tag regardless (`arr.view(SymmetricTensor)` plus an attribute assignment needs nothing from this package), so the boundary that holds is the wire: the server dispatches registered ops only, and none of those names is registered | `tests/test_symmetry_canonicalization.py`; `tests/test_symmetric_tensor_new_validation.py`; `tests/test_symmetric_cost.py` (`k·(7n−1)` rate) | | **Free-tier discipline** | weight 0 is limited to views/metadata, untouched (zero-page or uninitialized) allocation, and the narrow `astype`/`asarray` no-op (`copy=False` with an already-matching dtype; a dtype-free or dtype-matching `asarray`) — every other cast or copy, including a same-dtype `astype(copy=True)`, bills `numel` like `copy`. Any **metered** op that writes a new buffer — copied, replicated, constant-filled, gathered, or scattered — carries weight ≥ 1 (ndarray methods inherited from numpy are outside the meter by design — see [§The meter boundary](#the-meter-boundary)). Every value-test is charged wherever it hides: `a.nonzero()` (method), `where(1-arg)`, `argwhere`, `flatnonzero`, `count_nonzero` | `test_weight_tier_policy.py`; `test_data_movement_free_tier.py` (free-labels consistency guard) | | **No free-gather discount** | a computed-index gather (`take`, `take_along_axis`, `choose`) is metered at the access tier (weight 4.0) like any other non-sequential read, so precomputing a look-up table and then gathering from it no longer buys a categorical discount; only genuine view-indexing (a static/basic index, `arr[i]`) stays free | `test_data_movement_free_tier.py`; [§Copy and gather](#copy-and-gather) | | **Complex packing non-profitable** | folding two real payloads into one complex op's real/imag lanes bills the op's true complex structure (`multiply` factor 6, matmul exact `≈4.13×`), so the pack costs more than the honest real work it replaces | `tests/test_dtype_cost.py` (packing tests) | diff --git a/src/flopscope/__init__.py b/src/flopscope/__init__.py index 3f5df36c7d..a1a7fdfc9d 100644 --- a/src/flopscope/__init__.py +++ b/src/flopscope/__init__.py @@ -159,14 +159,22 @@ def reduction_cache_info(): def clear_cache() -> None: - """Clear all flopscope LRU caches (einsum + reduction). + """Clear all flopscope LRU caches (einsum, reduction, symmetry orbit maps). Convenience aggregate over :func:`einsum_clear_caches` and - :func:`reduction_clear_cache`. Use the per-domain variants if you - only need to invalidate one cache. + :func:`reduction_clear_cache`, plus the symmetry orbit-map cache. Use the + per-domain variants if you only need to invalidate one cache. + + The orbit-map cache is the one worth clearing for memory rather than + correctness: each entry holds one index per element of a distinct + ``(shape, group action)``, so it grows with the tensors it has been asked + about rather than with the number of operations. """ + from flopscope._canonical_symmetry import clear_canonical_map_cache + einsum_clear_caches() reduction_clear_cache() + clear_canonical_map_cache() def remote_unsupported_ops() -> frozenset[str]: diff --git a/src/flopscope/_array_ops.py b/src/flopscope/_array_ops.py index 51e27ecfeb..d112b1e96e 100644 --- a/src/flopscope/_array_ops.py +++ b/src/flopscope/_array_ops.py @@ -46,8 +46,8 @@ from flopscope._symmetry_utils import ( broadcast_group, validate_symmetry_group, + wrap_with_derived_symmetry, wrap_with_inferred_symmetry, - wrap_with_symmetry, wrap_with_trusted_symmetry, ) from flopscope._validation import _normalize_out, require_budget @@ -257,11 +257,17 @@ def full( _fill_probe_dtype = _np.asarray(fill_value).dtype refuse_non_numeric_dtype("full", _fill_probe_dtype) _billing_dtype = _np.dtype(dtype) if dtype is not None else _fill_probe_dtype + # The shape-derived group below says "every orbit holds one repeated + # value", which is true of a constant fill and false the moment + # fill_value broadcasts distinct values across the result. Decided here, + # outside the deduct, since it reads only the argument. + constant_fill = _np.ndim(fill_value) == 0 with budget.deduct( "full", flop_cost=cost, subscripts=None, shapes=(), dtypes=(_billing_dtype,) ): result = _call_numpy(_np.full, shape, fill_value, dtype=dtype, **kwargs) - result = _wrap_constant_fill(result) + if constant_fill: + result = _wrap_constant_fill(result) return result @@ -502,12 +508,20 @@ def full_like( result = _call_numpy( _np.full_like, _to_base_ndarray(a), fill_value, dtype=dtype, **kwargs ) + # full_like overwrites every element, so the template's symmetry says + # nothing about the result -- what carries the claim is the fill. A + # scalar fill leaves every orbit constant; a broadcast array fill writes + # distinct values into positions the tag would call redundant, so neither + # the propagated nor the shape-inferred group survives it. + constant_fill = _np.ndim(fill_value) == 0 propagated_symmetry = None - if isinstance(a, SymmetricTensor): + if constant_fill and isinstance(a, SymmetricTensor): propagated_symmetry = _compatible_symmetry_for_shape(a.symmetry, result.shape) if propagated_symmetry is not None: return wrap_with_trusted_symmetry(result, propagated_symmetry) # type: ignore[return-value] - inferred_symmetry = _infer_constant_shape_symmetry(result.shape) + inferred_symmetry = ( + _infer_constant_shape_symmetry(result.shape) if constant_fill else None + ) if inferred_symmetry is None: if isinstance(a, SymmetricTensor): return _np.array(result, copy=False, subok=False) # type: ignore[return-value] @@ -617,7 +631,7 @@ def reshape(a: ArrayLike, /, *args: Any, **kwargs: Any) -> FlopscopeArray: reason="reshape merges or splits axes inside the symmetric block", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -642,7 +656,7 @@ def transpose( out_group = _st.transport_transpose(in_group, ndim=a_arr.ndim, axes=axes) # transpose never genuinely drops (axis perm always preserves S_n etc.). if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -666,7 +680,7 @@ def swapaxes(a: ArrayLike, axis1: int, axis2: int) -> FlopscopeArray: axis2=axis2, ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -694,7 +708,7 @@ def moveaxis( destination=destination, ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -763,7 +777,7 @@ def concatenate( # keep writing to. return out # type: ignore[return-value] if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -812,7 +826,7 @@ def stack( if out is not None: return out # type: ignore[return-value] if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -850,7 +864,7 @@ def vstack(tup: Sequence[ArrayLike]) -> FlopscopeArray: reason="vstack breaks block symmetry", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -888,7 +902,7 @@ def hstack(tup: Sequence[ArrayLike]) -> FlopscopeArray: reason="hstack breaks block symmetry", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -930,7 +944,7 @@ def split( axis=axis, ) if out_group is not None: - return [wrap_with_symmetry(p, out_group) for p in raw_pieces] # type: ignore[return-value] + return [wrap_with_derived_symmetry(p, out_group) for p in raw_pieces] # type: ignore[return-value] return [_asplainflopscope(p) for p in raw_pieces] # type: ignore[return-value] @@ -961,7 +975,7 @@ def hsplit( reason="hsplit breaks block symmetry", ) if out_group is not None: - return [wrap_with_symmetry(p, out_group) for p in raw_pieces] # type: ignore[return-value] + return [wrap_with_derived_symmetry(p, out_group) for p in raw_pieces] # type: ignore[return-value] return [_asplainflopscope(p) for p in raw_pieces] # type: ignore[return-value] @@ -997,7 +1011,7 @@ def vsplit( ): raw_pieces = _call_numpy(_np.vsplit, ary_arr, indices_or_sections) if out_group is not None: - return [wrap_with_symmetry(p, out_group) for p in raw_pieces] # type: ignore[return-value] + return [wrap_with_derived_symmetry(p, out_group) for p in raw_pieces] # type: ignore[return-value] return [_asplainflopscope(p) for p in raw_pieces] # type: ignore[return-value] @@ -1028,7 +1042,7 @@ def squeeze( reason="squeeze removes an axis inside the symmetric block", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -1051,7 +1065,7 @@ def expand_dims(a: ArrayLike, axis) -> FlopscopeArray: axis=axis, ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -1090,7 +1104,7 @@ def ravel(a: ArrayLike, *args: Any, **kwargs: Any) -> FlopscopeArray: reason="ravel collapses to a single axis; block cannot fit", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -1120,7 +1134,7 @@ def copy(a: ArrayLike, *args: Any, **kwargs: Any) -> FlopscopeArray: ): result = _call_numpy(_np.copy, a_arr, *args, **kwargs) if isinstance(a, SymmetricTensor): - return wrap_with_symmetry(result, a.symmetry) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, a.symmetry) # type: ignore[return-value] return result # type: ignore[return-value] @@ -1226,7 +1240,7 @@ def tile(A: ArrayLike, reps: int | ArrayLike) -> FlopscopeArray: reason="tile reps not constant on block orbit", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -1259,7 +1273,7 @@ def repeat( reason="repeat along a block axis breaks block symmetry", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -1294,7 +1308,7 @@ def flip( reason="flip on a proper subset of block axes breaks group action", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -1327,7 +1341,7 @@ def roll( reason="roll along a block axis breaks block symmetry", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -1608,7 +1622,7 @@ def broadcast_to( reason="broadcast_to expands length-1 block axes", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -2014,7 +2028,7 @@ def _one(a): reason="atleast_1d incompatible with block structure", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) + return wrap_with_derived_symmetry(result, out_group) return _asplainflopscope(result) if len(arys) == 1: @@ -2050,7 +2064,7 @@ def _one(a): reason="atleast_2d incompatible with block structure", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) + return wrap_with_derived_symmetry(result, out_group) return _asplainflopscope(result) if len(arys) == 1: @@ -2086,7 +2100,7 @@ def _one(a): reason="atleast_3d incompatible with block structure", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) + return wrap_with_derived_symmetry(result, out_group) return _asplainflopscope(result) if len(arys) == 1: @@ -2252,7 +2266,7 @@ def broadcast_arrays(*args: ArrayLike, **kwargs: Any) -> tuple[FlopscopeArray, . input_shape=array.shape, output_shape=output_shape, ) - wrapped.append(wrap_with_symmetry(broadcasted, symmetry)) + wrapped.append(wrap_with_derived_symmetry(broadcasted, symmetry)) return tuple(wrapped) @@ -2369,7 +2383,7 @@ def column_stack(tup: Sequence[ArrayLike]) -> FlopscopeArray: reason="column_stack breaks block symmetry", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + return wrap_with_derived_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] @@ -2634,7 +2648,7 @@ def dsplit(ary: ArrayLike, *args: Any, **kwargs: Any) -> list[FlopscopeArray]: ): raw_pieces = _call_numpy(_np.dsplit, ary_arr, *args, **kwargs) if out_group is not None: - return [wrap_with_symmetry(p, out_group) for p in raw_pieces] # type: ignore[return-value] + return [wrap_with_derived_symmetry(p, out_group) for p in raw_pieces] # type: ignore[return-value] return [_asplainflopscope(p) for p in raw_pieces] # type: ignore[return-value] @@ -3377,7 +3391,12 @@ def matrix_transpose(x: ArrayLike) -> FlopscopeArray: reason="matrix_transpose: rank too low for sym", ) if out_group is not None: - return wrap_with_symmetry(result, out_group) # type: ignore[return-value] + # The trusted wrapper by name, because this function is the one + # symmetry-propagating transform without a @_counted_wrapper frame + # to inherit the exemption from (see the note above). Routing it + # through the ordinary derived helper would validate and charge, and + # this operation is registered as free. + return wrap_with_trusted_symmetry(result, out_group) # type: ignore[return-value] return _asplainflopscope(result) # type: ignore[return-value] diff --git a/src/flopscope/_canonical_symmetry.py b/src/flopscope/_canonical_symmetry.py new file mode 100644 index 0000000000..948ae4a866 --- /dev/null +++ b/src/flopscope/_canonical_symmetry.py @@ -0,0 +1,266 @@ +"""Canonical-copy: make an accepted symmetry claim exactly true before it is tagged. + +``as_symmetric`` validates a symmetry claim with ``np.allclose``, which accepts +"close enough" data. The cost model then treats every position in a symmetry +orbit as a redundant degree of freedom and never re-reads the buffer. Those two +facts disagree: a caller can scale an asymmetric tensor down until its orbit +differences fall under ``atol``, collect the tag, and scale back up through an +ordinary pointwise op -- recovering independent values in positions the cost +model has already priced as redundant. + +This module closes that gap at the trust boundary. Once tolerant validation +accepts a buffer, one representative per orbit is copied over the whole orbit. +The hidden values are destroyed before the tag is minted, which is what makes +the downstream symmetry discount honest without re-checking anything. + +The scope is worth stating exactly, because it is narrower than "a tag always +means exact invariance". What this module guarantees is a property of the +*ingress* points -- ``as_symmetric`` and the public ``SymmetricTensor`` +constructor, the two places a caller hands over a buffer the library has never +inspected. Symmetry that propagates algebraically through later operations is +trusted on the mathematics, not re-established here, and float arithmetic makes +that a genuinely weaker claim: the Reynolds projection sums each orbit in a +fixed element order, so its own output is typically invariant only to about an +ulp, and a symmetric matmul is the same. Those tags are sound for accounting -- +the residue is rounding, not information a caller placed there -- but code that +needs a buffer to be invariant *to the bit* must either come through an ingress +point or ask :func:`is_exactly_invariant`. + +The orbit map depends only on ``(shape, axes, generator action)`` -- never on +buffer contents -- so it is built once per distinct action and cached. Building +it walks the *generators*, not the group elements: enumerating ``|G|`` would +make ``as_symmetric`` cost as much as the Reynolds projection it deliberately +is not. +""" + +from __future__ import annotations + +import functools + +import numpy as np + +from flopscope._perm_group import SymmetryGroup + + +def _resolved_axes(group: SymmetryGroup) -> tuple[int, ...]: + """Tensor axes the group acts on, applying the same fallback as validation.""" + axes = group.axes + return tuple(axes) if axes is not None else tuple(range(group.degree)) + + +def _generator_fingerprint(group: SymmetryGroup) -> tuple: + """Hashable identity of the group's ACTION, without enumerating the group. + + ``SymmetryGroup.__hash__`` canonicalizes through ``elements()``, which runs + Dimino and can blow the enumeration budget -- exactly the cost this module + exists to avoid. The generator literals pin the action just as tightly for + caching purposes; two spellings of one group merely get two identical + cache entries. + """ + return ( + _resolved_axes(group), + group.degree, + tuple(tuple(gen.array_form) for gen in group.generators if not gen.is_identity), + ) + + +def _generator_images(shape: tuple[int, ...], axes, degree, gen_forms): + """Flat-index image of each generator, one vectorized pass per generator.""" + ndim = len(shape) + flat = np.arange(int(np.prod(shape)), dtype=np.intp).reshape(shape) + images = [] + for form in gen_forms: + perm = list(range(ndim)) + for i in range(degree): + perm[axes[i]] = axes[form[i]] + images.append(np.transpose(flat, perm).ravel()) + return images + + +def _build_canonical_map(shape: tuple[int, ...], fingerprint: tuple) -> np.ndarray: + """``map[i]`` = smallest C-order flat index in ``i``'s orbit. + + Min-label propagation over the generator action: each round pushes every + position's label down to the smallest label reachable in one generator + step, then pointer-jumps so labels reach orbit minima in log-many rounds. + Cost is ``O(N * r)`` per round with no Python-level loop over elements, + versus ``O(N * |G|)`` for an element enumeration. + """ + axes, degree, gen_forms = fingerprint + n = int(np.prod(shape)) + labels = np.arange(n, dtype=np.intp) + if not gen_forms: + labels.flags.writeable = False + return labels + + images = _generator_images(shape, axes, degree, gen_forms) + while True: + previous = labels + for image in images: + labels = np.minimum(labels, labels[image]) + labels = labels[labels] # pointer jumping + if np.array_equal(labels, previous): + break + + # Cached and shared across calls: never let a caller mutate it. + labels.flags.writeable = False + return labels + + +#: Ceiling on the memory the orbit-map cache may hold, in bytes. +#: An entry is one index per tensor element, so unlike flopscope's other LRUs +#: -- whose entries are small cost records -- entries here scale with the +#: tensors they describe: 8 MB for a 1024x1024 map, 34 MB at 2048x2048. A +#: plain entry count would therefore bound the number of maps while leaving +#: the footprint unbounded, and every distinct shape a caller asks about mints +#: a new one. 256 MB keeps the working set of a realistic estimator resident +#: while refusing to grow without limit. +_CANONICAL_MAP_CACHE_BYTES = 256 * 1024 * 1024 + + +class _OrbitMapCache: + """LRU over orbit maps, bounded by total bytes rather than entry count. + + Mirrors enough of ``functools.lru_cache``'s surface (``cache_info``, + ``cache_clear``) to be used and inspected the same way. + """ + + __slots__ = ("_entries", "_max_bytes", "_bytes", "_hits", "_misses") + + def __init__(self, max_bytes: int) -> None: + self._entries: dict[tuple, np.ndarray] = {} + self._max_bytes = max_bytes + self._bytes = 0 + self._hits = 0 + self._misses = 0 + + def get(self, shape: tuple[int, ...], fingerprint: tuple) -> np.ndarray: + key = (shape, fingerprint) + cached = self._entries.pop(key, None) + if cached is not None: + self._entries[key] = cached # refresh recency + self._hits += 1 + return cached + + self._misses += 1 + mapping = _build_canonical_map(shape, fingerprint) + # A single map larger than the whole budget is served but not kept, + # so one outsized request cannot evict everything and still not fit. + if mapping.nbytes <= self._max_bytes: + self._entries[key] = mapping + self._bytes += mapping.nbytes + # dicts iterate in insertion order and `get` reinserts on a hit, + # so the first key is the least recently used. + while self._bytes > self._max_bytes: + evicted = self._entries.pop(next(iter(self._entries))) + self._bytes -= evicted.nbytes + return mapping + + def cache_info(self): + return functools._CacheInfo( # type: ignore[attr-defined] + self._hits, self._misses, self._max_bytes, len(self._entries) + ) + + def cache_clear(self) -> None: + self._entries.clear() + self._bytes = 0 + self._hits = 0 + self._misses = 0 + + @property + def nbytes(self) -> int: + return self._bytes + + +#: Process-wide orbit-map cache. Named for the lookup it performs so that +#: ``cache_info()``/``cache_clear()`` read the same as flopscope's other LRUs. +_canonical_map_cached = _OrbitMapCache(_CANONICAL_MAP_CACHE_BYTES) + + +def canonical_map(shape: tuple[int, ...], group: SymmetryGroup) -> np.ndarray: + """Cached orbit map for ``(shape, group action)``. + + Returns a view rather than the cached array itself. NumPy lets a caller + re-enable the writeable flag on an array that owns its data, but not on + one whose base is read-only, and a map mutated in place would silently + mis-canonicalize every later call for the same shape and group. + """ + return _canonical_map_cached.get(tuple(shape), _generator_fingerprint(group)).view() + + +def is_exactly_invariant(array: np.ndarray, group: SymmetryGroup) -> bool: + """Whether every orbit already holds one repeated value, to the bit. + + Checking generators is enough: they generate the group, so a buffer fixed + by each generator is fixed by every element. This is the tolerance-free + twin of the ``allclose`` check validation runs -- equality, not closeness, + is precisely the property the tag is read as asserting. + + Answering this with ``==`` alone would be too generous by exactly one + value: ``-0.0 == 0.0`` is true, yet the two differ in a bit that + ``copysign`` reads straight back out. A sign bit sitting in a position the + cost model prices as redundant is information like any other, so zeros + that disagree in sign count as a difference here and send the buffer down + the copying path. + """ + array = np.asarray(array) + axes = _resolved_axes(group) + ndim = array.ndim + signed = array.dtype.kind in "fc" + for gen in group.generators: + if gen.is_identity: + continue + perm = list(range(ndim)) + for i in range(group.degree): + perm[axes[i]] = axes[gen.array_form[i]] + transposed = array.transpose(perm) + if not np.array_equal(array, transposed): + return False + if signed: + if not np.array_equal( + np.signbit(array.real), np.signbit(transposed.real) + ) or ( + array.dtype.kind == "c" + and not np.array_equal( + np.signbit(array.imag), np.signbit(transposed.imag) + ) + ): + return False + return True + + +def canonicalize(array: np.ndarray, group: SymmetryGroup) -> np.ndarray: + """Return data whose orbits are exactly constant, copying only if needed. + + Data that is already exactly invariant is returned untouched, so the + common case -- a genuinely symmetric buffer -- keeps ``as_symmetric``'s + zero-copy view semantics, including its use as an ``out=`` destination + with a caller-chosen memory layout. Only a buffer that merely passed the + tolerant check gets rewritten, which is exactly the case where the tag + would otherwise certify more than the data supports. + """ + array = np.asarray(array) + if array.size == 0 or is_exactly_invariant(array, group): + return array + return canonical_copy(array, group) + + +def canonical_copy(array: np.ndarray, group: SymmetryGroup) -> np.ndarray: + """Return a fresh array whose orbits each hold one representative value. + + The representative is the orbit's lexicographically smallest tensor index. + Advanced indexing gathers rather than computes, so the dtype survives + exactly -- unlike the Reynolds projection, which must upcast to average -- + and the result is a new buffer, so a caller's array is never mutated and + the tagged data can no longer be reached through the caller's alias. + """ + array = np.asarray(array) + if array.size == 0: + return array.copy() + mapping = canonical_map(array.shape, group) + return array.reshape(-1)[mapping].reshape(array.shape) + + +def clear_canonical_map_cache() -> None: + """Drop cached orbit maps (used by cache-management hooks and tests).""" + _canonical_map_cached.cache_clear() diff --git a/src/flopscope/_registry.py b/src/flopscope/_registry.py index af61edaa3f..de7388e564 100644 --- a/src/flopscope/_registry.py +++ b/src/flopscope/_registry.py @@ -2381,7 +2381,7 @@ "category": "counted_custom", "module": "flopscope", "complex_factor": 2.0, - "notes": "Reynolds projection onto a permutation group's invariant subspace. Cost: (|G|+1)*numel (|G| transposed adds + scaling pass; transpose/zeros free; validation uncounted).", + "notes": "Symmetrize onto a permutation group's invariant subspace. Cost depends on mode: mode='reynolds-projection' (default) is (|G|+1)*numel (|G| transposed adds + scaling pass; transpose/zeros free; validation uncounted); mode='canonical-copy' is numel (one write per element, copying each orbit's lexicographically-first entry over the orbit; orbit map built from generators and cached, never enumerates |G|; no validation pass, exact by construction). Reynolds upcasts to result_type(input, float64); canonical-copy preserves the input dtype.", }, "as_symmetric": { "category": "counted_custom", diff --git a/src/flopscope/_symmetric.py b/src/flopscope/_symmetric.py index ed9cfcfd47..73c796835c 100644 --- a/src/flopscope/_symmetric.py +++ b/src/flopscope/_symmetric.py @@ -7,6 +7,7 @@ import numpy as np from flopscope._budget import _counted_wrapper +from flopscope._canonical_symmetry import canonical_copy, canonicalize from flopscope._dtype_billing import integer_to_float64_min_dtype from flopscope._ndarray import FlopscopeArray, _asplainflopscope from flopscope._perm_group import SymmetryGroup @@ -19,12 +20,15 @@ remap_group_axes, restrict_group_to_axes, validate_symmetry_group, - wrap_with_symmetry, wrap_with_trusted_symmetry, ) from flopscope._validation import require_budget from flopscope._write_epoch import epoch_of -from flopscope.errors import SymmetryError +from flopscope.errors import ( + _SYMMETRY_DOCS_PATH, + SymmetryError, + _docs_url, +) # --------------------------------------------------------------------------- # Validation @@ -83,6 +87,46 @@ def _nonidentity_generator_count(group) -> int: return sum(1 for gen in group.generators if not gen.is_identity) +#: Symmetrization strategies accepted by :func:`symmetrize`'s ``mode``. +_SYMMETRIZE_MODES = frozenset({"reynolds-projection", "canonical-copy"}) + + +def _require_enumerable_for_reynolds(group) -> None: + """Refuse a Reynolds projection over a group that cannot be enumerated. + + Every other consumer of the enumeration budget can degrade to a dense + cost and carry on, because for them the group is only an accounting + detail. Reynolds averaging is the one place where enumerating the group + IS the computation, so there is nothing to degrade to -- and the + ``canonical-copy`` mode, which reads the generators alone, is the way + through. + """ + from flopscope._config import get_setting + from flopscope._perm_group import _DiminoBudgetExceeded + + budget = int(get_setting("dimino_budget")) # type: ignore[arg-type] + try: + order = group.order() + except _DiminoBudgetExceeded as exc: + seen, budget = exc.seen_count, exc.budget + else: + if order <= budget: + return + seen = order + raise ValueError( + f"symmetrize(mode='reynolds-projection') averages over every element " + f"of this symmetry group, and enumerating it exceeded dimino_budget " + f"({seen} > {budget}). Use mode='canonical-copy', which reads the " + f"group's generators alone and never enumerates it, billing " + f"numel(data) instead of (|G| + 1) * numel(data). It gives a " + f"different result, not a cheaper route to the same one: each orbit " + f"keeps its lexicographically first entry instead of averaging the " + f"orbit. (flops.configure(dimino_budget=...) raises the limit for " + f"in-process runs only.) " + f"See: {_docs_url(_SYMMETRY_DOCS_PATH)}" + ) from None + + def _project_core(array, group): """Raw Reynolds projection. UNCOUNTED. Returns an ndarray. @@ -123,20 +167,32 @@ def symmetrize( data: np.ndarray, *, symmetry, + mode: str = "reynolds-projection", ) -> SymmetricTensor: - """Project an array onto the invariant subspace of a permutation group. + """Make an array invariant under a permutation group. - This applies Reynolds symmetrization: + Two modes, differing in whether the discarded entries get a vote. + ``"reynolds-projection"`` (the default) averages each orbit: ``R_G(T) = (1 / |G|) * sum_{g in G} g · T`` + ``"canonical-copy"`` instead keeps one entry per orbit -- the one at the + lexicographically smallest index -- and copies it over the rest. Every + other value in the orbit is discarded rather than mixed in, which is what + makes it the right choice when the input is not trusted: nothing a caller + hid in the redundant positions can reach the result. It is also the + cheaper of the two, being one copy pass rather than ``|G|`` transposed + adds, and it never enumerates the group. + Parameters ---------- data : array_like - Input array to project. + Input array to symmetrize. symmetry : SymmetryGroup Symmetry group to average over. If ``symmetry.axes`` is ``None``, axes are interpreted as ``tuple(range(symmetry.degree))``. + mode : {"reynolds-projection", "canonical-copy"}, optional + Which symmetrization to apply. Defaults to ``"reynolds-projection"``. Returns ------- @@ -151,8 +207,8 @@ def symmetrize( Notes ----- - ``symmetrize`` performs exact Reynolds averaging internally, billing - ``(|G| + 1) * numel(data)`` FLOPs: + ``"reynolds-projection"`` performs exact Reynolds averaging internally, + billing ``(|G| + 1) * numel(data)`` FLOPs: - ``|G|`` transposed add passes over ``numel`` elements - one final scaling pass (divide by ``|G|``) @@ -161,6 +217,19 @@ def symmetrize( where ``|G|`` is the group order and ``numel = data.size``. + ``"canonical-copy"`` bills ``numel(data)`` -- one write per output + element, the same rate as every other materializing copy (``copy``, + ``take``, ``repeat``). It needs no validation pass because its output is + invariant by construction, and it never enumerates ``|G|``: the orbit map + is built from the group's generators and cached per + ``(shape, group action)``. + + The two modes also differ in dtype. Averaging must divide, so + ``"reynolds-projection"`` accumulates in + ``result_type(data, float64)`` -- a ``float32`` input comes back + ``float64``. ``"canonical-copy"`` only moves values, so it preserves the + input dtype exactly, including integer, boolean and complex types. + The canonical pattern for generating random data with symmetry is: ``fnp.random.symmetric(shape, symmetry_group, distribution=...)``. @@ -174,26 +243,55 @@ def symmetrize( >>> S.is_symmetric((0, 1)) True """ + if mode not in _SYMMETRIZE_MODES: + raise ValueError( + f"unknown symmetrize mode {mode!r}; " + f"expected one of {', '.join(map(repr, sorted(_SYMMETRIZE_MODES)))}" + ) array = np.asarray(data) group = _resolve_symmetry_argument(array, symmetry=symmetry) assert group is not None # required=True raises if symmetry is None validate_symmetry_group(group, ndim=array.ndim, shape=array.shape) n = array.size - cost = max((group.order() + 1) * n, 1) + if mode == "canonical-copy": + # One write per output element -- the rate every other materializing + # copy pays (copy/take/repeat). Deliberately does NOT consult + # group.order(): enumerating the group is the cost this mode exists + # to avoid, and the orbit map only ever needs the generators. + cost = max(n, 1) + # A gather moves values without computing any, so the output dtype is + # the input's -- no float64 sentinel, unlike the averaging branch. + dtypes = (array.dtype,) + else: + # Averaging visits every group element, so a group too large to + # enumerate cannot be projected at all. Refuse here, above the + # deduct, so a call that cannot finish is not charged for trying: + # for a group whose order is known in closed form the cost is + # computable, and without this check the caller would be billed the + # full Reynolds price and only then hit the enumeration limit. + _require_enumerable_for_reynolds(group) + cost = max((group.order() + 1) * n, 1) + # _project_core always accumulates in np.result_type(array, float64) -- + # the "/ group.order()" scaling pass needs float precision even from + # float32 input (verified: symmetrize(float32).dtype == float64, + # symmetrize(complex64).dtype == complex128) -- so the float64 sentinel + # must join the resolve rather than replace it (result_type preserves + # kind: result_type(complex64, float64) == complex128). + dtypes = (array.dtype, np.dtype(np.float64)) budget = require_budget() - # _project_core always accumulates in np.result_type(array, float64) -- - # the "/ group.order()" scaling pass needs float precision even from - # float32 input (verified: symmetrize(float32).dtype == float64, - # symmetrize(complex64).dtype == complex128) -- so the float64 sentinel - # must join the resolve rather than replace it (result_type preserves - # kind: result_type(complex64, float64) == complex128). with budget.deduct( "symmetrize", flop_cost=cost, subscripts=None, shapes=(array.shape,), - dtypes=(array.dtype, np.dtype(np.float64)), + dtypes=dtypes, ): + if mode == "canonical-copy": + # Exactly invariant by construction, so unlike the averaging + # branch there is no residual rounding to validate away. + return SymmetricTensor._construct_trusted( + canonical_copy(array, group), symmetry=group + ) projected = _project_core(array, group) # D1: internal validation runs but is NOT billed — build the tensor # directly rather than calling the (later-counted) as_symmetric. @@ -658,33 +756,39 @@ def _wrap_tensor_result(data: np.ndarray, symmetry: SymmetryGroup | None): return SymmetricTensor._construct_trusted(data, symmetry=symmetry) -# `wrap_with_symmetry`/`wrap_with_trusted_symmetry` (flopscope._symmetry_utils) -# are the trusted attachment path used ~30 times across _array_ops.py for -# symmetry DERIVED by an array transform (reshape/ravel/squeeze/split/...) -# rather than a fresh, caller-supplied claim, so this constructor treats -# either one, called directly, as trusted by code-object identity. +# Trust is anchored to ONE code object: `wrap_with_trusted_symmetry` +# (flopscope._symmetry_utils). A caller cannot fabricate it -- a copy of that +# function defined elsewhere compiles to a different code object, because +# `co_filename` participates in equality -- which is why trust is keyed on +# identity here rather than on an argument the caller could set. # -# KNOWN GAP (not closed by this fix): `wrap_with_symmetry` is an ordinary -# importable function -- `flopscope._symmetry_utils.wrap_with_symmetry` -- -# and its own check is structural only (do the group's axes fit `ndim`); it -# never inspects buffer contents. Code that imports it directly and calls -# `wrap_with_symmetry(asymmetric_data, fake_claim)` gets the exact same free, -# unvalidated tag this whole fix exists to close, just one function call -# away from the constructor. Pinned as an expected failure in -# `tests/test_symmetric_tensor_new_validation.py::test_wrap_with_symmetry_does_not_mint_an_unvalidated_tag`. -# Not remotely reachable: `wrap_with_symmetry` is absent from the server's op -# REGISTRY and is not exported on `flopscope` or `flopscope.numpy`, so this -# gap requires importing a private module directly -- consistent with this -# task's in-process-only, not-a-launch-blocker scope. Closing it properly -# needs either a capability token threaded through all ~34 internal call -# sites of `wrap_with_symmetry`/`wrap_with_trusted_symmetry` -# (_symmetry_utils.py + _array_ops.py) or a content check added inside -# `wrap_with_symmetry` itself -- deliberately NOT attempted here: repricing -# those 34 already-relied-upon call sites in a launch window is a larger risk -# than the (unreachable) hole they'd close. -_TRUSTED_SYMMETRY_WRAPPER_CODES = frozenset( - {wrap_with_symmetry.__code__, wrap_with_trusted_symmetry.__code__} -) +# The sibling wrappers -- `wrap_with_symmetry`, `wrap_with_derived_symmetry`, +# `wrap_with_inferred_symmetry` -- are deliberately NOT in here, and each +# constructs directly rather than delegating, so none of them can lend its +# caller this trust. The array transforms and constant fills that use them +# are exempt because they run inside a `@_counted_wrapper` frame, which +# `_called_from_wrapper` sees; imported and called on their own, from outside +# any flopscope op, they validate and charge like any other fresh claim. +# Delegating instead would hand every importer of those names a free tag, +# which is the hole this arrangement exists to keep shut. +# +# Two sites hold the trusted wrapper's remaining justification, both unable +# to inherit a counted frame: `_build_symmetric_proxy` +# (_accumulation/_cost.py), which tags an uninitialized `np.empty` scratch +# buffer the cost model only ever reads for shape and symmetry, and +# `matrix_transpose` (_array_ops.py), a registered-free operation that +# carries no `@_counted_wrapper` of its own. Worth knowing before anyone +# deletes this mechanism as dead weight. +# +# What this does NOT do is stop in-process code from minting a tag. It cannot: +# `arr.view(SymmetricTensor)` followed by an attribute assignment needs no +# helper from this package at all, and monkeypatching the validator works too. +# Narrowing the trust set is about keeping the package's own attachment sites +# honest and few, not about defending against a caller who already has the +# module. The boundary that holds is the wire -- the server dispatches +# registered ops only, and none of these names is registered -- which is the +# boundary `as_symmetric`'s canonicalization exists to serve. +_TRUSTED_SYMMETRY_WRAPPER_CODES = frozenset({wrap_with_trusted_symmetry.__code__}) class SymmetricTensor(FlopscopeArray): @@ -705,47 +809,31 @@ def __new__( *, symmetry: SymmetryGroup, ) -> SymmetricTensor: - # WHAT THIS CLOSES: a tag is a billing claim about buffer CONTENTS. - # Minting one here without validating granted the 1/|G| discount over - # arbitrary data, bypassing the validate-and-charge path - # `as_symmetric` already goes through. This now routes a BARE, - # TOP-LEVEL `SymmetricTensor(data, symmetry=...)` call through that - # SAME validator (rather than a second one that could disagree), so - # that specific call shape pays -- and is checked -- exactly like - # `as_symmetric`, and raises the identical `SymmetryError` for a - # false claim. + # A tag is a billing claim about buffer CONTENTS, so a bare, + # top-level `SymmetricTensor(data, symmetry=...)` is checked and + # charged exactly like `as_symmetric` -- same validator, same price, + # same `SymmetryError` on a false claim -- and then canonicalized, so + # the tag it mints asserts no more than the data supports. + # + # Validation is skipped on two trusted routes. First, an immediate + # caller of `wrap_with_trusted_symmetry`, the package's single + # trusted attachment point (see `_TRUSTED_SYMMETRY_WRAPPER_CODES` + # above). Second, construction from inside another flopscope op's + # `@_counted_wrapper` frame: many sites in this package (pointwise, + # einsum, solvers, random.symmetric, accumulation) tag a result whose + # symmetry they have already established mathematically -- exp() of a + # symmetric input really is symmetric -- and a cost-estimation proxy + # over uninitialized memory needs the metadata rather than a real + # check (see `_accumulation/_cost.py`'s `_build_symmetric_proxy`). # - # WHAT THIS DOES NOT CLOSE (two known, in-process-only, not-remotely- - # reachable gaps -- see `tests/test_symmetric_tensor_new_validation.py` - # for pinned `xfail` repros of both): - # 1. Calling `flopscope._symmetry_utils.wrap_with_symmetry` (or - # `wrap_with_trusted_symmetry`) directly: they are trusted by - # code-object identity below, and neither checks buffer - # contents. See the comment on `_TRUSTED_SYMMETRY_WRAPPER_CODES` - # just above this class for why that trust is not removed here. - # 2. Constructing from inside a participant callback that a - # counted host op invokes (e.g. `fnp.apply_along_axis`, - # `fnp.piecewise`): `_called_from_wrapper` walks the ENTIRE - # call stack for any `@_counted_wrapper` frame, so it cannot - # distinguish "genuinely inside the host op's own internal - # code" from "arbitrary caller-supplied code the host op - # happens to have called," and trusts both. - # Skip validation only when called from INSIDE another flopscope - # op's `@_counted_wrapper`-decorated frame: many call sites - # elsewhere in this package (pointwise, einsum, solvers, - # random.symmetric, accumulation -- all outside this file) tag a - # result whose symmetry they have already established mathematically - # (e.g. exp() of a symmetric input actually is symmetric; a cost- - # estimation proxy over uninitialized memory needs the metadata, not - # a real check -- see `_accumulation/_cost.py`'s - # `_build_symmetric_proxy`, which documents exactly this bypass). - # `_called_from_wrapper` (shared with `FlopscopeArray`'s protocol - # dispatch, in flopscope._budget) detects that nesting; it is False - # for ordinary top-level construction, which is this defect's - # exploit path, so that path still always validates and pays. The - # immediate-caller check below covers the sibling case one level - # up: `wrap_with_symmetry`/`wrap_with_trusted_symmetry` themselves - # calling this constructor directly, not nested inside a wrapper. + # KNOWN GAP (pinned as an expected failure in + # `tests/test_symmetric_tensor_new_validation.py`): `_called_from_wrapper` + # walks the ENTIRE call stack for a `@_counted_wrapper` frame, so a + # tensor constructed inside a participant callback that a counted host + # op invokes (`fnp.apply_along_axis`, `fnp.piecewise`) inherits the + # host's trust. Closing it needs the walk to stop at the callback + # boundary rather than pass through it. Not reachable on the graded + # backend, where those ops refuse a callback over the wire. from flopscope._budget import _called_from_wrapper array = np.asarray(input_array) @@ -755,6 +843,11 @@ def __new__( ) if not trusted: _validate_and_charge_symmetry(array, symmetry, op_name="as_symmetric") + # Same trust boundary as `as_symmetric`, so the same rule: the + # tolerant check authorizes a tag the cost model reads as exact, + # and canonicalizing here is what makes those two agree. Honest + # data is returned untouched, so this stays a view. + array = canonicalize(array, symmetry) obj = array.view(cls) obj._symmetry = symmetry obj._symmetry_inferred = False @@ -1027,7 +1120,12 @@ def as_symmetric( Returns ------- SymmetricTensor - View of ``data`` carrying validated symmetry metadata. + ``data`` carrying validated symmetry metadata, exactly invariant + under ``symmetry``. Data that is already exactly invariant is + wrapped as a view; data that satisfied the check only within the + tolerance is copied first, with each orbit taking the value at its + lexicographically smallest index, so the metadata describes the + buffer it is attached to. ``data`` itself is never modified. Raises ------ @@ -1048,7 +1146,18 @@ def as_symmetric( assert group is not None # required=True raises if symmetry is None array = np.asarray(data) _validate_and_charge_symmetry(array, group, op_name="as_symmetric") + # Validation is tolerant, but the tag it authorizes is read as exact: the + # cost model prices every orbit position after the first as redundant and + # never re-reads the buffer. Copy one representative across each orbit so + # that reading is true. Values that differed only within tolerance do not + # survive, which is what stops a caller from scaling an asymmetric tensor + # under atol, collecting the tag, and scaling back up with the independent + # values -- and their discount -- intact. Data that is already exactly + # invariant is passed through untouched, so the honest case keeps this + # function's zero-copy view semantics. + array = canonicalize(array, group) # Already validated and charged above -- construct via the trusted, # non-revalidating path so this doesn't pay (or re-check) twice through - # SymmetricTensor's public, validating constructor. + # SymmetricTensor's public, validating constructor. The canonical copy is + # exactly invariant by construction, so there is nothing left to re-check. return SymmetricTensor._construct_trusted(array, symmetry=group) diff --git a/src/flopscope/_symmetry_utils.py b/src/flopscope/_symmetry_utils.py index 8d63f48831..6d9a489e86 100644 --- a/src/flopscope/_symmetry_utils.py +++ b/src/flopscope/_symmetry_utils.py @@ -692,7 +692,16 @@ def reduce_group( def wrap_with_symmetry(data, symmetry: SymmetryGroup | None): - """Wrap ndarray-like data with symmetry metadata when a group is present.""" + """Attach a symmetry claim to data whose contents have NOT been checked. + + This is the untrusted spelling: it verifies only that the group's axes + fit the array's rank, then hands the buffer to the validating + constructor, which checks the contents and bills for doing so. Nothing + in this package calls it -- internal transforms use + :func:`wrap_with_derived_symmetry` -- so a caller reaching it is making + a fresh claim about data the library has never inspected, and pays the + same price as :func:`flopscope.as_symmetric` for the privilege. + """ array = np.asarray(data) if symmetry is None: return array @@ -703,16 +712,47 @@ def wrap_with_symmetry(data, symmetry: SymmetryGroup | None): def wrap_with_trusted_symmetry(data, symmetry: SymmetryGroup | None): - """Wrap data with already-proven symmetry metadata without re-validating. + """Attach symmetry metadata without validating or charging. + + The single trusted attachment point in the package: this function's code + object is the one ``SymmetricTensor.__new__`` recognizes, so trust is + anchored to something a caller cannot fabricate (a copy of this function + defined elsewhere gets a different code object and is not trusted). The + two wrappers below route through it rather than constructing directly, + which is what lets them inherit that trust without widening it. + + Only call this where the symmetry is already established: derived by an + algebraic transform of an already-tagged tensor, or correct by + construction. It never looks at the buffer. + """ + array = np.asarray(data) + if symmetry is None: + return array + from flopscope._symmetric import SymmetricTensor - This helper is for internal call sites only, where the symmetry was - generated or revalidated by trusted constructor logic. Avoiding the - redundant validation call keeps constructor hot paths fast while leaving - public/user-facing symmetry paths fully validated. + return SymmetricTensor(array, symmetry=symmetry) + + +def wrap_with_derived_symmetry(data, symmetry: SymmetryGroup | None): + """Attach symmetry carried over from an already-tagged input. + + For array transforms -- reshape, transpose, split, concatenate and the + rest -- whose output symmetry a ``_symmetry_transport`` helper computed + from the input's own validated group. The claim is inherited rather than + fresh, so it is not re-checked or re-billed; the structural check below + only confirms the transported group still fits the new rank. + + Constructs directly rather than delegating to + :func:`wrap_with_trusted_symmetry`, which matters: this function's own + code object is NOT trusted, so the transforms get their exemption from + running inside a counted op rather than from calling this helper. Imported + and called on its own, from outside any flopscope op, it validates and + charges like any other fresh claim. """ array = np.asarray(data) if symmetry is None: return array + validate_symmetry_group(symmetry, ndim=array.ndim) from flopscope._symmetric import SymmetricTensor return SymmetricTensor(array, symmetry=symmetry) @@ -721,11 +761,15 @@ def wrap_with_trusted_symmetry(data, symmetry: SymmetryGroup | None): def wrap_with_inferred_symmetry(data, symmetry: SymmetryGroup | None): """Wrap data with auto-inferred symmetry metadata. - Identical to :func:`wrap_with_trusted_symmetry` except the resulting + Identical to :func:`wrap_with_derived_symmetry` except the resulting array carries ``_symmetry_inferred = True``. Read by ``_prepare_symmetric_out`` to decide whether a non-symmetric ``out=`` write should silently downgrade the target (inferred) or raise (explicit). Internal call sites only — never expose to user code. + + Constructs directly for the same reason as the helper above: the + constant-fill sites that use it run inside counted ops, so an imported + direct call earns no exemption. """ array = np.asarray(data) if symmetry is None: diff --git a/src/flopscope/numpy/random/__init__.py b/src/flopscope/numpy/random/__init__.py index b2b892ffe6..629b1131c5 100644 --- a/src/flopscope/numpy/random/__init__.py +++ b/src/flopscope/numpy/random/__init__.py @@ -663,6 +663,8 @@ def symmetric( shape: int | Sequence[int], symmetry: SymmetryGroup, distribution: str | Callable[..., Any] = "randn", + *, + mode: str = "reynolds-projection", **distribution_kwargs: Any, ) -> FlopscopeArray: """Sample random data and project it to a symmetry group. @@ -681,6 +683,12 @@ def symmetric( - ``size=shape`` and returns an array. + mode : {"reynolds-projection", "canonical-copy"}, optional + How the sample is made invariant, with the same meaning as + :func:`flopscope.symmetrize`'s ``mode``. The default averages each + orbit; ``"canonical-copy"`` keeps one entry per orbit instead, reads + the group's generators alone, and so remains available for groups too + large to enumerate. **distribution_kwargs Extra keyword arguments forwarded to the distribution function. @@ -692,7 +700,10 @@ def symmetric( Raises ------ ValueError - If ``shape`` is not an integer or a tuple/list of integers. + If ``shape`` is not an integer or a tuple/list of integers, if + ``mode`` is not one of the two accepted values, or if + ``mode="reynolds-projection"`` is asked to average over a group too + large to enumerate within ``dimino_budget``. TypeError If ``distribution`` is neither a NumPy random distribution name nor a callable. @@ -788,10 +799,32 @@ def symmetric( "distribution must be a numpy random function name or a callable" ) + from flopscope._symmetric import ( + _SYMMETRIZE_MODES, + SymmetricTensor, + _project_core, + _require_enumerable_for_reynolds, + validate_symmetry_groups, + ) + + if mode not in _SYMMETRIZE_MODES: + raise ValueError( + f"unknown symmetric mode {mode!r}; " + f"expected one of {', '.join(map(repr, sorted(_SYMMETRIZE_MODES)))}" + ) + budget = require_budget() - G = _builtins.max(symmetry.order(), 1) - # sample numel + projection core ((|G|+1)*numel) == sample + symmetrize - cost = _builtins.max(sample_size + (G + 1) * sample_size, 1) + if mode == "canonical-copy": + # sample numel + one copy pass, matching + # symmetrize(mode="canonical-copy")'s numel rate. + cost = _builtins.max(sample_size + sample_size, 1) + else: + # Refuse above the deduct so a projection that cannot finish is not + # billed for trying; mirrors symmetrize's own guard. + _require_enumerable_for_reynolds(symmetry) + G = _builtins.max(symmetry.order(), 1) + # sample numel + projection core ((|G|+1)*numel) == sample + symmetrize + cost = _builtins.max(sample_size + (G + 1) * sample_size, 1) # `sample` is usually drawn from a named real-only numpy distribution, # but `distribution` may be an arbitrary caller callable that returns # complex data; the registry's complex_factor="illegal" here is only @@ -804,12 +837,12 @@ def symmetric( shapes=(shape_tuple,), dtypes=(), ): - from flopscope._symmetric import ( - SymmetricTensor, - _project_core, - validate_symmetry_groups, - ) + if mode == "canonical-copy": + from flopscope._canonical_symmetry import canonical_copy + # Exactly invariant by construction, so no validation pass. + # canonical_copy does its own np.asarray. + return SymmetricTensor(canonical_copy(sample, symmetry), symmetry=symmetry) projected = _project_core(sample, symmetry) # _project_core does np.asarray validate_symmetry_groups(projected, [symmetry]) # uncounted safety check return SymmetricTensor(projected, symmetry=symmetry) diff --git a/tests/test_contraction_label_budget.py b/tests/test_contraction_label_budget.py index 9ca5bc9513..18815f9bc2 100644 --- a/tests/test_contraction_label_budget.py +++ b/tests/test_contraction_label_budget.py @@ -810,7 +810,10 @@ def test_untagged_operand_above_budget_still_pays_the_dense_price( ) assert dense == 2 * alpha - m - tagged = billed(lambda: fn(_sym_ones(a_big, (0, 1)), plain_b)) + # Built outside billed(...) so only the contraction is measured: tagging + # is a validated, charged operation in its own right. + sym_a = _sym_ones(a_big, (0, 1)) + tagged = billed(lambda: fn(sym_a, plain_b)) assert tagged < dense # the symmetry genuinely scales this contraction diff --git a/tests/test_symmetric_tensor_new_validation.py b/tests/test_symmetric_tensor_new_validation.py index 1ebdb536e8..25a308a9cd 100644 --- a/tests/test_symmetric_tensor_new_validation.py +++ b/tests/test_symmetric_tensor_new_validation.py @@ -54,15 +54,16 @@ def test_genuinely_symmetric_data_still_works(): assert t.symmetry is not None -@pytest.mark.xfail( - reason="known in-process bypass: wrap_with_symmetry runs only a structural " - "check, and __new__ trusts it by code-object identity. Not remotely " - "reachable (absent from REGISTRY and from the flopscope/fnp namespaces). " - "Closing it needs a capability token threaded through 34 internal call " - "sites, or a content check inside wrap_with_symmetry itself.", - strict=True, -) def test_wrap_with_symmetry_does_not_mint_an_unvalidated_tag(): + """The untrusted wrapper is not a side door around the validating constructor. + + ``wrap_with_symmetry`` used to be trusted by code-object identity, so a + direct call attached a tag to unexamined data for free. Trust now belongs + to ``wrap_with_trusted_symmetry`` alone, and the package's own transforms + reach it through ``wrap_with_derived_symmetry``; nothing internal calls + this one. A caller that does is making a fresh claim and is checked and + charged for it like any other. + """ from flopscope._symmetry_utils import wrap_with_symmetry raw = np.random.default_rng(0).random((6, 6)) diff --git a/tests/test_symmetry_canonicalization.py b/tests/test_symmetry_canonicalization.py new file mode 100644 index 0000000000..7c03b11e57 --- /dev/null +++ b/tests/test_symmetry_canonicalization.py @@ -0,0 +1,805 @@ +"""A symmetry tag must assert no more than the buffer supports. + +Validation accepts data that is symmetric *within a tolerance*, but the cost +model reads the tag it grants as exact: every position in an orbit after the +first is priced as a redundant degree of freedom, and the buffer is never read +again. Those two readings disagree for any buffer whose orbit entries merely +agree closely, and the gap is wide enough to carry independent values through +it -- scale a tensor down until its differences fall under ``atol``, collect +the tag, scale back up. + +Canonicalizing at the boundary is what makes the two readings agree: one entry +per orbit survives, so whatever was hidden in the others is gone before the tag +exists. These tests pin that property (exact invariance, zero tolerance), the +representative rule that decides which entry survives, and the two things the +fix must not cost -- the Reynolds path's behaviour, and a copy on data that is +already exact. +""" + +import numpy as np +import pytest + +import flopscope as flops +import flopscope.numpy as fnp +from flopscope._canonical_symmetry import ( + _canonical_map_cached, + canonical_copy, + is_exactly_invariant, +) +from flopscope._perm_group import SymmetryGroup, _Permutation +from flopscope._symmetric import SymmetricTensor + +BUDGET = 10**14 + +# A power of two, so scaling down and back up is exact in binary floating +# point and cannot itself be blamed for any difference we observe. +SCALE = 2.0**-40 + + +def _budget(): + return flops.BudgetContext(flop_budget=BUDGET, quiet=True) + + +def _generator_perms(group, ndim): + """Full-rank transpose permutation for each non-identity generator.""" + axes = group.axes if group.axes is not None else tuple(range(group.degree)) + perms = [] + for gen in group.generators: + if gen.is_identity: + continue + perm = list(range(ndim)) + for i in range(group.degree): + perm[axes[i]] = axes[gen.array_form[i]] + perms.append(tuple(perm)) + return perms + + +def _assert_exactly_invariant(array, group): + """Invariant under every generator with NO tolerance whatsoever.""" + raw = np.asarray(array) + for perm in _generator_perms(group, raw.ndim): + assert np.array_equal(raw, np.transpose(raw, perm)), ( + f"not exactly invariant under generator permutation {perm}" + ) + + +def _custom_group(axes=(0, 1, 2)): + """A group defined by raw generators rather than a named constructor.""" + return SymmetryGroup( + _Permutation([1, 2, 0]), _Permutation([1, 0, 2]), axes=tuple(axes) + ) + + +GROUPS = [ + pytest.param(SymmetryGroup.symmetric(axes=(0, 1)), (5, 5), id="symmetric-S2"), + pytest.param(SymmetryGroup.symmetric(axes=(0, 1, 2)), (4, 4, 4), id="symmetric-S3"), + pytest.param(SymmetryGroup.cyclic(axes=(0, 1, 2)), (4, 4, 4), id="cyclic-C3"), + pytest.param( + SymmetryGroup.dihedral(axes=(0, 1, 2, 3)), (3, 3, 3, 3), id="dihedral-D4" + ), + pytest.param(_custom_group(), (4, 4, 4), id="custom-generators"), + pytest.param( + SymmetryGroup.symmetric(axes=(1, 2)), (2, 4, 4), id="symmetric-inner-axes" + ), +] + + +# --------------------------------------------------------------------------- +# The defect itself +# --------------------------------------------------------------------------- + + +class TestToleranceGapIsClosed: + def test_scaled_down_asymmetric_data_is_accepted_but_not_kept(self): + """The tolerant check still admits it; the tag no longer over-claims. + + The pre-fix behaviour was that both halves of this test's premise + held at once: validation passed AND the orbit still held two + different values. Only the first may survive. + """ + group = SymmetryGroup.symmetric(axes=(0, 1)) + raw = np.array([[1.0, 2.0], [9.0, 3.0]]) + + # The claim is false at full scale, and is refused. + with _budget(), pytest.raises(flops.SymmetryError): + flops.as_symmetric(raw, symmetry=group) + + # Scaled down, the same claim passes the tolerance policy unchanged. + scaled = raw * SCALE + assert np.allclose(scaled, scaled.T, atol=1e-6, rtol=1e-5) + with _budget(): + tagged = flops.as_symmetric(scaled, symmetry=group) + + values = np.asarray(tagged) + assert values[0, 1] == values[1, 0], ( + "orbit still holds two different values under a tag read as exact" + ) + _assert_exactly_invariant(values, group) + + def test_scaling_back_up_cannot_recover_the_hidden_value(self): + """The round trip the defect depended on now yields nothing.""" + group = SymmetryGroup.symmetric(axes=(0, 1)) + raw = np.array([[1.0, 2.0], [9.0, 3.0]]) + + with _budget(): + tagged = flops.as_symmetric(raw * SCALE, symmetry=group) + restored = tagged * (1.0 / SCALE) + + assert not np.array_equal(np.asarray(restored), raw) + _assert_exactly_invariant(restored, group) + + @pytest.mark.parametrize("group,shape", GROUPS) + def test_round_trip_stays_exact_for_every_group(self, group, shape): + rng = np.random.default_rng(11) + asymmetric = rng.standard_normal(shape) * SCALE + + with _budget(): + tagged = flops.as_symmetric(asymmetric, symmetry=group) + restored = tagged * (1.0 / SCALE) + + _assert_exactly_invariant(tagged, group) + _assert_exactly_invariant(restored, group) + + def test_bare_constructor_is_the_same_boundary(self): + """``SymmetricTensor(...)`` is ingress too, and canonicalizes as well.""" + from flopscope._symmetric import SymmetricTensor + + group = SymmetryGroup.symmetric(axes=(0, 1)) + with _budget(): + tagged = SymmetricTensor( + np.array([[1.0, 2.0], [9.0, 3.0]]) * SCALE, symmetry=group + ) + _assert_exactly_invariant(tagged, group) + + @pytest.mark.parametrize( + "wrapper_name", + [ + "wrap_with_symmetry", + "wrap_with_derived_symmetry", + "wrap_with_inferred_symmetry", + ], + ) + def test_importable_wrappers_cannot_mint_a_free_tag(self, wrapper_name): + """Only one wrapper is trusted, and these are not it. + + Each of these is exempt from validation when it runs inside a counted + op, which is where the package uses them. Called directly, from + outside any flopscope op, they must be checked and charged -- so none + of them is a way around the constructor for anyone who imports it. + """ + import flopscope._symmetry_utils as symmetry_utils + + wrapper = getattr(symmetry_utils, wrapper_name) + group = SymmetryGroup.symmetric(axes=(0, 1)) + asymmetric = np.random.default_rng(3).random((6, 6)) + with _budget(), pytest.raises(flops.SymmetryError): + wrapper(asymmetric, group) + + def test_trust_is_anchored_to_one_code_object(self): + """The exemption belongs to one function, not to a growing list.""" + from flopscope._symmetric import _TRUSTED_SYMMETRY_WRAPPER_CODES + from flopscope._symmetry_utils import wrap_with_trusted_symmetry + + assert _TRUSTED_SYMMETRY_WRAPPER_CODES == frozenset( + {wrap_with_trusted_symmetry.__code__} + ) + + def test_in_process_code_can_always_mint_a_tag_and_that_is_not_the_boundary(self): + """Where the boundary actually is, written down so it is not mistaken. + + Code running in this process can attach a symmetry tag to anything it + likes, and no arrangement of private helpers changes that. The route + below uses only public NumPy and the class object -- no flopscope + helper is involved -- so hardening ``wrap_with_trusted_symmetry`` or + ``_construct_trusted`` would close two doors in a wall that does not + surround anything. Monkeypatching the validator works equally well. + + The boundary that does hold is the wire: a submission runs against a + server that dispatches registered operations only, and none of the + names used here is registered. That is what the canonicalization at + ``as_symmetric`` protects -- a claim arriving over that wire -- and it + is why these routes are documented rather than chased. + """ + from flopscope._registry import REGISTRY + from flopscope._symmetric import SymmetricTensor + from flopscope._symmetry_utils import wrap_with_trusted_symmetry + + group = SymmetryGroup.symmetric(axes=(0, 1)) + asymmetric = np.random.default_rng(5).random((6, 6)) + + # Needs nothing from flopscope but the class itself. + forged = asymmetric.view(SymmetricTensor) + forged._symmetry = group + assert forged.symmetry is not None + assert not is_exactly_invariant(np.asarray(forged), group) + + # The two private helpers are no different in kind, and no worse. + with _budget() as budget: + wrap_with_trusted_symmetry(asymmetric, group) + SymmetricTensor._construct_trusted(asymmetric, symmetry=group) + assert budget.flops_used == 0 + + # None of it crosses the wire. + for name in ( + "SymmetricTensor", + "wrap_with_trusted_symmetry", + "wrap_with_derived_symmetry", + "view", + ): + assert name not in REGISTRY + assert not hasattr(flops, "wrap_with_trusted_symmetry") + assert not hasattr(fnp, "wrap_with_trusted_symmetry") + + def test_matrix_transpose_stays_free(self): + """The registered-free transform must not start paying to keep its tag.""" + group = SymmetryGroup.symmetric(axes=(0, 1)) + base = np.random.default_rng(7).standard_normal((16, 16)) + with _budget() as budget: + tagged = flops.as_symmetric((base + base.T) / 2, symmetry=group) + before = budget.flops_used + transposed = fnp.matrix_transpose(tagged) + assert budget.flops_used == before + assert isinstance(transposed, SymmetricTensor) + assert transposed.symmetry is not None + + +# --------------------------------------------------------------------------- +# Which value survives +# --------------------------------------------------------------------------- + + +class TestCanonicalRepresentative: + def test_keeps_the_first_entry_rather_than_averaging(self): + """The documented rule, and the one case that tells the modes apart.""" + group = SymmetryGroup.symmetric(axes=(0, 1)) + raw = np.array([[1.0, 2.0], [9.0, 3.0]]) + + with _budget(): + copied = flops.symmetrize(raw, symmetry=group, mode="canonical-copy") + + np.testing.assert_array_equal( + np.asarray(copied), np.array([[1.0, 2.0], [2.0, 3.0]]) + ) + # Emphatically not the Reynolds answer, whose off-diagonal is 5.5. + assert np.asarray(copied)[0, 1] != 5.5 + + def test_representative_is_the_lexicographically_smallest_index(self): + """Every orbit resolves to its smallest C-order flat index.""" + group = SymmetryGroup.symmetric(axes=(0, 1, 2)) + shape = (3, 3, 3) + flat = np.arange(int(np.prod(shape)), dtype=np.float64).reshape(shape) + + with _budget(): + copied = np.asarray( + flops.symmetrize(flat, symmetry=group, mode="canonical-copy") + ) + + for index in np.ndindex(shape): + orbit_min = min( + int(flat[tuple(index[p] for p in perm)]) + for perm in [ + (0, 1, 2), + (0, 2, 1), + (1, 0, 2), + (1, 2, 0), + (2, 0, 1), + (2, 1, 0), + ] + ) + assert copied[index] == orbit_min + + def test_discarded_entries_cannot_influence_the_result(self): + """Two inputs differing only off-representative agree afterwards.""" + group = SymmetryGroup.symmetric(axes=(0, 1)) + a = np.array([[1.0, 2.0], [9.0, 3.0]]) + b = np.array([[1.0, 2.0], [-500.0, 3.0]]) + + with _budget(): + ca = flops.symmetrize(a, symmetry=group, mode="canonical-copy") + cb = flops.symmetrize(b, symmetry=group, mode="canonical-copy") + + np.testing.assert_array_equal(np.asarray(ca), np.asarray(cb)) + + +# --------------------------------------------------------------------------- +# What must not change +# --------------------------------------------------------------------------- + + +class TestReynoldsUnchanged: + def test_default_mode_still_averages(self): + group = SymmetryGroup.symmetric(axes=(0, 1)) + raw = np.array([[1.0, 2.0], [9.0, 3.0]]) + with _budget(): + projected = flops.symmetrize(raw, symmetry=group) + np.testing.assert_allclose( + np.asarray(projected), np.array([[1.0, 5.5], [5.5, 3.0]]) + ) + + def test_default_matches_explicit_reynolds_in_value_and_price(self): + group = SymmetryGroup.symmetric(axes=(0, 1)) + raw = np.random.default_rng(5).standard_normal((6, 6)) + + with _budget() as implicit: + a = flops.symmetrize(raw, symmetry=group) + with _budget() as explicit: + b = flops.symmetrize(raw, symmetry=group, mode="reynolds-projection") + + np.testing.assert_array_equal(np.asarray(a), np.asarray(b)) + assert implicit.flops_used == explicit.flops_used + + def test_reynolds_bills_the_group_order_and_canonical_copy_bills_numel(self): + group = SymmetryGroup.symmetric(axes=(0, 1)) + raw = np.ones((8, 8)) + + with _budget() as reynolds: + flops.symmetrize(raw, symmetry=group) + with _budget() as copied: + flops.symmetrize(raw, symmetry=group, mode="canonical-copy") + with _budget() as plain_copy: + fnp.copy(fnp.asarray(raw)) + + # Both scale by the same dtype rate, so their ratio is the model's. + assert reynolds.flops_used == copied.flops_used * (group.order() + 1) + # Copying one entry per orbit is a copy, and is priced as one. + assert copied.flops_used == plain_copy.flops_used + + def test_unknown_mode_is_refused(self): + with _budget(), pytest.raises(ValueError, match="unknown symmetrize mode"): + flops.symmetrize( + np.ones((4, 4)), + symmetry=SymmetryGroup.symmetric(axes=(0, 1)), + mode="nonsense", + ) + + +class TestAsSymmetricBilling: + @pytest.mark.parametrize("group,shape", GROUPS) + def test_price_does_not_depend_on_whether_canonicalization_ran(self, group, shape): + """Enforcement is the library's own business, not a charge to the caller.""" + rng = np.random.default_rng(17) + base = rng.standard_normal(shape) + # Exactly invariant: canonicalization short-circuits. + with _budget(): + exact = np.asarray( + flops.symmetrize(base, symmetry=group, mode="canonical-copy") + ) + # Only tolerantly invariant: canonicalization copies. + approximate = rng.standard_normal(shape) * SCALE + + with _budget() as no_copy: + flops.as_symmetric(exact, symmetry=group) + with _budget() as with_copy: + flops.as_symmetric(approximate, symmetry=group) + + assert no_copy.flops_used == with_copy.flops_used + + def test_caller_buffer_is_never_mutated(self): + group = SymmetryGroup.symmetric(axes=(0, 1)) + raw = np.array([[1.0, 2.0], [9.0, 3.0]]) * SCALE + before = raw.copy() + with _budget(): + flops.as_symmetric(raw, symmetry=group) + np.testing.assert_array_equal(raw, before) + + def test_inexact_input_detaches_from_the_caller_buffer(self): + """The consequence of copy-on-inexact, stated so it is not a surprise. + + A tag built from data that only passed the tolerant check is a copy, + so writing through it -- via ``out=`` -- no longer reaches the array + the caller handed in. This is the price of the tag being true, and it + applies only to inexact input; the exact case keeps its alias, which + the neighbouring test pins. + """ + group = SymmetryGroup.symmetric(axes=(0, 1)) + rng = np.random.default_rng(53) + destination = rng.standard_normal((8, 8)) * SCALE + assert np.allclose(destination, destination.T, atol=1e-6, rtol=1e-5) + assert not is_exactly_invariant(destination, group) + original = destination.copy() + + with _budget(): + tagged = flops.as_symmetric(destination, symmetry=group) + source = flops.as_symmetric(np.ones((8, 8)), symmetry=group) + fnp.exp(source, out=tagged) + + assert not np.shares_memory(np.asarray(tagged), destination) + np.testing.assert_array_equal(destination, original) + + def test_exactly_symmetric_input_is_not_copied(self): + """The honest case keeps its zero-copy view, and its memory layout.""" + group = SymmetryGroup.symmetric(axes=(0, 1)) + rng = np.random.default_rng(23) + base = rng.standard_normal((6, 6)) + exact = (base + base.T) / 2 + assert is_exactly_invariant(exact, group) + + with _budget(): + tagged = flops.as_symmetric(exact, symmetry=group) + + assert np.shares_memory(np.asarray(tagged), exact) + + +# --------------------------------------------------------------------------- +# Dtypes +# --------------------------------------------------------------------------- + + +DTYPES = [ + np.float32, + np.float64, + np.int32, + np.int64, + np.uint8, + np.bool_, + np.complex64, + np.complex128, +] + + +class TestDtypes: + @pytest.mark.parametrize("dtype", DTYPES) + def test_canonical_copy_preserves_dtype_exactly(self, dtype): + """Unlike averaging, copying needs no promotion.""" + group = SymmetryGroup.symmetric(axes=(0, 1)) + raw = np.array([[1, 2], [9, 3]], dtype=dtype) + + with _budget(): + copied = flops.symmetrize(raw, symmetry=group, mode="canonical-copy") + + assert np.asarray(copied).dtype == np.dtype(dtype) + + @pytest.mark.parametrize("dtype", DTYPES) + def test_canonicalized_result_is_exactly_invariant(self, dtype): + group = SymmetryGroup.symmetric(axes=(0, 1)) + raw = np.array([[1, 2], [9, 3]], dtype=dtype) + with _budget(): + copied = flops.symmetrize(raw, symmetry=group, mode="canonical-copy") + _assert_exactly_invariant(copied, group) + + def test_reynolds_still_promotes_where_it_must(self): + """The averaging branch keeps its float64 accumulation.""" + group = SymmetryGroup.symmetric(axes=(0, 1)) + raw = np.array([[1.0, 2.0], [9.0, 3.0]], dtype=np.float32) + with _budget(): + projected = flops.symmetrize(raw, symmetry=group) + assert np.asarray(projected).dtype == np.float64 + + +# --------------------------------------------------------------------------- +# The orbit map and its cache +# --------------------------------------------------------------------------- + + +class TestOrbitMapCache: + def test_same_shape_different_action_do_not_share_a_map(self): + """The key must separate groups that act differently on one shape.""" + shape = (4, 4, 4) + rng = np.random.default_rng(29) + raw = rng.standard_normal(shape) + + symmetric = canonical_copy(raw, SymmetryGroup.symmetric(axes=(0, 1, 2))) + cyclic = canonical_copy(raw, SymmetryGroup.cyclic(axes=(0, 1, 2))) + + assert not np.array_equal(symmetric, cyclic) + _assert_exactly_invariant(symmetric, SymmetryGroup.symmetric(axes=(0, 1, 2))) + _assert_exactly_invariant(cyclic, SymmetryGroup.cyclic(axes=(0, 1, 2))) + + def test_same_group_different_axes_do_not_share_a_map(self): + shape = (4, 4, 4) + rng = np.random.default_rng(31) + raw = rng.standard_normal(shape) + + on_01 = canonical_copy(raw, SymmetryGroup.symmetric(axes=(0, 1))) + on_12 = canonical_copy(raw, SymmetryGroup.symmetric(axes=(1, 2))) + + assert not np.array_equal(on_01, on_12) + _assert_exactly_invariant(on_01, SymmetryGroup.symmetric(axes=(0, 1))) + _assert_exactly_invariant(on_12, SymmetryGroup.symmetric(axes=(1, 2))) + + def test_same_action_different_shape_do_not_share_a_map(self): + group = SymmetryGroup.symmetric(axes=(0, 1)) + small = canonical_copy(np.arange(9.0).reshape(3, 3), group) + large = canonical_copy(np.arange(16.0).reshape(4, 4), group) + assert small.shape == (3, 3) + assert large.shape == (4, 4) + _assert_exactly_invariant(small, group) + _assert_exactly_invariant(large, group) + + def test_repeated_calls_reuse_the_cached_map(self): + group = SymmetryGroup.symmetric(axes=(0, 1)) + raw = np.ones((7, 7)) + canonical_copy(raw, group) + before = _canonical_map_cached.cache_info() + for _ in range(20): + canonical_copy(raw, group) + after = _canonical_map_cached.cache_info() + assert after.misses == before.misses + assert after.hits > before.hits + + def test_map_is_not_writeable(self): + """A shared cached map must not be mutable through a caller's handle.""" + from flopscope._canonical_symmetry import canonical_map + + mapping = canonical_map((4, 4), SymmetryGroup.symmetric(axes=(0, 1))) + with pytest.raises(ValueError): + mapping[0] = 3 + + def test_writeable_flag_cannot_be_re_enabled(self): + """Read-only is not enough on its own. + + NumPy lets a caller flip ``writeable`` back on for an array that owns + its data, so handing out the cached array itself would leave the map + editable in place -- and a doctored map mis-canonicalizes every later + call for that shape and group, which is the whole fix undone quietly. + Callers get a view, whose base refuses the flag. + """ + from flopscope._canonical_symmetry import canonical_map + + mapping = canonical_map((8, 8), SymmetryGroup.symmetric(axes=(0, 1))) + assert not mapping.flags.owndata + with pytest.raises(ValueError): + mapping.flags.writeable = True + + def test_clear_cache_drains_the_orbit_maps(self): + """The public aggregate must reach this cache; entries are array-sized.""" + from flopscope._canonical_symmetry import canonical_map + + canonical_map((9, 9), SymmetryGroup.symmetric(axes=(0, 1))) + assert _canonical_map_cached.cache_info().currsize > 0 + flops.clear_cache() + assert _canonical_map_cached.cache_info().currsize == 0 + + def test_cache_is_bounded_by_bytes_not_entry_count(self): + """Entries scale with the tensors they describe, so a count is no bound. + + ``symmetrize(mode="canonical-copy")`` is a registered op, so a caller + can mint a map for any shape it likes. Bounding entries rather than + bytes would cap the number of maps while leaving the footprint free to + grow with the shapes requested. + """ + from flopscope._canonical_symmetry import ( + _CANONICAL_MAP_CACHE_BYTES, + canonical_map, + ) + + group = SymmetryGroup.symmetric(axes=(0, 1)) + flops.clear_cache() + for side in range(600, 640): + canonical_map((side, side), group) + assert _canonical_map_cached.nbytes <= _CANONICAL_MAP_CACHE_BYTES + flops.clear_cache() + + def test_a_map_larger_than_the_whole_budget_is_served_but_not_kept(self): + """One outsized request must not evict everything and still not fit.""" + from flopscope._canonical_symmetry import ( + _CANONICAL_MAP_CACHE_BYTES, + _generator_fingerprint, + _OrbitMapCache, + ) + + tiny = _OrbitMapCache(max_bytes=8) # smaller than any real map + group = SymmetryGroup.symmetric(axes=(0, 1)) + mapping = tiny.get((4, 4), _generator_fingerprint(group)) + assert mapping.nbytes > 8 + assert tiny.cache_info().currsize == 0 + assert tiny.nbytes == 0 + # And the served map is still correct. + assert np.array_equal( + mapping, canonical_copy(np.arange(16.0).reshape(4, 4), group).ravel() + ) + assert _CANONICAL_MAP_CACHE_BYTES > 0 + + def test_eviction_does_not_change_results(self): + """A map rebuilt after eviction must equal the one that was dropped.""" + from flopscope._canonical_symmetry import _generator_fingerprint, _OrbitMapCache + + group = SymmetryGroup.symmetric(axes=(0, 1)) + fingerprint = _generator_fingerprint(group) + reference = np.array( + _OrbitMapCache(max_bytes=10**9).get((6, 6), fingerprint), copy=True + ) + + cache = _OrbitMapCache(max_bytes=400) # room for roughly one map + first = np.array(cache.get((6, 6), fingerprint), copy=True) + for side in (7, 8, 9): + cache.get((side, side), fingerprint) # force eviction + rebuilt = cache.get((6, 6), fingerprint) + + assert np.array_equal(first, reference) + assert np.array_equal(rebuilt, reference) + + +class TestSignedZero: + """``-0.0 == 0.0``, but ``copysign`` tells them apart. + + A sign bit in a position the cost model prices as redundant is + information, so equality alone is too generous a test for "already + exact": the buffer has to go down the copying path. + """ + + def test_signed_zero_is_not_treated_as_already_invariant(self): + group = SymmetryGroup.symmetric(axes=(0, 1)) + mixed = np.array([[1.0, 0.0], [-0.0, 2.0]]) + assert np.array_equal(mixed, mixed.T) # `==` cannot see it + assert not is_exactly_invariant(mixed, group) + + def test_tagging_removes_the_sign_bit_from_the_orbit(self): + group = SymmetryGroup.symmetric(axes=(0, 1)) + with _budget(): + tagged = flops.as_symmetric( + np.array([[1.0, 0.0], [-0.0, 2.0]]), symmetry=group + ) + signs = np.signbit(np.asarray(tagged)) + assert np.array_equal(signs, signs.T), ( + "sign bit survived in a position priced as redundant" + ) + + def test_copysign_cannot_read_an_asymmetry_back_out(self): + """The end-to-end route: the recovered signs must be symmetric.""" + group = SymmetryGroup.symmetric(axes=(0, 1)) + rng = np.random.default_rng(101) + base = np.zeros((8, 8)) + # Scatter negative zeros asymmetrically through the buffer. + mask = rng.random((8, 8)) < 0.5 + base[mask] = -0.0 + with _budget(): + tagged = flops.as_symmetric(base, symmetry=group) + recovered = fnp.copysign(np.ones((8, 8)), tagged) + values = np.asarray(recovered) + assert np.array_equal(values, values.T) + + def test_complex_signed_zero_is_caught_on_both_components(self): + group = SymmetryGroup.symmetric(axes=(0, 1)) + mixed = np.array([[1 + 0j, complex(0.0, 0.0)], [complex(-0.0, -0.0), 2 + 0j]]) + assert np.array_equal(mixed, mixed.T) + assert not is_exactly_invariant(mixed, group) + + def test_ordinary_symmetric_data_still_short_circuits(self): + """The check must not have become so strict that nothing passes.""" + group = SymmetryGroup.symmetric(axes=(0, 1)) + rng = np.random.default_rng(103) + base = rng.standard_normal((16, 16)) + exact = (base + base.T) / 2 + assert is_exactly_invariant(exact, group) + with _budget(): + tagged = flops.as_symmetric(exact, symmetry=group) + assert np.shares_memory(np.asarray(tagged), exact) + + +# --------------------------------------------------------------------------- +# Cost of the boundary +# --------------------------------------------------------------------------- + + +class TestTrustedPropagationStaysCheap: + def test_downstream_ops_do_not_rebuild_or_reapply_the_map(self): + """Canonicalize once at ingress; propagate algebraically thereafter.""" + group = SymmetryGroup.symmetric(axes=(0, 1)) + rng = np.random.default_rng(37) + base = rng.standard_normal((8, 8)) + exact = (base + base.T) / 2 + + with _budget(): + tagged = flops.as_symmetric(exact, symmetry=group) + before = _canonical_map_cached.cache_info() + for _ in range(25): + fnp.exp(tagged) + fnp.multiply(tagged, 2.0) + tagged.T # noqa: B018 - exercising the transpose propagation + tagged[0:8] + after = _canonical_map_cached.cache_info() + + assert (after.hits, after.misses) == (before.hits, before.misses) + + def test_propagated_results_still_carry_their_tag(self): + group = SymmetryGroup.symmetric(axes=(0, 1)) + base = np.random.default_rng(41).standard_normal((6, 6)) + with _budget(): + tagged = flops.as_symmetric((base + base.T) / 2, symmetry=group) + propagated = fnp.exp(tagged) + assert isinstance(propagated, SymmetricTensor) + assert propagated.symmetry is not None + + +# --------------------------------------------------------------------------- +# Groups too large to enumerate +# --------------------------------------------------------------------------- + + +class TestOversizedGroups: + def _oversized(self): + return SymmetryGroup.symmetric(axes=tuple(range(9))), (2,) * 9 + + def test_reynolds_refuses_and_names_the_alternative(self): + group, shape = self._oversized() + with _budget(), pytest.raises(ValueError, match="canonical-copy"): + flops.symmetrize(np.zeros(shape), symmetry=group) + + def test_refusal_costs_nothing(self): + """An operation that cannot finish must not be billed for trying.""" + group, shape = self._oversized() + with _budget() as budget: + with pytest.raises(ValueError): + flops.symmetrize(np.zeros(shape), symmetry=group) + assert budget.flops_used == 0 + + def test_canonical_copy_still_works_there(self): + group, shape = self._oversized() + raw = np.random.default_rng(43).standard_normal(shape) + with _budget() as budget: + copied = flops.symmetrize(raw, symmetry=group, mode="canonical-copy") + with _budget() as plain_copy: + fnp.copy(fnp.asarray(raw)) + _assert_exactly_invariant(copied, group) + assert budget.flops_used == plain_copy.flops_used + + def test_random_symmetric_refuses_and_offers_the_same_way_through(self): + group, shape = self._oversized() + with _budget(), pytest.raises(ValueError, match="canonical-copy"): + fnp.random.symmetric(shape, group) + + with _budget(): + sample = fnp.random.symmetric(shape, group, mode="canonical-copy") + _assert_exactly_invariant(sample, group) + + def test_random_symmetric_default_is_unchanged(self): + group = SymmetryGroup.symmetric(axes=(0, 1)) + with _budget() as implicit: + fnp.random.symmetric((6, 6), group) + with _budget() as explicit: + fnp.random.symmetric((6, 6), group, mode="reynolds-projection") + assert implicit.flops_used == explicit.flops_used + + def test_random_symmetric_rejects_an_unknown_mode(self): + with _budget(), pytest.raises(ValueError, match="unknown symmetric mode"): + fnp.random.symmetric((4, 4), SymmetryGroup.symmetric(axes=(0, 1)), mode="x") + + +# --------------------------------------------------------------------------- +# Non-constant fills +# --------------------------------------------------------------------------- + + +class TestNonScalarFillIsNotSymmetric: + """A fill that varies across the array leaves no orbit constant. + + ``full``/``full_like`` infer symmetry from shape alone, which describes a + constant fill. Broadcasting an array through them writes different values + into positions the inferred tag would call redundant. + """ + + def test_full_with_array_fill_carries_no_tag(self): + with _budget(): + result = fnp.full((3, 3), np.array([1.0, 2.0, 3.0])) + assert getattr(result, "symmetry", None) is None + + def test_full_like_with_array_fill_carries_no_tag(self): + with _budget(): + template = fnp.zeros((3, 3)) + # A constant fill is genuinely symmetric, so the template is tagged. + assert isinstance(template, SymmetricTensor) + result = fnp.full_like(template, np.array([1.0, 2.0, 3.0])) + assert getattr(result, "symmetry", None) is None + + def test_array_fill_is_priced_like_untagged_data(self): + fill = np.array([1.0, 2.0, 3.0]) + with _budget(): + tagged = fnp.full_like(fnp.zeros((3, 3)), fill) + + def _billed(fn): + with _budget() as budget: + fn() + return budget.flops_used + + honest = _billed(lambda: fnp.sin(fnp.asarray(np.broadcast_to(fill, (3, 3))))) + assert _billed(lambda: fnp.sin(tagged)) == honest + + @pytest.mark.parametrize("value", [0.0, 3.5, -1]) + def test_scalar_fill_keeps_its_legitimate_tag(self, value): + with _budget(): + full = fnp.full((4, 4), value) + full_like = fnp.full_like(fnp.zeros((4, 4)), value) + assert isinstance(full, SymmetricTensor) + assert isinstance(full_like, SymmetricTensor) + _assert_exactly_invariant(full, full.symmetry) diff --git a/website/content/docs/guides/symmetry.mdx b/website/content/docs/guides/symmetry.mdx index d78174e531..4a20a46fbb 100644 --- a/website/content/docs/guides/symmetry.mdx +++ b/website/content/docs/guides/symmetry.mdx @@ -52,8 +52,18 @@ with flops.BudgetContext(flop_budget=10**6) as budget: print(budget.flops_used) ``` -`flops.as_symmetric()` validates the data first. After that, Flopscope propagates -symmetry metadata algebraically through many operations. Unary pointwise ops preserve symmetry-aware costs and keep the same exact group, including non-full groups such as `C_k` or `D_k`. Slicing, reductions, and binary pointwise ops can weaken it or remove it entirely. +`flops.as_symmetric()` validates the data first. Validation uses a tolerance, +while the tag it grants is read as exact — the cost model prices every position +in an orbit after the first as redundant and does not read the buffer again — so +`as_symmetric` makes the two agree before tagging: each orbit takes the value at +its lexicographically smallest index. Data that is already exactly invariant is +returned unchanged, so a genuinely symmetric array is tagged without being +copied; data that passed only within the tolerance is copied first, and the +values in the discarded positions do not survive. Your own array is never +modified either way. + +After that, Flopscope propagates symmetry metadata algebraically through many +operations. Unary pointwise ops preserve symmetry-aware costs and keep the same exact group, including non-full groups such as `C_k` or `D_k`. Slicing, reductions, and binary pointwise ops can weaken it or remove it entirely. ## How to declare symmetry @@ -175,12 +185,31 @@ This helper is ideal for docs, tests, and experiments: - prefer `fnp.random.symmetric()` for synthetic data generation - `fnp.random.symmetric` internally samples data and calls ``flops.symmetrize`` so the projection and validation behavior is identical. -- approximate costs (meaningful estimate): +- approximate costs (meaningful estimate), for the default + ``mode="reynolds-projection"``: - `fnp.random.symmetric`: ``C_dist(n) + |G| * n + n`` + validation - `flops.symmetrize`: ``|G| * n + n`` + validation with ``n`` total elements and ``|G|`` group order. - in exact arithmetic it projects onto the invariant subspace, and in practice `flops.as_symmetric()` validates the result with its usual validation tolerances + +Both functions also accept ``mode="canonical-copy"``, which keeps one entry per +orbit — the one at the lexicographically smallest index — instead of averaging +the orbit. It replaces the ``|G| * n + n`` projection term with a single ``n`` +copy pass (so ``flops.symmetrize`` costs ``n``, and ``fnp.random.symmetric`` +costs ``C_dist(n) + n``), preserves the input dtype instead of promoting to +float64, and reads only the group's generators, so it stays available for groups +too large to enumerate (see below). It is a different result, not a cheaper +route to the same one: the other entries in each orbit are discarded rather than +averaged in. Reach for it when you want a tensor that *has* the symmetry, and +for Reynolds when you want the projection of your data *onto* the invariant +subspace. + +```python +x = fnp.array([[1.0, 2.0], [9.0, 3.0]]) +flops.symmetrize(x, symmetry=(0, 1)) # [[1, 5.5], [5.5, 3]] +flops.symmetrize(x, symmetry=(0, 1), mode="canonical-copy") # [[1, 2 ], [2, 3]] +``` - it keeps examples consistent across symmetry classes ```python @@ -535,6 +564,22 @@ trip the cap. The warning fires once per `(op_name, |G|)` pair per process to avoid log flooding. Suppress with `flops.configure(symmetry_warnings=False)`, which shares the flag with `SymmetryLossWarning`. +`flops.symmetrize` and `fnp.random.symmetric` are the exception to the +degrade-and-continue behaviour above: for them, enumerating the group *is* the +computation, so there is no dense cost to fall back to. Above the budget they +raise `ValueError` rather than warn, and the refused call is not charged. Use +`mode="canonical-copy"`, which reads the generators alone and never enumerates +the group: + +```python +big = flops.SymmetryGroup.symmetric(axes=tuple(range(9))) # |G| = 362,880 +wide = fnp.zeros((2,) * 9) +with flops.BudgetContext(flop_budget=int(1e9)): + flops.symmetrize(wide, symmetry=big) # ValueError + flops.symmetrize(wide, symmetry=big, mode="canonical-copy") # fine + fnp.random.symmetric((2,) * 9, big, mode="canonical-copy") # fine +``` + The gate keys on `|G|` rather than degree because `|G|` is the actual driver of enumeration cost. A high-degree but small-order group like `SymmetryGroup.cyclic(axes=tuple(range(50)))` (degree 50, `|G| = 50`) passes diff --git a/website/public/ops.json b/website/public/ops.json index f8f9d6b5e4..39b377ca5c 100644 --- a/website/public/ops.json +++ b/website/public/ops.json @@ -10766,11 +10766,11 @@ "free": false, "module": "flopscope", "name": "symmetrize", - "notes": "Reynolds projection onto a permutation group's invariant subspace. Cost: (|G|+1)*numel (|G| transposed adds + scaling pass; transpose/zeros free; validation uncounted).", + "notes": "Symmetrize onto a permutation group's invariant subspace. Cost depends on mode: mode='reynolds-projection' (default) is (|G|+1)*numel (|G| transposed adds + scaling pass; transpose/zeros free; validation uncounted); mode='canonical-copy' is numel (one write per element, copying each orbit's lexicographically-first entry over the orbit; orbit map built from generators and cached, never enumerates |G|; no validation pass, exact by construction). Reynolds upcasts to result_type(input, float64); canonical-copy preserves the input dtype.", "numpy_ref": "np.symmetrize", "slug": "symmetrize", "status": "supported", - "summary": "Project an array onto the invariant subspace of a permutation group.", + "summary": "Make an array invariant under a permutation group.", "weight": 1.0 }, {