Symmetry-aware einsum cost rewrite + de-vendor opt_einsum

.gitlint

@@ -8,4 +8,4 @@ ignore-merge-commits=true
 types=feat,fix,docs,style,refactor,perf,test,build,ci,chore,revert
 
 [title-max-length]
-line-length=80
+line-length=90

CHANGELOG.md

@@ -2,6 +2,133 @@
 
 ## Unreleased
 
+### Changed (BREAKING)
+
+- **Vendored opt_einsum replaced with runtime dependency.** flopscope now
+  depends on `opt_einsum>=3.3.0,<4.0.0` instead of vendoring its source.
+  The remaining `flopscope._opt_einsum` is a slim ~830-line shim that
+  adapts upstream's `PathInfo` to flopscope's expected shape and recomputes
+  per-step FLOP costs using flopscope's FMA convention.
+  - `flopscope._opt_einsum.contract_path` still returns a flopscope
+    `PathInfo` with the same essential fields (`path`, `steps`,
+    `optimized_cost`, etc.).
+  - `StepInfo` no longer carries the dead symmetry-related fields
+    (`input_groups`, `output_group`, `inner_group`, `inner_applied`,
+    `dense_flop_cost`, `symmetry_savings`). It now has 4-5 fields:
+    `subscript`, `flop_count`, `input_shapes`, `output_shape` (plus
+    `merged_subset` if still used by display).
+
+- **Einsum cost model rewritten** to mirror the JS Symmetry-Aware Einsum
+  Contractions explorer's α/M direct-event model. The charged FLOP cost is
+  now path-independent: `(k - 1) · ∏ M_a + ∏ α_a` summed across components.
+  - `path_info.optimized_cost` returns the new whole-expression cost. For
+    expressions with declared symmetry, this number differs from the old
+    per-step `cost · unique/total` formula. See migration notes below.
+  - Path optimization no longer uses symmetry; `opt_einsum.contract_path`
+    behaves like upstream stock opt_einsum.
+  - Per-step `path_info.steps[i].flop_count` reverts to dense (no symmetry
+    adjustment per step).
+
+- **Reduction cost model rewritten** to the orbit-mapping model. Tier 1
+  (`np.ufunc.reduce` ops — sum, prod, max, min, all, any, bitwise_or/and/xor,
+  logical_or/and):
+  ```
+  cost = op_factor × (α - num_output_orbits) + extra_ops
+  ```
+  where α is the per-output-orbit input-orbit count summed across output
+  orbits. The `α - num_output_orbits` correction fixes #56's off-by-one
+  (dense `sum(n)` charges `n - 1`, was `n`). For symmetric inputs the new
+  model charges more than the legacy `unique_elements_for_shape` formula —
+  see #56 for the architectural shift.
+- **`np.median`, `np.percentile`, `np.quantile`** now use a Tier-2
+  output-discounted formula:
+  ```
+  cost = num_output_orbits × dense_per_output_cost
+  ```
+  For median/percentile/quantile, `dense_per_output_cost = axis_dim` (one
+  partition pass per output cell).
+- **`np.mean`** charges `sum_cost + num_output_orbits` (one divide per
+  output orbit; orbit-shared output values share the divisor).
+- `flopscope._flops.analytical_reduction_cost` body replaced with a
+  delegating call to `compute_reduction_accumulation_cost`. Signature
+  unchanged; numbers change.
+- `flopscope.accounting.reduction_cost` returns different numbers for both
+  dense and symmetric inputs (via the body change above).
+
+### Fixed
+
+- `flopscope.numpy.einsum_path` cache now keys on `fma_cost()` in addition
+  to `(subscripts, shapes, optimize)`. Previously, toggling
+  `flopscope.configure(fma_cost=2)` after a path was cached would return a
+  stale `PathInfo` whose per-step `flop_count` values were computed under
+  the old FMA convention.
+
+### Added
+
+- `flopscope.fma_cost()` — top-level re-export of the FMA-convention reader
+  (was `flopscope._cost_model.fma_cost`). Returns the current value of the
+  `fma_cost` setting (1 or 2).
+- `flopscope.einsum_clear_caches()` — clears the einsum path cache and the
+  einsum accumulation-cost cache together. Useful for cold-call benchmarks.
+  (The narrower `fnp.clear_einsum_cache()` continues to clear only the path
+  cache.)
+- `flopscope.einsum_cache_info()` — returns
+  `{"path": CacheInfo, "accumulation": CacheInfo}` so callers can inspect
+  both einsum caches in one call.
+- `flopscope.reduction_clear_cache()` and `flopscope.reduction_cache_info()`
+  — same pattern for the reduction accumulation-cost cache used by
+  `fnp.sum` / `fnp.mean` / `fnp.median` / etc.
+- `flopscope.clear_cache()` — convenience aggregate that clears both
+  einsum and reduction caches in one call.
+- `flopscope.tier2_reduction_cost(a, axis=None, *, dense_per_output_cost=None)`
+  — public inspection function for selection-style reductions
+  (`np.median` / `np.percentile` / `np.quantile`). Hides the
+  `op_factor=0, extra_ops=…` invocation pattern; `dense_per_output_cost`
+  defaults to the product of the reduced axes' lengths.
+- `flopscope.einsum_accumulation_cost(subscripts, *operands, partition_budget=None)`
+  — public inspection function returning the new `AccumulationCost` decomposition
+  (path-independent, per-component breakdown, regime trace).
+- `flopscope.AccumulationCost`, `flopscope.ComponentCost`, `flopscope.RegimeStep`
+  — public dataclasses.
+- New settings:
+  - `partition_budget` (default 100 000): per-component typed-partition cap.
+  - `dimino_budget` (default 500 000): whole-expression `G_pt` closure cap.
+- `CostFallbackWarning` now also fires when a partition counter exceeds its
+  budget; total falls back to `k · dense_baseline` (the no-symmetry direct-
+  event count).
+- New configurable setting `fma_cost` (default 1). Counts a fused
+  multiply-add as 1 op (hardware convention). Set to 2 to get the
+  textbook / opt_einsum convention.
+- `flopscope._cost_model.fma_cost()` function replaces the
+  `FMA_COST` constant. The constant is removed.
+- `flopscope.reduction_accumulation_cost(a, axis=None, *, op_factor=1, extra_ops=0)`
+  — public inspection function returning an `AccumulationCost` for a
+  reduction. Parallel to `einsum_accumulation_cost`.
+- Internal `_accumulation/_reduction.py`: `compute_reduction_accumulation_cost`
+  orchestrator, `output_discounted_reduction_cost` (Tier 2), and
+  `_normalize_axis` / `_num_output_orbits` helpers.
+- `_accumulation/_cache.py`: `_reduction_cache` + `get_reduction_cost_cached`
+  (LRU 4,096).
+- `_accumulation/_cost.py:aggregate_reduction` body implemented (was a
+  signature-locked `NotImplementedError`).
+
+### Removed
+
+- `flopscope._opt_einsum._paths` — now upstream
+- `flopscope._opt_einsum._path_random` — now upstream
+- `flopscope._opt_einsum._parser` — now upstream (re-exported via shim)
+- `flopscope._opt_einsum._blas` — was unused dead code
+- `flopscope._opt_einsum._testing` — was unused dead code
+- `flopscope._opt_einsum._typing` — now upstream
+- `flopscope._cost_model.FMA_COST` constant — replaced by `fma_cost()`
+- `flopscope._opt_einsum._subgraph_symmetry` — internal module deleted.
+- `flopscope._opt_einsum._symmetry` — internal module deleted (was mostly
+  `symmetric_flop_count`, `unique_elements`, `SubsetSymmetry` — all only used
+  by the deleted oracle).
+- `use_inner_symmetry` setting — was a knob on the deleted oracle.
+
+---
+
 ### BREAKING
 
 - `BudgetContext.untracked_time` and `summary_dict()["untracked_time_s"]` now

Makefile

@@ -103,6 +103,13 @@ install:  ## Install all deps (dev + docs) and set up git hooks
 	uv sync --all-extras
 	git config core.hooksPath .githooks
 
+# ---------------------------------------------------------------------------
+# Benchmarks
+# ---------------------------------------------------------------------------
+.PHONY: bench-accumulation
+bench-accumulation:  ## Cold + warm latency benchmark for einsum_accumulation_cost
+	$(UV) python benchmarks/accumulation/bench_cost_compute.py
+
 # ---------------------------------------------------------------------------
 # Help
 # ---------------------------------------------------------------------------

README.md

@@ -73,7 +73,7 @@ for i, W in enumerate(weights):
     h = fnp.einsum('ij,j->i', W, h)
     if i < depth - 1:
         h = fnp.maximum(h, 0)
-flops.budget_summary()  # 984,321 FLOPs
+flops.budget_summary()  # 6,231,041 FLOPs
 ```
 
 </td>
@@ -87,6 +87,7 @@ flops.budget_summary()  # 984,321 FLOPs
 - **Budget enforcement** -- operations are checked before execution; exceeding the budget raises a clear error
 - **Symmetry-aware einsum** -- automatic FLOP savings for repeated operands and declared symmetry groups
 - **Transparent diagnostics** -- inspect per-operation costs, cumulative budget usage, and detailed summaries at any time
+- **Inspect costs analytically** -- `flops.einsum_accumulation_cost(...)` and `flops.reduction_accumulation_cost(...)` return the exact FLOP count for an einsum or reduction without running the op
 - **Truncated SVD** -- top-k singular value decomposition with `O(m * n * k)` cost
 
 ## What's Supported
@@ -143,16 +144,16 @@ flops.budget_summary()           # accumulated session/global summary
 ```
 flopscope FLOP Budget Summary
 =============================
-  Total budget:     100,000,000
-  Used:                 984,321  (1.0%)
-  Remaining:         99,015,679  (99.0%)
+  Total budget:             100,000,000
+  Used:                       6,231,041  (6.2%)
+  Remaining:                 93,768,959  (93.8%)
 
   By operation:
-    random.randn          327,936  ( 33.3%)  [6 calls]
-    multiply              327,680  ( 33.3%)  [5 calls]
-    einsum                327,680  ( 33.3%)  [5 calls]
-    maximum                 1,024  (  0.1%)  [4 calls]
-    sqrt                        1  (  0.0%)  [1 call]
+    random.randn            5,246,976  ( 84.2%)  [6 calls]
+    einsum                    655,360  ( 10.5%)  [5 calls]
+    multiply                  327,680  (  5.3%)  [5 calls]
+    maximum                     1,024  (  0.0%)  [4 calls]
+    sqrt                            1  (  0.0%)  [1 call]
 
   Total Wall Time:     ...s
   Flopscope Backend:   ...s  (...%)
@@ -164,11 +165,26 @@ flopscope FLOP Budget Summary
 
 ```python
 # Query FLOP costs without running anything (no BudgetContext needed)
-cost = flops.einsum_cost('ij,jk->ik', shapes=[(256, 256), (256, 256)])
-print(f"Matmul cost: {cost:,}")  # 16,777,216
+cost = flops.accounting.einsum_cost('ij,jk->ik', shapes=[(256, 256), (256, 256)])
+print(f"Matmul cost: {cost:,}")  # 33,554,432
 
-cost = flops.svd_cost(m=256, n=256, k=10)
-print(f"SVD cost: {cost:,}")     # 655,360
+cost = flops.accounting.svd_cost(m=256, n=256, k=10)
+print(f"SVD cost: {cost:,}")     # 2,621,440
+```
+
+For symmetry-aware inspection that takes actual array inputs (and reflects
+declared symmetry), use the accumulation cost APIs:
+
+```python
+import numpy as np
+A = np.zeros((256, 256))
+B = np.zeros((256, 256))
+# Returns an AccumulationCost decomposition without running the op
+cost = flops.einsum_accumulation_cost('ij,jk->ik', A, B)
+print(cost.total)                # 33,554,432
+
+cost = flops.reduction_accumulation_cost(A, op_factor=1)
+print(cost.total)                # 65,535
 ```
 
 ### Symmetry Savings

benchmarks/_contractions.py

@@ -27,7 +27,7 @@
     "outer": "M*N",
     "tensordot": "product of free * contracted dims",
     "kron": "numel(output)",
-    "einsum": "product of index dims (FMA=1)",
+    "einsum": "α/M model (FMA=2 textbook)",
 }
 
 _BENCHMARK_SIZE_STRINGS: dict[str, str] = {
@@ -59,27 +59,28 @@ def _analytical_cost(op: str, **kwargs: int) -> int:
         Analytical FLOP count.
     """
     costs: dict[str, int] = {
-        # dot: 2D matrix multiply A(512,512) @ B(512,512), FMA=1 op
+        # dot: 2D matrix multiply A(512,512) @ B(512,512), FMA=2 textbook
         "dot": 512 * 512 * 512,
         # matmul: identical to dot for 2D
         "matmul": 512 * 512 * 512,
         # inner: dot product of two 1M-element vectors.
-        # Runtime charges a.size — matches flopscope's convention (FMA=1).
+        # Runtime charges a.size — matches flopscope's convention (FMA=2, but
+        # a.size is pointwise-shaped so the FMA off-by-one doesn't apply here).
         "inner": 1_000_000,
         # vdot: same as inner for 1D real inputs.
-        # Runtime charges a.size (FMA=1).
+        # Runtime charges a.size (FMA=2, pointwise-shaped — no off-by-one).
         "vdot": 1_000_000,
         # vecdot: batched dot product A(1000,512) . B(1000,512)
         # Output (1000,) with contracted axis 512.
-        # Runtime charges result.size * contracted = 1000 * 512 (FMA=1).
+        # Runtime charges result.size * contracted = 1000 * 512 (FMA=2 textbook).
         "vecdot": 1000 * 512,
         # outer: outer product of two 5000-element vectors
         "outer": 5000 * 5000,
         # tensordot: A(64,64,64) . B(64,64,64) axes=1 -> contract last of A with first of B
         "tensordot": 64**5,
         # kron: Kronecker product A(64,64) x B(64,64)
         "kron": 64**4,
-        # einsum: 'ij,jk->ik' is matrix multiply (512,512)x(512,512), FMA=1
+        # einsum: 'ij,jk->ik' is matrix multiply (512,512)x(512,512), FMA=2 textbook
         "einsum": 512 * 512 * 512,
     }
     return costs[op]

benchmarks/_linalg_delegates.py

@@ -79,19 +79,24 @@ def _analytical_cost(op_name: str) -> int:
     costs: dict[str, int] = {
         "cond": 512 * 512 * 512,  # m*n*min(m,n) via SVD
         "cross": 6 * 1_000_000,  # 6*n
-        "matmul": 512 * 512 * 512,  # M*N*K (FMA=1)
+        "matmul": 2 * 512 * 512 * 512
+        - 512 * 512,  # 2*M*N*K - M*N (FMA=2); = 268,173,312
         "matrix_norm": 2 * 512 * 512,  # 2*numel (Frobenius)
         "matrix_power": 3 * 64**3,  # 3 matmuls for n=5
         "matrix_rank": 512 * 512 * 512,  # m*n*min(m,n) via SVD
-        "multi_dot": 128 * 64 * 128 + 128 * 128 * 64,  # FMA=1
-        "norm": 10_000_000,  # numel (L2)
+        "multi_dot": 128 * 64 * 128
+        + 128
+        * 128
+        * 64,  # optimal chain (FMA=2); coincidentally same as FMA=1 value = 2,097,152
+        "norm": 2 * 10_000_000,  # 2*numel (FMA=2, vector L2)
         "outer": 5000 * 5000,  # M*N
-        "tensordot": 64**5,  # d^5 (FMA=1)
+        "tensordot": 64
+        ** 5,  # d^5 (FMA=2 textbook; matches flopscope charge = 1,073,741,824)
         "tensorinv": 64**3,  # n^3 after reshape
         "tensorsolve": 64**3,  # n^3 after reshape
         "trace": 10_000,  # min(m,n)
         "vecdot": 1000 * 512,  # batch*K
-        "vector_norm": 10_000_000,  # numel (L2)
+        "vector_norm": 2 * 10_000_000,  # 2*numel (FMA=2)
     }
     return costs[short]
 

benchmarks/_polynomial.py

@@ -20,7 +20,7 @@
 ]
 
 _FORMULA_STRINGS: dict[str, str] = {
-    "polyval": "n * degree (FMA=1)",
+    "polyval": "2 * n * degree (FMA=2)",
     "polyfit": "2 * n * (degree+1)^2",
     "roots": "degree^3",
     "polymul": "(degree+1)^2",
@@ -40,7 +40,9 @@ def _analytical_cost(op: str, n: int, degree: int) -> int:
     benchmark denominator and the budget deduction use the same formula.
     """
     if op == "polyval":
-        return n * degree
+        return (
+            2 * n * degree
+        )  # Updated for FMA=2 unification (spec 2026-05-20): polyval formula doubled m*deg → 2*m*deg.
     elif op == "polyfit":
         return 2 * n * (degree + 1) ** 2
     elif op == "roots":