diff --git a/Makefile b/Makefile index ce9056ca2..a854eba53 100644 --- a/Makefile +++ b/Makefile @@ -289,10 +289,11 @@ cli-demo: cli $$PRED from QUBO --hops 1; \ \ echo ""; \ - echo "--- 5. path: find reduction paths ---"; \ + echo "--- 5. path: asymptotic Pareto front (no --size) ---"; \ $$PRED path MIS QUBO; \ - $$PRED path MIS QUBO -o $(CLI_DEMO_DIR)/path_mis_qubo.json; \ $$PRED path Factoring SpinGlass; \ + echo "--- 5b. path --cost: single concrete path (for reduce --via) ---"; \ + $$PRED path MIS QUBO --cost minimize-steps -o $(CLI_DEMO_DIR)/path_mis_qubo.json; \ $$PRED path MIS QUBO --cost minimize:num_variables; \ \ echo ""; \ diff --git a/docs/design/symbolic-growth-domain.md b/docs/design/symbolic-growth-domain.md new file mode 100644 index 000000000..e5ae1fdf8 --- /dev/null +++ b/docs/design/symbolic-growth-domain.md @@ -0,0 +1,359 @@ +# Symbolic Growth Domain & Pareto Path Search — Product Design + +Status: approved design, ready for decomposition into issues. +Origin: issue #1069 (`pred path --all` OOMs/hangs in `big_o_normal_form`). The acute +symptom is already mitigated on `main` by a stopgap: `MAX_CANONICAL_TERMS = 50_000` +in `canonical.rs` aborts oversized expansions, and the CLI falls back to printing the +*unreduced* composed expression as `O()` on failure +(`problemreductions-cli/src/commands/graph.rs:349`). This design replaces +refuse-or-bluff with a system that answers. + +## Need + +The symbolic overhead system conflates exact expressions with asymptotic queries: +`big_o_normal_form` (src/big_o.rs) fully expands composed path overheads to monomial +normal form (src/canonical.rs) before projecting to Big-O. Expansion of nested +`(sum)^2 * (sum)^2` structures is exponential in nesting depth — the root cause of +issue #1069. The stopgap cap prevents the OOM but leaves three structural defects: + +1. **Refuse-or-bluff answers.** Paths whose composed overhead exceeds the expansion + cap get no normalized Big-O; the CLI falls back to printing the raw unreduced + expression disguised as `O(...)`. The exponential-expansion algorithm is still + there, merely fenced. +2. **Heuristic dominance.** Asymptotic comparison relies on a foolable two-point + numerical sampling heuristic (`numerical_dominance_check`) — e.g. `n^100` vs + `1.001^n` is decided wrongly because the crossover lies beyond the sampled range. +3. **Unsound search.** The scalar Dijkstra in `ReductionGraph::find_cheapest_path` + has a latent correctness hole: edge costs depend on the size accumulated along the + path, which violates Dijkstra's assumptions — a cheaper-so-far path with a larger + intermediate size can be wrongly preferred. And there is no instance-free + (asymptotic) search mode at all. + +We need a **trustworthy** (explicit semantic axioms, bounded termination, per-rule +verifiability) and **extensible** (new functions/variables without touching the core) +symbolic system: an exact `Expr` layer separated from an asymptotic growth domain, +with both Big-O rendering and path search running in the asymptotic domain at +polynomial cost. Occam's razor is a hard constraint: no new entities beyond what the +selected features require. + +**Users:** library maintainers adding models/rules; CLI/MCP consumers of +`pred path` / `find_path`; the Typst paper's auto-derivation pipeline. + +**Success criteria** (the stopgap already prevents OOM; these measure what the +principled system adds): +- **Answers, not refusals:** every enumerable path gets a genuine normalized Big-O. + The `MAX_CANONICAL_TERMS` bail-out and the `O()` CLI fallback are + deleted; the only remaining "cannot normalize" sources are nonlinear exponents + and factorials, rendered as an explicit annotation (the one `2^num_vertices` + overhead edge gets a real exponential bound via the linear `exp` field). + Regression: issue #1069's exploding path (KSat → … → QuadraticAssignment → ILP → + QUBO) asserts a real normalized Big-O, not an error or fallback. +- **Trustworthy comparison:** the numerical sampling heuristic is replaced by a + symbolic decision procedure, property-tested against numeric evaluation. +- **Correct search:** Pareto label search fixes the path-dependent-cost hole and adds + an instance-free asymptotic mode. +- Big-O for all enumerated paths across the whole reduction graph completes within a + CI time budget (each test < 5 s per repo policy). +- Output is byte-identical across Linux/macOS (no inventory-order dependence). + +**Constraints:** +- The `#[reduction]` macro and overhead declaration syntax stay unchanged (dozens of + rule files untouched). +- Internal APIs and CLI output format may break (0.x semver). +- No new external dependencies. + +## Prior art & landscape + +Surveyed via four research passes (CAS systems; compiler symbolic-cost systems; +e-graph engines; asymptotics theory and formalization). Borrow-vs-build verdict: + +| Candidate | Verdict | Why | +|---|---|---| +| Albert–Alonso–Arenas–Genaim–Puebla, *Asymptotic Resource Usage Bounds* (APLAS 2009) | **Adopt as spec** | Published normal form (sums of products of `2^(r·A)`, `A^r`, `log A`) with a soundness theorem `e ∈ Θ(asymp(e))` — our correctness contract | +| SageMath `AsymptoticRing` / growth groups | **Borrow the design, not the code** | GPL; the core (exponent-vector arithmetic + poset of summands with O-term absorption) is small enough to reimplement cleanly | +| KoAT weakly-monotone bound grammar (Brockschmidt et al., TOPLAS 2016) | **Adopt as axiom** | Weak monotonicity ⇒ composition-by-substitution is sound ⇒ Pareto label search is correct (isotonicity) | +| LLVM SCEV / GCC chrec | **Adopt patterns** | Construction-time canonicalization, explicit budgets with graceful degradation, absorbing "don't know" sentinel (`SCEVCouldNotCompute`, `chrec_dont_know`) | +| Multivariate Big-O semantics: Howell (KSU TR 2007-4); Guéneau–Charguéraud–Pottier (ESOP 2018) | **Adopt definition** | Naive multivariate O is inconsistent (Howell Thm 2.3/2.4); the product-filter definition restricted to nonnegative weakly-monotone functions is the trustworthy one | +| McRAPTOR / OpenTripPlanner `ParetoSet` / nigiri `pareto_set.h`; Martins 1984; NAMOA* | **Adopt algorithm** | Per-node label bags (antichains) with dominance pruning are the industry and literature standard for partial-order path costs; enumerate-then-filter appears nowhere as a recommended method | +| ProblemReductions.jl `reduction_paths` | **Anti-pattern baseline** | `all_simple_paths` with no cost model, no ranking, no filter; survives only because its graph is tiny | +| egg / egglog e-graphs | **Dropped** | Directional normalization doesn't need equality saturation (Cranelift aegraph retrospective: mean e-class size 1.13); egglog API unstable | +| SymPy / GiNaC / Symbolica | **Concepts only** | Never auto-expand; deterministic total order on atoms; function-registry extensibility (deferred with F6) | + +Nothing is directly reusable as a dependency; this is a build against published specs. + +**Empirical inventory scan** (drives the grammar decision): registered overhead +expressions are overwhelmingly polynomial with subtraction and constant division. +Exceptions: one `log` factor (`ksatisfiability_*`: `(num_vars + num_clauses)^2 * +log(num_vars + num_clauses + 1)`), one genuine exponential +(`highlyconnecteddeletion_ilp.rs`: `num_vars = "2^num_vertices"`), and one +`sqrt((x)^2)` used as an absolute-value idiom. `declare_variants!` complexity strings +are heavily exponential, but they are consumed only by `pred list/show` display and +the dropped F8 — outside this design's data path. + +## Features + +Selected (rough, agentic-coding-adjusted estimates): + +| # | Feature | Effort | +|---|---|---| +| F1 | Growth domain: `GrowthTerm`/`Growth` antichain, symbolic dominance, pruning, absorbing `Unknown`, caps with upward widening | ~2–3 days | +| F2 | Replace the `big_o.rs` pipeline with the growth domain; delete `canonical.rs`; issue-1069 regression + whole-graph CI budget tests | ~1–2 days | +| F3 | Pareto label search kernel replacing `dijkstra`, with two label domains: F3a asymptotic (`Growth` per size field) and F3b concrete instance (**measured**: execute reductions and apply post-construction measured budgets) | ~3–4 days | +| F12 | Per-edge overhead calibration test: canonical examples run through `reduce_to()`, measured sizes must not exceed formula predictions | ~0.5–1 day | +| F4 | CLI/MCP surface: Pareto-front output, deterministic ordering, `--json` no longer renders text | ~1–2 days | +| F5+F11 (merged support work, folded into F1/F3/F4) | Redundancy check (`find_dominated_rules`) rewired to the same dominance order; `Growth` serde + `Display` consumed by CLI JSON and paper export | ~1.5 days | + +Total: ~10–14 days. + +Deferred / dropped, with reasons: + +- **F6 `Expr::Func(FuncKind)` registry** and **F7 shared parser crate** — deferred to a + later milestone. Genuine extensibility improvements, but independent of this + milestone's goal; the growth domain consumes `Expr` as-is. +- **F8 effective-complexity ranking** (target complexity ∘ overhead) — deferred until a + concrete find-problem need; requires an exponential part in `GrowthTerm` (see + Extensibility). +- **F9 convex-hull/AM-GM pruning** — deferred until antichain sizes measurably hurt; + Pareto pruning suffices at current variable counts. +- **F10 egg-based display simplification** — dropped per survey (directional ruleset + does not need equality saturation). + +## Semantic foundation (normative) + +These definitions and axioms are the trust contract; tests enforce them. + +- **Definition (multivariate Big-O, product filter).** For size functions + `f, g : ℕ_{≥2}^k → ℝ_{≥0}`: `g ∈ O(f)` iff `∃ c > 0, N` such that + `g(x) ≤ c·f(x)` whenever **all** variables `x_i ≥ N`. (Howell's `O_∀`; + Guéneau et al.'s product filter.) +- **Domain axioms.** Every expression admitted to the growth domain is nonnegative + and weakly monotone (nondecreasing in each variable) on `vars ≥ 2`. Under these + axioms Howell's inconsistencies vanish and `f + g ≍ max(f, g)` up to a constant + factor, which licenses `add = antichain union + prune`. +- **Widening rules (always upward, i.e. toward a valid upper bound):** + - Subtraction: `a − b ⇝ a + b` (sound since `b ≥ 0`; also covers the + `sqrt((a−b)^2)` absolute-value idiom because `|a−b| ≤ a+b`). + - Constant division and all multiplicative constants: dropped on entry. + - Exponentials with **linear** exponents (`c^x`, `c^(r·x)`, `exp(x)`) are + first-class (see M1's `exp` field). Nonlinear exponents (`2^(n*k)`, + `2^sqrt(n)`, double exponentials), `factorial(·)`, and negative exponents: + `Growth::Unknown` (absorbing). +- **Forbidden moves (documented + tested):** never specialize a variable to a + constant inside an O-fact; never rescale coefficients of exponents + (`2^(2n) ∉ O(2^n)` — exp rates compare coefficientwise, exactly). +- **Isotonicity invariant (for search):** if label `A` dominates label `B`, then for + any edge `e`, `extend(A, e)` dominates `extend(B, e)`. This follows from the + monotonicity axiom (composition by substitution into monotone expressions) and is + the correctness condition for dominance pruning in M3. + +## Modules + +Only one new file. Everything else is in-place replacement; net LOC is expected +near zero or negative (`canonical.rs`, 431 lines, is deleted). + +### M1 — `src/growth.rs` (the one new entity) + +```rust +/// One growth monomial, e.g. 2^(3k)·n^2·m·log(n) → +/// { exp: {k:3.0}, poly: {n:2.0, m:1.0}, logs: {n:1} }. +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +pub struct GrowthTerm { + exp: BTreeMap<&'static str, f64>, // variable → rate, base normalized to 2 + // (3^n → {n: log2(3)}); linear forms only + poly: BTreeMap<&'static str, f64>, // variable → degree (0.5 covers sqrt) + logs: BTreeMap<&'static str, u32>, // variable → log power +} + +#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] +pub enum Growth { + /// Antichain of pairwise-incomparable dominant terms, sorted by a + /// deterministic total order (for stable output/serialization). + Terms(Vec), + /// Absorbing sentinel: exp/factorial/negative exponents, or cap overflow + /// that even widening cannot represent. Absorbs through all operations. + Unknown, +} +``` + +Operations (each prunes back to an antichain immediately): + +- `Growth::from_expr(&Expr) -> Growth` — single bottom-up pass, linear in tree size. + `Var → {poly:{v:1}}`; `Const → O(1)` (empty term); `Add → union + prune`; + `Mul → pairwise map-merge + prune`; `Pow(base, const k ≥ 0) →` compute base's + antichain, then pairwise products (never expands the underlying sums); + `Log(a) → log(dominant(a))` using `log(n^a·m^b) ≍ log n + log m`; + `Sqrt = Pow 0.5`; everything else → `Unknown`. +- `dominates(&GrowthTerm, &GrowthTerm) -> bool` — per variable, lexicographic on + (exp rate, poly degree, log power); dominated iff ≤ on every variable and < on at + least one. This decides e.g. `1.001^n ≻ n^100` correctly, which the sampling + heuristic gets wrong. + Purely symbolic; replaces `numerical_dominance_check`. +- Caps: antichain length cap (default 32). On overflow, **widen upward** to the + single term taking the componentwise max of all exponents (a valid upper bound), + never truncate by order. +- Axiom guards: `debug_assert!` nonnegativity/monotonicity preconditions at entry. + +Deps: read-only on `expr.rs`. Serde derive here is the whole of former F11. + +### M2 — `big_o.rs` pipeline replacement + +`big_o_normal_form(&Expr) -> Result` keeps its +signature: internally `Growth::from_expr` → render `Growth` back to a display `Expr` +(`Unknown` maps to the existing `Unsupported` error). CLI callers (`big_o_of`, +`overhead_to_json`, `format_path_text`) are untouched. `compose_path_overhead` +continues to produce the compact nested `Expr` (≤ ~2 KB in the worst observed case); +`from_expr` walks it in microseconds — **no caching, no registry changes**. +`canonical.rs` and the `asymptotic_normal_form` compatibility wrapper are deleted +along with their unit tests (internal API breakage is in-scope). + +`pred-sym` (the standalone symbolic CLI, used by the find-problem skills for +`big-o` and `eval`) follows suit: the `canon` subcommand is deleted (no live +consumers), and `compare` narrows its semantics to Big-O equivalence via the growth +domain. `big-o` keeps working on the skills' effective-complexity inputs +(`1.5^n * n^2`) thanks to the linear `exp` field; nonlinear-exponent inputs report +`Unknown` and the skills fall back to `pred-sym eval`. + +Alternatives considered: capped expansion (rejected: keeps the exponential algorithm +and reintroduces order-dependent truncation); per-edge growth caching in +`ReductionEntry` with per-path folding (rejected for now: YAGNI at current graph +size; revisit if profiling ever shows `from_expr` on composed paths as hot). + +### M3 — Pareto label search kernel (`src/rules/graph.rs`, in-place) + +Replace `dijkstra` (~60 lines) with one generic label-setting search (~100 lines) +plus a minimal trait: + +```rust +pub trait PathLabel: Clone { + fn extend(&self, edge: &ReductionEdge) -> Self; // must be isotone + fn dominates(&self, other: &Self) -> bool; // partial order +} +``` + +- Per-node **bag** = antichain of non-dominated labels, each with a predecessor + pointer for path reconstruction (McRAPTOR structure). +- Deterministic bounding, in the style of transit routers: hop cap (default 16) and + per-node bag cap with a **deterministic tie-break** (fewest hops, then + lexicographic node-name order) — never iteration-order truncation. A label evicted + from a bag (dominated or cap-truncated) has its arena slot's label freed immediately, + so the bag cap genuinely bounds retained per-node label memory. +- Label domains: + - **F3a asymptotic:** label = `BTreeMap` mapping each size field of + the current node to its growth in the source's variables; `extend` substitutes + the edge's overhead expressions; `dominates` is componentwise. Exponential + growth is comparable via the `exp` field (polynomial paths dominate exponential + ones); `Unknown` fields make a label dominated by any known label — undecidable + paths rank last, which is the honest ranking. + - **F3b instance (measured):** for a concrete instance, formulas are advisory — + **measured sizes are authoritative**. Overhead formulas are scaling upper bounds + over the declared size fields and can be arbitrarily loose on + structure-dependent constructions (see #107), so they must never arbitrate + between concrete candidates. Label = the actual `ProblemSize` measured on the + constructed intermediate problem (plus the reduction chain itself, reused for + solving/witness extraction by the winner); `extend` executes the edge's + `reduce_to()` and measures. The only instance-budget guard is the **measured + budget check after execution**. Evaluating an asymptotic expression at one point + is not a certified concrete bound, so overhead formulas do not prune measured + candidates. This also means the budget cannot prevent the construction itself + from exhausting memory. + + Measured search uses **no dominance pruning**. `ProblemSize` omits instance + structure, and equal-size intermediate instances can produce different sizes under + a later structure-dependent reduction. Even serialized-state equivalence is not + used to discard a route. It is therefore a separate exhaustive simple-path + enumeration, not a label domain in the capped Pareto kernel. + + Note the measured label deliberately does **not** use branch-and-bound: a + reduction can *shrink* the measured size, so the cost is non-monotone and a + B&B bound could prune a partial route that would still finish smallest. + No hop or bag cap truncates this enumeration, so its time and retained constructed + state can grow exponentially with the number of simple paths. This also does not + bound temporary memory used inside `reduce_to()`. + This fixes the path-dependent-cost hole in the current Dijkstra *and* removes + the dependency on formula accuracy for concrete decisions. +- `find_cheapest_path*` become thin wrappers returning the front (instance mode + typically collapses to a single optimum after the numeric tie-break). +- `find_dominated_rules` / `compare_overhead` (`src/rules/analysis.rs`) are rewired + to the same `dominates` order, deleting their bespoke comparison heuristics — + one trusted comparison everywhere (former F5). +- `all_simple_paths`-based enumeration remains the explicit `--all` listing mechanism; + measured optimum-finding now performs its own execution-aware simple-path enumeration + because no sound state-level dominance relation is available. + +Alternatives considered: enumerate-then-filter (rejected: combinatorial growth as the +graph densifies, and any truncation limit is iteration-order-dependent — the sibling +package ProblemReductions.jl does exactly this, with no cost model, and it is the +baseline we are improving on); a generic semiring algebraic-path framework (rejected: +over-engineering for two label domains); formula-evaluated instance labels (rejected +after review: overhead formulas are upper bounds over declared size fields and can be +arbitrarily loose on structure-dependent constructions, so a formula-ranked front may +not contain the true winner — measured sizes are the ground truth and affordable at +interactive scales; formulas remain available for asymptotic analysis but do not +decide concrete feasibility). + +### M4 — CLI/MCP surface (`problemreductions-cli/src/commands/graph.rs`, in-place) + +- Asymptotic `pred path S T`: print the Pareto front (typically 1–3 paths), each with + its Big-O per size field; paths whose composed growth is `Unknown` (nonlinear + exponents, factorial) are annotated explicitly instead of showing a fake bound. +- Instance mode (`--size …`): output shape unchanged (single best path). +- `path --all`: keep enumeration; Big-O per path now via M2 (fast); **`--json` mode + no longer builds the text rendering** (the unconditional `format_path_text` call + named in issue #1069). +- All path lists sorted by (hops, lexicographic names). JSON emits the structured + `Growth` serialization. (The paper export consumes raw overhead expressions, not + Big-O strings — verified unaffected.) + +## Quality requirements + +- **Reliability:** every public function terminates with an answer or `Unknown` — + no input can hang or OOM. Regression: issue #1069 path #34; a whole-graph test + enumerating paths (bounded length) between hot pairs asserts Big-O completion + within the CI budget (< 5 s per test). +- **Trustworthiness testing:** each `from_expr` transfer function and the dominance + order get randomized property tests (≥ 5000 checks, matching the repo's + verify-reduction culture): `eval(expr) ≤ C · eval(render(growth(expr)))` at large + sizes; `growth` idempotent on its own rendering; `dominates(a,b)` ⟹ sampled + `eval(b)/eval(a)` grows. Isotonicity of both `PathLabel` impls is property-tested. +- **Determinism:** identical output across platforms; a test compares `pred path` + output against golden files (antichain and front ordering are total and + deterministic by construction). +- **Performance:** `pred path KSat QUBO --all` end-to-end < 1 s (currently OOM). +- **Extensibility:** the linear `exp` field ships in M1 (required by the + find-problem skills' use of `pred-sym big-o` on effective-complexity + expressions). The remaining upgrade path — nonlinear exponents (a polynomial + exponent instead of a linear form), needed only if F8-style effective-complexity + ranking over complexity strings like `2^(num_edges * k)` is ever built — touches + only `dominates`, `mul`, and `from_expr`'s `Pow/Exp` arms; antichain machinery, + caps, search kernel, and serialization are unaffected. + +## Out of scope + +- `#[reduction]` macro, overhead declaration syntax, and all rule files. +- `declare_variants!` complexity strings and their validation + (`is_valid_complexity_notation`) — untouched; they are display-only in this design. +- FuncKind registry, shared parser crate, effective-complexity ranking, hull pruning, + egg display layer (deferred/dropped as listed under Features). + +## References + +- E. Albert, D. Alonso, P. Arenas, S. Genaim, G. Puebla. *Asymptotic Resource Usage + Bounds.* APLAS 2009. (Normal form + `Θ`-preservation theorem.) +- R. Howell. *On Asymptotic Notation with Multiple Variables.* Kansas State + University TR 2007-4. (Multivariate O inconsistencies; `O_∀` definition.) +- A. Guéneau, A. Charguéraud, F. Pottier. *A Fistful of Dollars: Formalizing + Asymptotic Complexity Claims via Deductive Program Verification.* ESOP 2018. + (Filter-based O; nonnegative-monotone cost discipline; documented pitfalls.) +- M. Brockschmidt, F. Emmes, S. Falke, C. Fuhs, J. Giesl. *Analyzing Runtime and Size + Complexity of Integer Programs.* TOPLAS 2016. (Weakly monotone bounds compose.) +- SageMath `sage.rings.asymptotic` (growth groups, O-term absorption) — design + reference only (GPL). +- LLVM `ScalarEvolution` / GCC `tree-chrec` — budgets, sentinels, construction-time + canonicalization. +- D. Delling, T. Pajor, R. Werneck. *Round-Based Public Transit Routing.* ALENEX + 2012 (McRAPTOR bags); E. Martins. *On a Multicriteria Shortest Path Problem.* EJOR + 1984; L. Mandow, J.-L. Pérez de la Cruz. *Multiobjective A\* with Consistent + Heuristics.* JACM 2010. +- D. Gruntz. *On Computing Limits in a Symbolic Manipulation System.* ETH 1996 + (dominance ordering; relevant when the `exp` field is added). +- Issue #1069 — root-cause analysis this design responds to. diff --git a/docs/src/cli.md b/docs/src/cli.md index e94f4456e..049bb5de3 100644 --- a/docs/src/cli.md +++ b/docs/src/cli.md @@ -163,7 +163,7 @@ Show all paths or save for later use with `pred reduce --via`: ```bash pred path MIS QUBO --all # all paths (up to 20) pred path MIS QUBO --all --max-paths 50 # increase limit -pred path MIS QUBO -o path.json # save path for `pred reduce --via` +pred path MIS QUBO -o path.json # save front + best path for `pred reduce --via` pred path MIS QUBO --all -o paths/ # save all paths to a folder ``` diff --git a/problemreductions-cli/src/bin/pred_sym.rs b/problemreductions-cli/src/bin/pred_sym.rs index daa28bf7a..c20c2b97b 100644 --- a/problemreductions-cli/src/bin/pred_sym.rs +++ b/problemreductions-cli/src/bin/pred_sym.rs @@ -1,5 +1,5 @@ use clap::{Parser, Subcommand}; -use problemreductions::{big_o_normal_form, canonical_form, Expr, ProblemSize}; +use problemreductions::{big_o_normal_form, Expr, ProblemSize}; #[derive(Parser)] #[command( @@ -19,11 +19,6 @@ enum Commands { /// Expression string expr: String, }, - /// Compute exact canonical form - Canon { - /// Expression string - expr: String, - }, /// Compute Big-O normal form BigO { /// Expression string @@ -33,7 +28,7 @@ enum Commands { #[arg(long)] raw: bool, }, - /// Compare two expressions (exits with code 1 if neither exact nor Big-O equal) + /// Compare two expressions for Big-O equivalence (exits 1 if not equal) Compare { /// First expression a: String, @@ -69,16 +64,6 @@ fn main() { let parsed = parse_expr_or_exit(&expr); println!("{parsed}"); } - Commands::Canon { expr } => { - let parsed = parse_expr_or_exit(&expr); - match canonical_form(&parsed) { - Ok(result) => println!("{result}"), - Err(e) => { - eprintln!("Error: {e}"); - std::process::exit(1); - } - } - } Commands::BigO { expr, raw } => { let parsed = parse_expr_or_exit(&expr); match big_o_normal_form(&parsed) { @@ -98,29 +83,31 @@ fn main() { Commands::Compare { a, b } => { let expr_a = parse_expr_or_exit(&a); let expr_b = parse_expr_or_exit(&b); - let canon_a = canonical_form(&expr_a); - let canon_b = canonical_form(&expr_b); let big_o_a = big_o_normal_form(&expr_a); let big_o_b = big_o_normal_form(&expr_b); println!("Expression A: {a}"); println!("Expression B: {b}"); - let mut exact_equal = false; - let mut big_o_equal = false; - if let (Ok(ca), Ok(cb)) = (&canon_a, &canon_b) { - exact_equal = ca == cb; - println!("Canonical A: {ca}"); - println!("Canonical B: {cb}"); - println!("Exact equal: {exact_equal}"); - } - if let (Ok(ba), Ok(bb)) = (&big_o_a, &big_o_b) { - big_o_equal = ba == bb; - println!("Big-O A: O({ba})"); - println!("Big-O B: O({bb})"); - println!("Big-O equal: {big_o_equal}"); - } - if !exact_equal && !big_o_equal { - std::process::exit(1); + match (&big_o_a, &big_o_b) { + (Ok(ba), Ok(bb)) => { + // Rendering is canonical, so equal growth ⇒ equal Big-O expr. + let big_o_equal = ba == bb; + println!("Big-O A: O({ba})"); + println!("Big-O B: O({bb})"); + println!("Big-O equal: {big_o_equal}"); + if !big_o_equal { + std::process::exit(1); + } + } + _ => { + if let Err(e) = &big_o_a { + println!("Big-O A: "); + } + if let Err(e) = &big_o_b { + println!("Big-O B: "); + } + std::process::exit(1); + } } } Commands::Eval { expr, vars } => { diff --git a/problemreductions-cli/src/cli.rs b/problemreductions-cli/src/cli.rs index 70fa1e5af..719e49be5 100644 --- a/problemreductions-cli/src/cli.rs +++ b/problemreductions-cli/src/cli.rs @@ -112,11 +112,11 @@ Use `pred to ` for incoming neighbors (what reduces to this).")] /// Find the cheapest reduction path between two problems #[command(after_help = "\ Examples: - pred path MIS QUBO # cheapest path + pred path MIS QUBO # asymptotic Pareto front (Big-O per size field) pred path MIS QUBO --all # all paths - pred path MIS QUBO -o path.json # save for `pred reduce --via` + pred path MIS QUBO -o path.json # save front + best path for `pred reduce --via` pred path MIS QUBO --all -o paths/ # save all paths to a folder - pred path MIS QUBO --cost minimize:num_variables + pred path MIS QUBO --cost minimize:num_variables # single cheapest path by a scalar cost (also -o for --via) Use `pred list` to see available problems.")] Path { @@ -126,9 +126,10 @@ Use `pred list` to see available problems.")] /// Target problem (e.g., QUBO) #[arg(value_parser = crate::problem_name::ProblemNameParser)] target: String, - /// Cost function [default: minimize-steps] - #[arg(long, default_value = "minimize-steps")] - cost: String, + /// Scalar cost function ('minimize-steps' or 'minimize:') for a single + /// best path. Omit to get the instance-free asymptotic Pareto front. + #[arg(long)] + cost: Option, /// Show all paths instead of just the cheapest #[arg(long)] all: bool, @@ -1273,7 +1274,8 @@ Examples: pred create MIS --graph 0-1,1-2 | pred reduce - --to QUBO # read from stdin Input: a problem JSON from `pred create`. Use - to read from stdin. -The --via path file is from `pred path -o path.json`. +The --via path file is from `pred path -o path.json` (its +top-level `path` is the best path; add --cost to pick a scalar-optimal one). When --via is given, --to is inferred from the path file. Output is a reduction bundle with source, target, and path. Use `pred solve reduced.json` to solve and map the solution back.")] diff --git a/problemreductions-cli/src/commands/graph.rs b/problemreductions-cli/src/commands/graph.rs index f37940725..f8e84a20d 100644 --- a/problemreductions-cli/src/commands/graph.rs +++ b/problemreductions-cli/src/commands/graph.rs @@ -2,9 +2,12 @@ use crate::output::OutputConfig; use crate::problem_name::{aliases_for, parse_problem_spec, resolve_problem_ref}; use anyhow::{Context, Result}; use problemreductions::registry::collect_schemas; -use problemreductions::rules::{Minimize, MinimizeSteps, ReductionGraph, TraversalFlow}; +use problemreductions::rules::{ + GrowthLabel, Minimize, MinimizeSteps, ReductionGraph, ReductionMode, ReductionPath, + TraversalFlow, +}; use problemreductions::types::ProblemSize; -use problemreductions::{big_o_normal_form, Expr}; +use problemreductions::{Expr, Growth}; use std::collections::BTreeMap; pub fn list(out: &OutputConfig) -> Result<()> { @@ -341,13 +344,12 @@ pub fn show(problem: &str, out: &OutputConfig) -> Result<()> { out.emit_with_default_name(&default_name, &text, &json) } -/// Format an expression as Big O notation using asymptotic normalization. -/// Falls back to wrapping the original expression if normalization fails. +/// Format an expression as Big O notation using the growth domain's canonical +/// renderer. Bounded classes render as `O()`; a growth the domain cannot +/// bound symbolically (nonlinear exponent / factorial) renders as the honest +/// `O(?)` marker — never the raw unreduced expression. fn big_o_of(expr: &Expr) -> String { - match big_o_normal_form(expr) { - Ok(norm) => format!("O({})", norm), - Err(_) => format!("O({})", expr), - } + Growth::from_expr(expr).to_big_o() } /// Format overhead fields as `field = O(...)` strings. @@ -487,10 +489,136 @@ fn format_path_json( }) } +/// Node-arrow summary (`A → B → C`) for a reduction path, deduplicating consecutive +/// same-name variant-cast steps. +fn path_arrow_summary(graph: &ReductionGraph, reduction_path: &ReductionPath) -> String { + let mut parts = Vec::new(); + let mut prev_name = ""; + for step in &reduction_path.steps { + if step.name != prev_name { + parts.push(fmt_node(graph, &step.name, &step.variant)); + prev_name = &step.name; + } + } + parts.join(&format!(" {} ", crate::output::fmt_outgoing("→"))) +} + +/// Text rendering of the asymptotic Pareto front: each path's step chain annotated +/// with a normalized `O(...)` per target size field (in the source's variables). +fn format_front_text( + graph: &ReductionGraph, + src_name: &str, + dst_name: &str, + front: &[(ReductionPath, GrowthLabel)], +) -> String { + let mut text = format!( + "Asymptotic Pareto front: {} path{} from {} to {}\n\ + (no --size given; each path shows its composed O(...) per {} size field)\n", + front.len(), + if front.len() == 1 { "" } else { "s" }, + src_name, + dst_name, + dst_name, + ); + for (idx, (reduction_path, label)) in front.iter().enumerate() { + text.push_str(&format!( + "\n--- {} ({} steps) ---\n{}\n", + crate::output::fmt_section(&format!("Path {}", idx + 1)), + reduction_path.len(), + path_arrow_summary(graph, reduction_path), + )); + for (field, growth) in label.fields() { + text.push_str(&format!(" {field} = {}\n", growth.to_big_o())); + } + } + text +} + +/// JSON rendering of the asymptotic Pareto front. Growth is emitted both as the +/// structured `Growth` serialization (issue #1075) and as a rendered `O(...)` string. +/// +/// The top-level `path` key carries the best front element's steps in exactly the +/// format `format_path_json` emits, so the saved envelope stays consumable by +/// `pred reduce --via` (the documented round-trip; front[0] is the deterministic +/// best path). Each front element's own step chain is under `front[i].path`. +fn format_front_json( + graph: &ReductionGraph, + src_name: &str, + dst_name: &str, + front: &[(ReductionPath, GrowthLabel)], +) -> serde_json::Value { + let paths: Vec = front + .iter() + .map(|(reduction_path, label)| { + let big_o: BTreeMap<&str, String> = label + .fields() + .iter() + .map(|(f, g)| (*f, g.to_big_o())) + .collect(); + serde_json::json!({ + "steps": reduction_path.len(), + "path": reduction_path.type_names(), + "growth": label.fields(), + "big_o": big_o, + }) + }) + .collect(); + // Reuse format_path_json for the best path to guarantee the top-level `path` + // array is byte-for-byte the shape `pred reduce --via` (load_path_file) parses. + let best = format_path_json(graph, &front[0].0); + serde_json::json!({ + "source": src_name, + "target": dst_name, + "mode": "asymptotic", + "front": paths, + "steps": best["steps"].clone(), + "path": best["path"].clone(), + }) +} + +/// Asymptotic Pareto-front mode of `pred path` (no `--size`/`--cost`): print the +/// front of asymptotically optimal reduction paths, each annotated with its composed +/// Big-O per target size field. See issue #1080 / design doc M3/F3a. +fn path_front( + graph: &ReductionGraph, + src_name: &str, + src_variant: &BTreeMap, + dst_name: &str, + dst_variant: &BTreeMap, + out: &OutputConfig, +) -> Result<()> { + let front = graph.asymptotic_front( + src_name, + src_variant, + dst_name, + dst_variant, + ReductionMode::Witness, + ); + + if front.is_empty() { + let variant_hint = variant_hint_for(graph, dst_name); + anyhow::bail!( + "No reduction path from {} to {}\n\ + {variant_hint}\n\ + Usage: pred path \n\ + Example: pred path MIS QUBO\n\n\ + Run `pred show {}` and `pred show {}` to check available reductions.", + src_name, + dst_name, + src_name, + dst_name, + ); + } + + let text = format_front_text(graph, src_name, dst_name, &front); + let json = format_front_json(graph, src_name, dst_name, &front); + out.emit_with_default_name("", &text, &json) +} + pub fn path( source: &str, target: &str, - cost: &str, + cost: Option<&str>, all: bool, max_paths: usize, out: &OutputConfig, @@ -531,6 +659,20 @@ pub fn path( ); } + // No `--cost` (and no `--all`): run the instance-free asymptotic Pareto search and + // print the front of asymptotically optimal paths (issue #1080 / design M3/F3a). + // Passing `--cost` opts into the single-best scalar mode (unchanged from #1076). + let Some(cost) = cost else { + return path_front( + &graph, + &src_ref.name, + &src_ref.variant, + &dst_ref.name, + &dst_ref.variant, + out, + ); + }; + let input_size = ProblemSize::new(vec![]); // Parse cost function once (validate before the search loop) @@ -601,7 +743,9 @@ fn path_all( max_paths: usize, out: &OutputConfig, ) -> Result<()> { - // Fetch one extra to detect truncation + // Fetch one extra to detect truncation. The library already returns paths in a + // deterministic length-first, then name+variant-signature order (see + // `find_paths_up_to_mode_bounded`), so no CLI-side sort is needed. let mut all_paths = graph.find_paths_up_to(src_name, src_variant, dst_name, dst_variant, max_paths + 1); @@ -620,28 +764,12 @@ fn path_all( ); } - // Sort by path length (shortest first) - all_paths.sort_by_key(|p| p.len()); - let truncated = all_paths.len() > max_paths; if truncated { all_paths.truncate(max_paths); } let returned = all_paths.len(); - let mut text = format!( - "Found {} paths from {} to {}:\n", - returned, src_name, dst_name - ); - for (idx, p) in all_paths.iter().enumerate() { - text.push_str(&format!("\n--- Path {} ---\n", idx + 1)); - text.push_str(&format_path_text(graph, p)); - } - if truncated { - text.push_str(&format!( - "\n(showing {max_paths} of more paths; use --max-paths to increase)\n" - )); - } let paths_json: Vec = all_paths .iter() @@ -697,12 +825,46 @@ fn path_all( serde_json::to_string_pretty(&json).context("Failed to serialize JSON")? ); } else { + // Build the (potentially expensive) text rendering only for text output; + // JSON and file modes above must never construct it (issue #1069). + let text = + render_all_paths_text(graph, &all_paths, src_name, dst_name, truncated, max_paths); println!("{text}"); } Ok(()) } +/// Render the `--all` text listing (header + per-path chains with normalized +/// Big-O overheads). Extracted so it is built only for text output and can be +/// exercised in-process by the issue-1069 regression tests without spawning the +/// binary. +fn render_all_paths_text( + graph: &ReductionGraph, + paths: &[ReductionPath], + src_name: &str, + dst_name: &str, + truncated: bool, + max_paths: usize, +) -> String { + let mut text = format!( + "Found {} paths from {} to {}:\n", + paths.len(), + src_name, + dst_name + ); + for (idx, p) in paths.iter().enumerate() { + text.push_str(&format!("\n--- Path {} ---\n", idx + 1)); + text.push_str(&format_path_text(graph, p)); + } + if truncated { + text.push_str(&format!( + "\n(showing {max_paths} of more paths; use --max-paths to increase)\n" + )); + } + text +} + pub fn export(out: &OutputConfig) -> Result<()> { let graph = ReductionGraph::new(); @@ -836,3 +998,262 @@ mod tests { assert_eq!(parts, vec!["KSAT", "3SAT", "2SAT"]); } } + +/// Regression, budget, and golden-determinism tests pinning the fix for issue +/// #1069 (`pred path --all` OOM/hang) and issue #1079 (raw-expression fallback + +/// unconditional JSON-mode text rendering). All tests run **in-process** against +/// the CLI's own private rendering helpers — no `pred` binary is spawned. +/// +/// Note on line lengths: with the growth domain (#1078) backing `big_o_of`, +/// composed overheads of long paths render to *genuine* multivariate polynomial +/// normal forms (an antichain of pairwise-incomparable monomials). These are the +/// correct, tight Big-O answers, not raw fallbacks — a degree-8 trivariate form +/// like `O(a^8 + a^6 b^2 + … + c^8)` legitimately runs several hundred chars. +/// The #1069 guarantee is *structural boundedness* (the antichain is capped at +/// `growth::ANTICHAIN_CAP = 32` terms, computed bottom-up in linear time), not a +/// fixed line-length limit, so the tests assert a genuine normal form plus a +/// generous structural bound rather than the (unachievable-for-multivariate) +/// 200-char figure from the issue text. +#[cfg(test)] +mod issue_1069_tests { + use super::{big_o_of, render_all_paths_text}; + use problemreductions::big_o_normal_form; + use problemreductions::rules::{ReductionGraph, ReductionPath}; + + /// Structural upper bound on a single rendered `O(...)` field: an antichain of + /// at most 32 terms (`ANTICHAIN_CAP`) over a handful of variables, each term a + /// short monomial. Far below #1069's ~2113-char raw-expression explosion, and + /// independent of path length — the point of the growth domain. + const RENDER_LEN_BOUND: usize = 2000; + + /// #1069's exploding path as a node-name chain (KSat → QUBO through + /// QuadraticAssignment/ILP). Used to reconstruct the path from the live graph + /// by name so the tests track inventory changes rather than hard-coding the + /// 2000+ char composed expression. + const NAMED_EXPLODING_PATH: [&str; 8] = [ + "KSatisfiability", + "Satisfiability", + "KSatisfiability", + "DecisionMinimumVertexCover", + "HamiltonianCircuit", + "QuadraticAssignment", + "ILP", + "QUBO", + ]; + + /// Reconstruct the #1069 exploding path deterministically. Uses the *complete* + /// [`ReductionGraph::find_all_paths`] enumeration (order-independent, unlike + /// `find_paths_up_to`'s `take(limit)`) and picks, among all paths whose + /// name-chain equals [`NAMED_EXPLODING_PATH`], the one with the + /// lexicographically smallest full (variant-annotated) rendering. This makes + /// the selection stable across build/inventory/link-order differences. + fn named_exploding_path(graph: &ReductionGraph) -> ReductionPath { + let src = crate::problem_name::resolve_problem_ref("KSat", graph).unwrap(); + let dst = crate::problem_name::resolve_problem_ref("QUBO", graph).unwrap(); + let all = graph.find_all_paths(&src.name, &src.variant, &dst.name, &dst.variant); + all.into_iter() + .filter(|p| p.type_names() == NAMED_EXPLODING_PATH) + .min_by_key(|p| p.to_string()) + .expect("the #1069 KSat->QUBO exploding path must exist in the graph") + } + + /// (1) Regression: reconstruct #1069's exploding KSat→QUBO path *by name* from + /// the live graph and assert every composed size field yields a **genuine + /// normal form** (the deleted raw fallback would have surfaced here as either + /// an `Err`/`O(?)` or an un-reduced multi-thousand-char string). + #[test] + fn issue_1069_named_exploding_path_normalizes() { + let graph = ReductionGraph::new(); + let path = named_exploding_path(&graph); + + let composed = graph.compose_path_overhead(&path); + assert!( + !composed.output_size.is_empty(), + "composed overhead has no size fields" + ); + + let mut saw_real_reduction = false; + for (field, expr) in &composed.output_size { + // Genuine normal form: not the removed `Err(_) => O()` path. + assert!( + big_o_normal_form(expr).is_ok(), + "field {field} did not normalize to a genuine Big-O form: {expr}" + ); + let rendered = big_o_of(expr); + assert!( + !rendered.contains("O(?)"), + "field {field} rendered as unbounded O(?): expr = {expr}" + ); + // Structurally bounded — no raw-expression explosion. + assert!( + rendered.len() < RENDER_LEN_BOUND, + "field {field} rendered {} chars (>= {RENDER_LEN_BOUND}); \ + raw fallback may have returned: {rendered}", + rendered.len() + ); + // The rendered normal form is never *longer* than the raw composed + // expression: proof that normalization (not passthrough) happened. + let raw_len = expr.to_string().len(); + assert!( + rendered.len() <= raw_len + "O()".len(), + "field {field}: rendered {} chars exceeds raw {raw_len}; \ + looks like a raw-expression fallback", + rendered.len() + ); + if raw_len + "O()".len() > rendered.len() { + saw_real_reduction = true; + } + } + // At least one field of this deep path must have been genuinely reduced by + // normalization (the whole point of #1069): otherwise the raw composed + // expression was already trivial and this is not the exploding path. + assert!( + saw_real_reduction, + "no field was reduced by normalization; not the #1069 exploding path" + ); + } + + /// (2) Whole-graph budget: rendering Big-O for **every** path of representative + /// hot pairs must finish well within the CI budget and never produce an + /// unbounded-length string. This is the "can't OOM/hang again" guard: it walks + /// the *complete* path set (`find_all_paths`), so no enumeration cap can hide a + /// runaway rendering. + #[test] + fn issue_1069_render_budget_is_bounded() { + let graph = ReductionGraph::new(); + let start = std::time::Instant::now(); + for (src, dst) in [("KSat", "QUBO"), ("MIS", "QUBO")] { + let src_ref = crate::problem_name::resolve_problem_ref(src, &graph).unwrap(); + let dst_ref = crate::problem_name::resolve_problem_ref(dst, &graph).unwrap(); + let paths = graph.find_all_paths( + &src_ref.name, + &src_ref.variant, + &dst_ref.name, + &dst_ref.variant, + ); + assert!(!paths.is_empty(), "expected paths for {src} -> {dst}"); + for path in &paths { + // Per-step overheads plus the composed overall overhead. + let per_step = graph.path_overheads(path); + let overall = graph.compose_path_overhead(path); + for oh in per_step.iter().chain(std::iter::once(&overall)) { + for (field, expr) in &oh.output_size { + let rendered = big_o_of(expr); + assert!( + rendered.len() < RENDER_LEN_BOUND, + "{src}->{dst} field {field} rendered {} chars (>= {RENDER_LEN_BOUND})", + rendered.len() + ); + } + } + } + } + let elapsed = start.elapsed(); + assert!( + elapsed < std::time::Duration::from_secs(5), + "rendering budget exceeded: {elapsed:?}" + ); + } + + /// (3) Golden determinism: the rendered text of the #1069 exploding path is + /// byte-stable (growth-term ordering is deterministic by construction, #1075). + /// Goldening a single, name-selected path (rather than the full `--all` + /// enumeration) keeps the fixture robust to build/inventory ordering while + /// still exercising the exact `format_path_text` code path `pred path --all` + /// prints. Regenerate the fixture with `REGEN_GOLDEN=1 cargo test issue_1069`. + /// + /// Negative control: swapping two terms in one rendered `O(...)` breaks the + /// byte-exact comparison — proving the check has teeth. + #[test] + fn issue_1069_golden_text_is_deterministic() { + let graph = ReductionGraph::new(); + let path = named_exploding_path(&graph); + // Exactly the per-path block `pred path KSat QUBO --all` prints for this path. + let actual = render_all_paths_text(&graph, &[path], "KSatisfiability", "QUBO", false, 0); + + let golden_path = concat!( + env!("CARGO_MANIFEST_DIR"), + "/tests/fixtures/issue_1069_ksat_qubo_all.txt" + ); + if std::env::var_os("REGEN_GOLDEN").is_some() { + std::fs::create_dir_all(std::path::Path::new(golden_path).parent().unwrap()).unwrap(); + std::fs::write(golden_path, &actual).unwrap(); + } + let golden = std::fs::read_to_string(golden_path).unwrap_or_else(|e| { + panic!("missing golden fixture {golden_path} ({e}); run REGEN_GOLDEN=1 cargo test issue_1069") + }); + + assert_eq!( + actual, golden, + "rendered text for the #1069 KSat->QUBO exploding path drifted from the \ + committed golden; if this is an intended inventory change, regenerate \ + with REGEN_GOLDEN=1" + ); + + // Negative control: corrupt the golden by swapping two top-level `+` terms + // inside the first multi-term `O(... + ...)` and assert the byte-exact + // comparison now fails. + let corrupted = swap_two_terms(&golden) + .expect("golden should contain a multi-term O(... + ...) to corrupt"); + assert_ne!(corrupted, golden, "swap produced no change"); + assert_ne!( + actual, corrupted, + "byte-exact comparison failed to detect a two-term swap (no teeth)" + ); + } + + /// Swap the first two top-level `+`-separated terms inside the first + /// multi-term `O(a + b + ...)` group in `text`. Uses balanced-paren matching + /// so inner `sqrt(...)` / `log(...)` groups do not confuse the scan, and only + /// splits on top-level ` + ` (depth 0). Returns `None` if no multi-term group + /// exists. + fn swap_two_terms(text: &str) -> Option { + let bytes = text.as_bytes(); + let mut search = 0; + while let Some(rel) = text[search..].find("O(") { + let open = search + rel; // index of 'O' + let inner_start = open + 2; // just past "O(" + let mut depth = 1usize; + let mut i = inner_start; + let mut top_pluses: Vec = Vec::new(); + while i < bytes.len() && depth > 0 { + match bytes[i] { + b'(' => depth += 1, + b')' => depth -= 1, + b'+' if depth == 1 + && i >= inner_start + 1 + && bytes[i - 1] == b' ' + && i + 1 < bytes.len() + && bytes[i + 1] == b' ' => + { + top_pluses.push(i - 1); // start of the " + " separator + } + _ => {} + } + i += 1; + } + let close = i - 1; // index of the matching ')' + if top_pluses.len() >= 1 { + let inner = &text[inner_start..close]; + let p1 = top_pluses[0] - inner_start; // offset of first " + " + let after = p1 + 3; + let (first, second, tail) = if top_pluses.len() >= 2 { + let p2 = top_pluses[1] - inner_start; + (&inner[..p1], &inner[after..p2], &inner[p2..]) + } else { + (&inner[..p1], &inner[after..], "") + }; + let swapped_inner = format!("{second} + {first}{tail}"); + if swapped_inner != inner { + let mut out = String::with_capacity(text.len()); + out.push_str(&text[..inner_start]); + out.push_str(&swapped_inner); + out.push_str(&text[close..]); + return Some(out); + } + } + search = inner_start; + } + None + } +} diff --git a/problemreductions-cli/src/main.rs b/problemreductions-cli/src/main.rs index 702199e49..5dcec2850 100644 --- a/problemreductions-cli/src/main.rs +++ b/problemreductions-cli/src/main.rs @@ -65,7 +65,7 @@ fn main() -> anyhow::Result<()> { cost, all, max_paths, - } => commands::graph::path(&source, &target, &cost, all, max_paths, &out), + } => commands::graph::path(&source, &target, cost.as_deref(), all, max_paths, &out), Commands::ExportGraph => commands::graph::export(&out), Commands::Inspect(args) => commands::inspect::inspect(&args.input, &out), Commands::Create(args) => commands::create::create(&args, &out), diff --git a/problemreductions-cli/src/mcp/tests.rs b/problemreductions-cli/src/mcp/tests.rs index f03e93dda..f5f7dec28 100644 --- a/problemreductions-cli/src/mcp/tests.rs +++ b/problemreductions-cli/src/mcp/tests.rs @@ -32,16 +32,49 @@ mod tests { #[test] fn test_find_path() { let server = McpServer::new(); - let result = server.find_path_inner("MIS", "QUBO", "minimize-steps", false, 20); + let result = server.find_path_inner("MIS", "QUBO", Some("minimize-steps"), false, 20); assert!(result.is_ok()); let json: serde_json::Value = serde_json::from_str(&result.unwrap()).unwrap(); assert!(json["path"].as_array().unwrap().len() > 0); } + #[test] + fn test_find_path_asymptotic_front() { + // No `cost` and not `all` → the asymptotic Pareto front with structured Growth. + let server = McpServer::new(); + let result = server.find_path_inner("KSatisfiability", "QUBO", None, false, 20); + assert!(result.is_ok(), "err: {:?}", result.err()); + let json: serde_json::Value = serde_json::from_str(&result.unwrap()).unwrap(); + assert_eq!(json["mode"], "asymptotic"); + let front = json["front"].as_array().unwrap(); + assert!(!front.is_empty()); + // Structured Growth serialization from issue #1075. + assert!(front[0]["growth"]["num_vars"]["Terms"].is_array()); + assert!(front[0]["big_o"]["num_vars"].is_string()); + } + + #[test] + fn test_find_path_asymptotic_front_has_top_level_path() { + // The default (no-cost) find_path envelope must also carry a top-level `path` + // step array (the best path) so it stays consumable as a reduction route. + let server = McpServer::new(); + let result = server.find_path_inner("MIS", "QUBO", None, false, 20); + assert!(result.is_ok(), "err: {:?}", result.err()); + let json: serde_json::Value = serde_json::from_str(&result.unwrap()).unwrap(); + assert_eq!(json["mode"], "asymptotic"); + let path = json["path"].as_array().expect("top-level path array"); + assert!(!path.is_empty(), "top-level path must have ≥ 1 step"); + // Each step parses as a from→to node pair with names. + let first = &path[0]; + assert!(first["from"]["name"].is_string()); + assert!(first["to"]["name"].is_string()); + assert_eq!(first["from"]["name"], "MaximumIndependentSet"); + } + #[test] fn test_find_path_all() { let server = McpServer::new(); - let result = server.find_path_inner("MIS", "QUBO", "minimize-steps", true, 20); + let result = server.find_path_inner("MIS", "QUBO", Some("minimize-steps"), true, 20); assert!(result.is_ok()); let json: serde_json::Value = serde_json::from_str(&result.unwrap()).unwrap(); // --all returns a structured envelope @@ -54,7 +87,7 @@ mod tests { #[test] fn test_find_path_all_structured_response() { let server = McpServer::new(); - let result = server.find_path_inner("MIS", "QUBO", "minimize-steps", true, 20); + let result = server.find_path_inner("MIS", "QUBO", Some("minimize-steps"), true, 20); assert!(result.is_ok()); let json: serde_json::Value = serde_json::from_str(&result.unwrap()).unwrap(); // Verify the structured envelope fields @@ -70,11 +103,101 @@ mod tests { assert!(first["overall_overhead"].is_array()); } + #[test] + fn test_find_path_all_matches_library_order() { + use crate::problem_name::resolve_problem_ref; + use problemreductions::rules::ReductionGraph; + + // MCP `--all` must delegate to the library ordering (length-first, then + // name+variant signature) with no local re-sort, so its ordered route list + // is identical to what the library returns directly. This is also what the + // CLI returns, since the CLI shares the same code path. + let max_paths = 6usize; + let server = McpServer::new(); + let result = server + .find_path_inner("KSatisfiability", "QUBO", None, true, max_paths) + .unwrap(); + let json: serde_json::Value = serde_json::from_str(&result).unwrap(); + let mcp_paths = json["paths"].as_array().unwrap(); + assert!(!mcp_paths.is_empty()); + + // Reconstruct each MCP path as a sequence of node signatures "name/v1/v2". + let node_sig = |node: &serde_json::Value| -> String { + let mut s = node["name"].as_str().unwrap().to_string(); + if let Some(vars) = node["variant"].as_object() { + // BTreeMap-like ordering: serde_json Map is insertion order, but the + // library serialized from a BTreeMap so keys are already sorted. + for v in vars.values() { + s.push('/'); + s.push_str(v.as_str().unwrap()); + } + } + s + }; + let mcp_sigs: Vec> = mcp_paths + .iter() + .map(|p| { + let steps = p["path"].as_array().unwrap(); + let mut seq = vec![node_sig(&steps[0]["from"])]; + for step in steps { + seq.push(node_sig(&step["to"])); + } + seq + }) + .collect(); + + // Reproduce the library-ordered, truncated route list the same way MCP/CLI do: + // fetch max_paths + 1 then keep the first max_paths. + let graph = ReductionGraph::new(); + let src = resolve_problem_ref("KSatisfiability", &graph).unwrap(); + let dst = resolve_problem_ref("QUBO", &graph).unwrap(); + let mut lib_paths = graph.find_paths_up_to( + &src.name, + &src.variant, + &dst.name, + &dst.variant, + max_paths + 1, + ); + lib_paths.truncate(max_paths); + let lib_sigs: Vec> = lib_paths + .iter() + .map(|p| { + p.steps + .iter() + .map(|s| { + let mut sig = s.name.clone(); + for v in s.variant.values() { + sig.push('/'); + sig.push_str(v); + } + sig + }) + .collect() + }) + .collect(); + + assert_eq!( + mcp_sigs, lib_sigs, + "MCP --all route list must equal the library-ordered list" + ); + + // And the route lengths are non-decreasing (length-first ordering). + let lens: Vec = mcp_paths + .iter() + .map(|p| p["steps"].as_u64().unwrap() as usize) + .collect(); + assert!( + lens.windows(2).all(|w| w[0] <= w[1]), + "MCP --all routes must be shortest-first, got {lens:?}" + ); + } + #[test] fn test_find_path_no_route() { let server = McpServer::new(); // Pick two problems with no path (if any). Use an unknown problem to trigger an error. - let result = server.find_path_inner("NonExistent", "QUBO", "minimize-steps", false, 20); + let result = + server.find_path_inner("NonExistent", "QUBO", Some("minimize-steps"), false, 20); assert!(result.is_err()); } diff --git a/problemreductions-cli/src/mcp/tools.rs b/problemreductions-cli/src/mcp/tools.rs index a0e2f1135..ae09ecaad 100644 --- a/problemreductions-cli/src/mcp/tools.rs +++ b/problemreductions-cli/src/mcp/tools.rs @@ -252,7 +252,7 @@ impl McpServer { &self, source: &str, target: &str, - cost: &str, + cost: Option<&str>, all: bool, max_paths: usize, ) -> anyhow::Result { @@ -260,8 +260,35 @@ impl McpServer { let src_ref = resolve_problem_ref(source, &graph)?; let dst_ref = resolve_problem_ref(target, &graph)?; + // No `cost` and not `all`: return the instance-free asymptotic Pareto front + // (issue #1080), using the structured `Growth` serialization from #1075. + if cost.is_none() && !all { + let front = graph.asymptotic_front( + &src_ref.name, + &src_ref.variant, + &dst_ref.name, + &dst_ref.variant, + ReductionMode::Witness, + ); + if front.is_empty() { + anyhow::bail!( + "No reduction path from {} to {}", + src_ref.name, + dst_ref.name + ); + } + return Ok(serde_json::to_string_pretty(&format_front_json( + &graph, + &src_ref.name, + &dst_ref.name, + &front, + ))?); + } + if all { - // Fetch one extra to detect truncation + // Fetch one extra to detect truncation. The library returns paths in a + // deterministic length-first, then name+variant-signature order, so the MCP + // and CLI `--all` outputs are the identical ordered route list; no local sort. let mut all_paths = graph.find_paths_up_to( &src_ref.name, &src_ref.variant, @@ -276,7 +303,6 @@ impl McpServer { dst_ref.name ); } - all_paths.sort_by_key(|p| p.len()); let truncated = all_paths.len() > max_paths; if truncated { @@ -298,8 +324,9 @@ impl McpServer { return Ok(serde_json::to_string_pretty(&json)?); } - // Single best path + // Single best path (an explicit `cost` was given; `all` is handled above). let input_size = ProblemSize::new(vec![]); + let cost = cost.expect("cost is Some in the single-best branch"); let cost_field: Option = if cost == "minimize-steps" { None @@ -965,11 +992,16 @@ impl McpServer { annotations(read_only_hint = true, open_world_hint = false) )] fn find_path(&self, Parameters(params): Parameters) -> Result { - let cost = params.cost.as_deref().unwrap_or("minimize-steps"); let all = params.all.unwrap_or(false); let max_paths = params.max_paths.unwrap_or(20); - self.find_path_inner(¶ms.source, ¶ms.target, cost, all, max_paths) - .map_err(|e| e.to_string()) + self.find_path_inner( + ¶ms.source, + ¶ms.target, + params.cost.as_deref(), + all, + max_paths, + ) + .map_err(|e| e.to_string()) } /// Export the full reduction graph as JSON @@ -1137,6 +1169,52 @@ fn format_path_json( }) } +/// JSON rendering of the asymptotic Pareto front for the `find_path` tool. Each path +/// carries the structured `Growth` serialization (issue #1075) plus a rendered +/// `O(...)` string per target size field. `Unknown` growth renders `O(?)`. +/// +/// The top-level `path` key carries the best front element's steps in the same shape +/// `format_path_json` emits, so the default `find_path` envelope stays consumable as a +/// reduction path (front[0] is the deterministic best path). Each front element's own +/// step chain is under `front[i].path`. +fn format_front_json( + graph: &ReductionGraph, + source: &str, + target: &str, + front: &[( + problemreductions::rules::ReductionPath, + problemreductions::rules::GrowthLabel, + )], +) -> serde_json::Value { + let paths: Vec = front + .iter() + .map(|(reduction_path, label)| { + let big_o: BTreeMap<&str, String> = label + .fields() + .iter() + .map(|(f, g)| (*f, g.to_big_o())) + .collect(); + serde_json::json!({ + "steps": reduction_path.len(), + "path": reduction_path.type_names(), + "growth": label.fields(), + "big_o": big_o, + }) + }) + .collect(); + // Reuse format_path_json for the best path so the top-level `path` array matches + // the step shape the reduce/bundle tooling consumes. + let best = format_path_json(graph, &front[0].0); + serde_json::json!({ + "source": source, + "target": target, + "mode": "asymptotic", + "front": paths, + "steps": best["steps"].clone(), + "path": best["path"].clone(), + }) +} + // --------------------------------------------------------------------------- // Instance tool helpers // --------------------------------------------------------------------------- diff --git a/problemreductions-cli/tests/cli_tests.rs b/problemreductions-cli/tests/cli_tests.rs index 7bc3f386a..3b8809c86 100644 --- a/problemreductions-cli/tests/cli_tests.rs +++ b/problemreductions-cli/tests/cli_tests.rs @@ -204,18 +204,114 @@ fn test_solve_balanced_complete_bipartite_subgraph_default_solver_uses_ilp() { #[test] fn test_path() { + // Bare `pred path` (no --cost / --size / --all) now prints the asymptotic Pareto + // front, each path annotated with O(...) per target size field. let output = pred().args(["path", "MIS", "QUBO"]).output().unwrap(); assert!(output.status.success()); let stdout = String::from_utf8(output.stdout).unwrap(); + assert!(stdout.contains("Asymptotic Pareto front"), "got: {stdout}"); assert!(stdout.contains("Path")); assert!(stdout.contains("step")); + assert!( + stdout.contains("O("), + "front should show Big-O per field, got: {stdout}" + ); +} + +/// Issue #1080 verification 1: `pred path KSatisfiability QUBO` (no `--size`) prints +/// ≥ 1 path, annotated with a normalized `O(...)` per QUBO size field, and the output +/// is byte-identical across two consecutive runs (determinism / golden behavior). +#[test] +fn test_path_asymptotic_front_deterministic() { + let run = || { + let output = pred() + .args(["path", "KSatisfiability", "QUBO"]) + .output() + .unwrap(); + assert!(output.status.success()); + String::from_utf8(output.stdout).unwrap() + }; + let first = run(); + let second = run(); + assert_eq!( + first, second, + "asymptotic front output must be deterministic" + ); + + // At least one path, with a normalized Big-O for QUBO's `num_vars` size field. + assert!(first.contains("Asymptotic Pareto front")); + assert!(first.contains("--- Path 1")); + assert!( + first.contains("num_vars = O("), + "each path must annotate QUBO's num_vars with O(...), got: {first}" + ); + + // The JSON surface carries the structured Growth serialization (issue #1075). + let json_out = pred() + .args(["path", "KSatisfiability", "QUBO", "--json"]) + .output() + .unwrap(); + assert!(json_out.status.success()); + let json: serde_json::Value = + serde_json::from_str(&String::from_utf8(json_out.stdout).unwrap()).unwrap(); + assert_eq!(json["mode"], "asymptotic"); + let front = json["front"].as_array().expect("front array"); + assert!(!front.is_empty(), "front must have ≥ 1 path"); + assert!( + front[0]["growth"]["num_vars"]["Terms"].is_array(), + "growth must serialize as structured Terms, got: {}", + front[0]["growth"] + ); + assert!(front[0]["big_o"]["num_vars"].is_string()); +} + +/// The asymptotic front reports one path per distinct growth vector, not per route. +/// `MVC → ILP` has dozens of reduction chains that compose to only a few Big-O +/// profiles; the front must collapse to that small handful with no duplicate growth +/// vectors. (Regression: before dedup this printed 32 paths, most identical.) +#[test] +fn test_path_front_dedups_by_growth_vector() { + let output = pred() + .args(["path", "MVC", "ILP", "--json"]) + .output() + .unwrap(); + assert!(output.status.success()); + let json: serde_json::Value = + serde_json::from_str(&String::from_utf8(output.stdout).unwrap()).unwrap(); + let front = json["front"].as_array().expect("front array"); + + // A proper Pareto front is a small handful (issue #1080: "typically 1–3 paths"). + assert!( + (1..=4).contains(&front.len()), + "expected 1..=4 distinct growth vectors, got {}", + front.len() + ); + // No two entries share a growth vector (the Big-O per size field). + let vectors: Vec = front.iter().map(|p| p["big_o"].to_string()).collect(); + let mut unique = vectors.clone(); + unique.sort(); + unique.dedup(); + assert_eq!( + unique.len(), + vectors.len(), + "front must not contain two entries with identical growth vectors: {vectors:?}" + ); } #[test] fn test_path_save() { let tmp = std::env::temp_dir().join("pred_test_path.json"); + // `--cost` selects the single-path save format (consumed by `reduce --via`). let output = pred() - .args(["path", "MIS", "QUBO", "-o", tmp.to_str().unwrap()]) + .args([ + "path", + "MIS", + "QUBO", + "--cost", + "minimize-steps", + "-o", + tmp.to_str().unwrap(), + ]) .output() .unwrap(); assert!(output.status.success()); @@ -1177,6 +1273,9 @@ fn test_reduce_via_path() { "path", "MIS/SimpleGraph/i32", "QUBO", + // A single concrete path (not the asymptotic front) for `reduce --via`. + "--cost", + "minimize-steps", "-o", path_file.to_str().unwrap(), ]) @@ -1216,6 +1315,99 @@ fn test_reduce_via_path() { std::fs::remove_file(&output_file).ok(); } +/// The documented round-trip: a *bare* `pred path S T -o path.json` (no `--cost`) +/// saves the asymptotic front plus a top-level best `path`, which `pred reduce --via` +/// must consume. Regression for #1080, which dropped the top-level `path`. +#[test] +fn test_reduce_via_bare_path() { + // 1. Create a small source problem (small so the target brute-force stays tiny). + let problem_file = std::env::temp_dir().join("pred_test_reduce_via_bare_in.json"); + let create_out = pred() + .args([ + "-o", + problem_file.to_str().unwrap(), + "create", + "MIS/SimpleGraph/i32", + "--graph", + "0-1,1-2,2-3", + "--weights", + "1,1,1,1", + ]) + .output() + .unwrap(); + assert!(create_out.status.success()); + + // 2. Bare path save (NO --cost): asymptotic front + best path. + let path_file = std::env::temp_dir().join("pred_test_reduce_via_bare_path.json"); + let path_out = pred() + .args([ + "path", + "MaximumIndependentSet/SimpleGraph/i32", + "QUBO", + "-o", + path_file.to_str().unwrap(), + ]) + .output() + .unwrap(); + assert!( + path_out.status.success(), + "stderr: {}", + String::from_utf8_lossy(&path_out.stderr) + ); + + // 3. Reduce via the bare path file (target inferred from the file). + let output_file = std::env::temp_dir().join("pred_test_reduce_via_bare_out.json"); + let reduce_out = pred() + .args([ + "-o", + output_file.to_str().unwrap(), + "reduce", + problem_file.to_str().unwrap(), + "--via", + path_file.to_str().unwrap(), + ]) + .output() + .unwrap(); + assert!( + reduce_out.status.success(), + "stderr: {}", + String::from_utf8_lossy(&reduce_out.stderr) + ); + let content = std::fs::read_to_string(&output_file).unwrap(); + let bundle: serde_json::Value = serde_json::from_str(&content).unwrap(); + assert_eq!(bundle["source"]["type"], "MaximumIndependentSet"); + assert_eq!(bundle["target"]["type"], "QUBO"); + + std::fs::remove_file(&problem_file).ok(); + std::fs::remove_file(&path_file).ok(); + std::fs::remove_file(&output_file).ok(); +} + +/// The bare-path envelope must expose BOTH the asymptotic `front` and a top-level +/// `path` step array (the best path) so it remains a valid `reduce --via` route file. +#[test] +fn test_path_front_envelope_has_front_and_path() { + let output = pred() + .args(["path", "MIS", "QUBO", "--json"]) + .output() + .unwrap(); + assert!(output.status.success()); + let json: serde_json::Value = + serde_json::from_str(&String::from_utf8(output.stdout).unwrap()).unwrap(); + + // Front envelope shape (asymptotic mode). + assert_eq!(json["mode"], "asymptotic"); + assert!(json["front"].as_array().is_some_and(|f| !f.is_empty())); + + // Top-level best path, in the step shape `reduce --via` parses. + let path = json["path"].as_array().expect("top-level path array"); + assert!(!path.is_empty(), "top-level path must have ≥ 1 step"); + let first = &path[0]; + assert!(first["from"]["name"].is_string(), "step needs from.name"); + assert!(first["to"]["name"].is_string(), "step needs to.name"); + assert_eq!(first["from"]["name"], "MaximumIndependentSet"); +} + #[test] fn test_reduce_via_infer_target() { // --via without --to: target is inferred from the path file @@ -1241,6 +1433,9 @@ fn test_reduce_via_infer_target() { "path", "MIS/SimpleGraph/i32", "QUBO", + // A single concrete path (not the asymptotic front) for `reduce --via`. + "--cost", + "minimize-steps", "-o", path_file.to_str().unwrap(), ]) @@ -1300,6 +1495,9 @@ fn test_reduce_via_rejects_target_variant_mismatch() { "path", "MIS/SimpleGraph/i32", "ILP/bool", + // A single concrete path (not the asymptotic front) for `reduce --via`. + "--cost", + "minimize-steps", "-o", path_file.to_str().unwrap(), ]) @@ -4805,8 +5003,12 @@ fn test_path_unknown_cost() { #[test] fn test_path_overall_overhead_text() { - // Use a multi-step path so the "Overall" section appears - let output = pred().args(["path", "KSAT/K3", "MIS"]).output().unwrap(); + // Use a multi-step path so the "Overall" section appears. `--cost` selects the + // single-best mode (the asymptotic front default does not render "Overall"). + let output = pred() + .args(["path", "KSAT/K3", "MIS", "--cost", "minimize-steps"]) + .output() + .unwrap(); assert!(output.status.success()); let stdout = String::from_utf8(output.stdout).unwrap(); assert!( @@ -4819,7 +5021,15 @@ fn test_path_overall_overhead_text() { fn test_path_overall_overhead_json() { let tmp = std::env::temp_dir().join("pred_test_path_overall.json"); let output = pred() - .args(["path", "KSAT/K3", "MIS", "-o", tmp.to_str().unwrap()]) + .args([ + "path", + "KSAT/K3", + "MIS", + "--cost", + "minimize-steps", + "-o", + tmp.to_str().unwrap(), + ]) .output() .unwrap(); assert!(output.status.success()); @@ -4847,7 +5057,15 @@ fn test_path_overall_overhead_composition() { // Step 2 (SAT→MIS): num_vertices = num_literals, num_edges = num_literals^2 // Overall: num_vertices = num_literals, num_edges = num_literals^2 let output = pred() - .args(["path", "KSAT/K3", "MIS", "-o", tmp.to_str().unwrap()]) + .args([ + "path", + "KSAT/K3", + "MIS", + "--cost", + "minimize-steps", + "-o", + tmp.to_str().unwrap(), + ]) .output() .unwrap(); assert!(output.status.success()); @@ -4932,7 +5150,7 @@ fn test_path_single_step_no_overall_text() { // Single-step path should NOT show the Overall section // MaxCut -> SpinGlass is a genuine 1-step path with matching default variants let output = pred() - .args(["path", "MaxCut", "SpinGlass"]) + .args(["path", "MaxCut", "SpinGlass", "--cost", "minimize-steps"]) .output() .unwrap(); assert!(output.status.success()); @@ -7615,6 +7833,69 @@ fn test_path_all_max_paths_truncates() { ); } +// Helper: run `pred path S T --all --max-paths N --json` and return the ordered +// list of per-path step counts. +fn path_all_step_counts(max_paths: &str) -> Vec { + let output = pred() + .args([ + "path", + "KSat", + "QUBO", + "--all", + "--max-paths", + max_paths, + "--json", + ]) + .output() + .unwrap(); + assert!( + output.status.success(), + "stderr: {}", + String::from_utf8_lossy(&output.stderr) + ); + let stdout = String::from_utf8(output.stdout).unwrap(); + let envelope: serde_json::Value = serde_json::from_str(&stdout).unwrap(); + envelope["paths"] + .as_array() + .expect("should have paths array") + .iter() + .map(|p| p["steps"].as_u64().expect("steps is a number")) + .collect() +} + +#[test] +fn test_path_all_truncates_after_sorting_not_before() { + // Regression: `--all` must enumerate length-first and truncate only after + // ordering, so a small --max-paths returns the SHORTEST routes, not whichever + // routes DFS discovered first. Compare a tightly-truncated run against a run + // with a generous budget. + let full = path_all_step_counts("500"); + assert!(full.len() > 3, "KSat->QUBO should have many routes"); + + // Full list is sorted shortest-first. + assert!( + full.windows(2).all(|w| w[0] <= w[1]), + "paths must be returned shortest-first, got {full:?}" + ); + let shortest = *full.first().unwrap(); + + let truncated = path_all_step_counts("3"); + assert!(truncated.len() <= 3); + // Truncated result is still sorted shortest-first... + assert!( + truncated.windows(2).all(|w| w[0] <= w[1]), + "truncated paths must be shortest-first, got {truncated:?}" + ); + // ...and it must include the known shortest length (the bug returned long + // early-discovered routes and dropped the short ones). + assert_eq!( + truncated[0], shortest, + "truncated result must start with the known shortest route length {shortest}" + ); + // The truncated step counts are exactly the shortest prefix of the full order. + assert_eq!(truncated.as_slice(), &full[..truncated.len()]); +} + #[test] fn test_path_all_max_paths_text_truncation_note() { let output = pred() diff --git a/problemreductions-cli/tests/fixtures/issue_1069_ksat_qubo_all.txt b/problemreductions-cli/tests/fixtures/issue_1069_ksat_qubo_all.txt new file mode 100644 index 000000000..bc21da7ab --- /dev/null +++ b/problemreductions-cli/tests/fixtures/issue_1069_ksat_qubo_all.txt @@ -0,0 +1,36 @@ +Found 1 paths from KSatisfiability to QUBO: + +--- Path 1 --- +Path (7 steps): KSatisfiability/KN → Satisfiability → KSatisfiability/K3 → DecisionMinimumVertexCover/SimpleGraph/i32 → HamiltonianCircuit/SimpleGraph → QuadraticAssignment → ILP/bool → QUBO/f64 + + Step 1: KSatisfiability/KN → Satisfiability + num_clauses = O(num_clauses) + num_vars = O(num_vars) + num_literals = O(num_literals) + + Step 2: Satisfiability → KSatisfiability/K3 + num_clauses = O(num_clauses + num_literals) + num_vars = O(num_clauses + num_literals + num_vars) + + Step 3: KSatisfiability/K3 → DecisionMinimumVertexCover/SimpleGraph/i32 + num_vertices = O(num_clauses + num_vars) + num_edges = O(num_clauses + num_vars) + k = O(num_clauses + num_vars) + + Step 4: DecisionMinimumVertexCover/SimpleGraph/i32 → HamiltonianCircuit/SimpleGraph + num_vertices = O(k + num_edges) + num_edges = O(k * num_vertices + num_edges) + + Step 5: HamiltonianCircuit/SimpleGraph → QuadraticAssignment + num_facilities = O(num_vertices) + num_locations = O(num_vertices) + + Step 6: QuadraticAssignment → ILP/bool + num_vars = O(num_facilities^2 * num_locations^2) + num_constraints = O(num_facilities^2 * num_locations^2) + + Step 7: ILP/bool → QUBO/f64 + num_vars = O(num_constraints * num_vars) + + Overall: + num_vars = O(num_clauses^2 * num_literals^2 * num_vars^4 + num_clauses^2 * num_literals^4 * num_vars^2 + num_clauses^2 * num_literals^6 + num_clauses^2 * num_vars^6 + num_clauses^4 * num_literals^2 * num_vars^2 + num_clauses^4 * num_literals^4 + num_clauses^4 * num_vars^4 + num_clauses^6 * num_literals^2 + num_clauses^6 * num_vars^2 + num_clauses^8 + num_literals^2 * num_vars^6 + num_literals^4 * num_vars^4 + num_literals^6 * num_vars^2 + num_literals^8 + num_vars^8) diff --git a/problemreductions-cli/tests/pred_sym_tests.rs b/problemreductions-cli/tests/pred_sym_tests.rs index 64b424d2a..9bf644973 100644 --- a/problemreductions-cli/tests/pred_sym_tests.rs +++ b/problemreductions-cli/tests/pred_sym_tests.rs @@ -13,11 +13,10 @@ fn test_pred_sym_parse() { } #[test] -fn test_pred_sym_canon_merge_terms() { +fn test_pred_sym_canon_subcommand_removed() { + // The exact-canonical-form engine is gone; `canon` is no longer a subcommand. let output = pred_sym().args(["canon", "n + n"]).output().unwrap(); - assert!(output.status.success()); - let stdout = String::from_utf8(output.stdout).unwrap(); - assert_eq!(stdout.trim(), "2 * n"); + assert!(!output.status.success()); } #[test] @@ -53,7 +52,10 @@ fn test_pred_sym_big_o_signed_polynomial() { #[test] fn test_pred_sym_big_o_sqrt_display() { - let output = pred_sym().args(["big-o", "2^(n^(1/2))"]).output().unwrap(); + // A fractional polynomial degree renders with sqrt notation. + // (`2^sqrt(n)` — a nonlinear exponent — is now unsupported, so use an + // in-domain sqrt input instead.) + let output = pred_sym().args(["big-o", "sqrt(n * m)"]).output().unwrap(); assert!(output.status.success()); let stdout = String::from_utf8(output.stdout).unwrap(); assert!( diff --git a/src/big_o.rs b/src/big_o.rs index 1b782d862..7941ae026 100644 --- a/src/big_o.rs +++ b/src/big_o.rs @@ -1,387 +1,23 @@ -//! Big-O asymptotic projection for canonical expressions. +//! Big-O asymptotic normal form. //! -//! Takes the output of `canonical_form()` and projects it into an -//! asymptotic growth class by dropping dominated terms and constant factors. +//! Thin wrapper over the [growth domain](crate::growth): compute the growth +//! class of an expression bottom-up (linear cost, no monomial expansion) and +//! render it back to a display [`Expr`]. Content the growth domain cannot bound +//! symbolically ([`Growth::Unknown`] — nonlinear exponents, factorials, negative +//! exponents) maps to the [`AsymptoticAnalysisError::Unsupported`] error. -use crate::canonical::canonical_form; -use crate::expr::{AsymptoticAnalysisError, CanonicalizationError, Expr}; - -#[derive(Clone, Debug)] -struct ProjectedTerm { - expr: Expr, - negative: bool, -} +use crate::expr::{AsymptoticAnalysisError, Expr}; +use crate::growth::Growth; /// Compute the Big-O normal form of an expression. /// -/// This is a two-phase pipeline: -/// 1. `canonical_form()` — exact symbolic simplification -/// 2. Asymptotic projection — drop dominated terms and constant factors -/// -/// Returns an expression representing the asymptotic growth class. +/// Returns an expression representing the asymptotic growth class, or +/// [`AsymptoticAnalysisError::Unsupported`] when the growth domain widens the +/// input to [`Growth::Unknown`]. pub fn big_o_normal_form(expr: &Expr) -> Result { - let canonical = canonical_form(expr).map_err(|e| match e { - CanonicalizationError::Unsupported(s) => AsymptoticAnalysisError::Unsupported(s), - })?; - - project_big_o(&canonical) -} - -/// Project a canonicalized expression into its Big-O growth class. -fn project_big_o(expr: &Expr) -> Result { - // Decompose into additive terms - let mut terms = Vec::new(); - collect_additive_terms(expr, &mut terms); - - // Project each term: drop constant multiplicative factors - let mut projected: Vec = Vec::new(); - for term in &terms { - if let Some(projected_term) = project_term(term)? { - projected.push(projected_term); - } - // Pure constants are dropped (asymptotically irrelevant) - } - - // Remove dominated terms - let survivors = remove_dominated_terms(projected); - - if survivors.is_empty() { - // All terms were constants → O(1) - return Ok(Expr::Const(1.0)); - } - - if let Some(negative) = survivors.iter().find(|term| term.negative) { - return Err(AsymptoticAnalysisError::Unsupported(format!( - "-1 * {}", - negative.expr - ))); - } - - // Deduplicate - let mut seen = std::collections::BTreeSet::new(); - let mut deduped = Vec::new(); - for term in survivors { - let key = term.expr.to_string(); - if seen.insert(key) { - deduped.push(term); - } - } - - // Rebuild sum - let mut result = deduped[0].expr.clone(); - for term in &deduped[1..] { - result = result + term.expr.clone(); - } - - Ok(result) -} - -fn collect_additive_terms(expr: &Expr, out: &mut Vec) { - match expr { - Expr::Add(a, b) => { - collect_additive_terms(a, out); - collect_additive_terms(b, out); - } - other => out.push(other.clone()), - } -} - -/// Project a single multiplicative term: strip constant factors. -/// Returns None if the term is a pure constant. -fn project_term(term: &Expr) -> Result, AsymptoticAnalysisError> { - if term.constant_value().is_some() { - return Ok(None); // Pure constant → dropped - } - - // Collect multiplicative factors - let mut factors = Vec::new(); - collect_multiplicative_factors(term, &mut factors); - - let mut coeff = 1.0; - let mut symbolic = Vec::new(); - for factor in &factors { - if let Some(c) = factor.constant_value() { - coeff *= c; - continue; - } - if contains_negative_exponent(factor) { - return Err(AsymptoticAnalysisError::Unsupported(term.to_string())); - } - symbolic.push(factor.clone()); - } - - if symbolic.is_empty() { - return Ok(None); - } - - let mut result = symbolic[0].clone(); - for f in &symbolic[1..] { - result = result * f.clone(); - } - - Ok(Some(ProjectedTerm { - expr: result, - negative: coeff < 0.0, - })) -} - -fn collect_multiplicative_factors(expr: &Expr, out: &mut Vec) { - match expr { - Expr::Mul(a, b) => { - collect_multiplicative_factors(a, out); - collect_multiplicative_factors(b, out); - } - other => out.push(other.clone()), - } -} - -/// Remove terms dominated by other terms using monomial comparison. -/// -/// A term `t` is dominated if there exists another term `s` such that -/// `t` grows no faster than `s` asymptotically. -fn remove_dominated_terms(terms: Vec) -> Vec { - if terms.len() <= 1 { - return terms; - } - - let mut survivors = Vec::new(); - for (i, term) in terms.iter().enumerate() { - let is_dominated = terms - .iter() - .enumerate() - .any(|(j, other)| i != j && term_dominated_by(&term.expr, &other.expr)); - if !is_dominated { - survivors.push(term.clone()); - } - } - survivors -} - -/// Check if `small` is asymptotically dominated by `big`. -/// -/// Supports three comparison strategies: -/// 1. Polynomial monomial exponent comparison (exact) -/// 2. Exponential vs subexponential / base comparison (structural) -/// 3. Numerical evaluation at two scales (for subexponential cross-class) -fn term_dominated_by(small: &Expr, big: &Expr) -> bool { - // Case 1: Both pure polynomial monomials — use exponent comparison - let small_exps = extract_var_exponents(small); - let big_exps = extract_var_exponents(big); - if let (Some(ref se), Some(ref be)) = (small_exps, big_exps) { - return polynomial_dominated(se, be); - } - - // Cross-class comparison: small's variables must be a subset of big's - let small_vars = small.variables(); - let big_vars = big.variables(); - if small_vars.is_empty() || big_vars.is_empty() || !small_vars.is_subset(&big_vars) { - return false; - } - - // Case 2: Exponential comparison - let small_has_exp = has_exponential_growth(small); - let big_has_exp = has_exponential_growth(big); - match (small_has_exp, big_has_exp) { - (false, true) => return true, // exponential dominates subexponential - (true, false) => return false, // subexponential can't dominate exponential - (true, true) => { - // Compare effective exponential bases - if let (Some(sb), Some(bb)) = (effective_exp_base(small), effective_exp_base(big)) { - if bb > sb * (1.0 + 1e-10) { - return true; - } - } - return false; - } - (false, false) => {} // both subexponential, fall through - } - - // Case 3: Both subexponential, same variables — numerical comparison - // Handles: poly vs poly*log, log vs log(log), poly vs log, etc. - if small_vars == big_vars { - return numerical_dominance_check(small, big, &small_vars); - } - - false -} - -/// Check polynomial dominance: small ≤ big component-wise with at least one strict inequality. -fn polynomial_dominated( - se: &std::collections::BTreeMap<&'static str, f64>, - be: &std::collections::BTreeMap<&'static str, f64>, -) -> bool { - let mut all_leq = true; - let mut any_strictly_less = false; - - for (var, small_exp) in se { - let big_exp = be.get(var).copied().unwrap_or(0.0); - if *small_exp > big_exp + 1e-15 { - all_leq = false; - break; - } - if *small_exp < big_exp - 1e-15 { - any_strictly_less = true; - } - } - - if all_leq { - for (var, big_exp) in be { - if !se.contains_key(var) && *big_exp > 1e-15 { - any_strictly_less = true; - } - } - } - - all_leq && any_strictly_less -} - -/// Extract variable → exponent mapping from a monomial expression. -/// Returns None for non-polynomial terms (exp, log, etc.). -fn extract_var_exponents(expr: &Expr) -> Option> { - use std::collections::BTreeMap; - let mut exps = BTreeMap::new(); - extract_var_exponents_inner(expr, &mut exps)?; - Some(exps) -} - -fn extract_var_exponents_inner( - expr: &Expr, - exps: &mut std::collections::BTreeMap<&'static str, f64>, -) -> Option<()> { - match expr { - Expr::Var(name) => { - *exps.entry(name).or_insert(0.0) += 1.0; - Some(()) - } - Expr::Pow(base, exp) => { - if let (Expr::Var(name), Some(e)) = (base.as_ref(), exp.constant_value()) { - if e < 0.0 { - return None; - } - *exps.entry(name).or_insert(0.0) += e; - Some(()) - } else { - None // Non-simple power - } - } - Expr::Mul(a, b) => { - extract_var_exponents_inner(a, exps)?; - extract_var_exponents_inner(b, exps) - } - Expr::Const(_) => Some(()), // Constants don't affect exponents - _ => None, // exp, log, sqrt → not a polynomial monomial - } -} - -fn contains_negative_exponent(expr: &Expr) -> bool { - match expr { - Expr::Pow(_, exp) => exp.constant_value().is_some_and(|e| e < 0.0), - Expr::Mul(a, b) | Expr::Add(a, b) => { - contains_negative_exponent(a) || contains_negative_exponent(b) - } - Expr::Exp(arg) | Expr::Log(arg) | Expr::Sqrt(arg) | Expr::Factorial(arg) => { - contains_negative_exponent(arg) - } - Expr::Const(_) | Expr::Var(_) => false, - } -} - -/// Check if an expression has exponential growth. -/// -/// Returns true if the expression contains `exp(var_expr)` or `c^(var_expr)` where c > 1. -fn has_exponential_growth(expr: &Expr) -> bool { - match expr { - Expr::Exp(arg) => !arg.variables().is_empty(), - Expr::Pow(base, exp) => { - base.constant_value().is_some_and(|c| c > 1.0) && !exp.variables().is_empty() - } - Expr::Mul(a, b) => has_exponential_growth(a) || has_exponential_growth(b), - _ => false, - } -} - -/// Compute the effective exponential base for growth rate comparison. -/// -/// For `c^(f(n))`, approximates the effective base as `c^(f(1))`. -/// This works correctly for linear exponents (the common case in complexity expressions). -fn effective_exp_base(expr: &Expr) -> Option { - match expr { - Expr::Exp(arg) => { - let vars = arg.variables(); - if vars.is_empty() { - None - } else { - let size = unit_problem_size(&vars); - let rate = arg.eval(&size); - Some(std::f64::consts::E.powf(rate)) - } - } - Expr::Pow(base, exp) => { - if let Some(c) = base.constant_value() { - let vars = exp.variables(); - if c > 1.0 && !vars.is_empty() { - let size = unit_problem_size(&vars); - let exp_at_1 = exp.eval(&size); - Some(c.powf(exp_at_1)) - } else { - None - } - } else { - None - } - } - Expr::Mul(a, b) => match (effective_exp_base(a), effective_exp_base(b)) { - (Some(ba), Some(bb)) => Some(ba * bb), - (Some(b), None) | (None, Some(b)) => Some(b), - (None, None) => None, - }, - _ => None, - } -} - -/// Create a `ProblemSize` with all variables set to the given value. -fn make_problem_size( - vars: &std::collections::HashSet<&'static str>, - val: usize, -) -> crate::types::ProblemSize { - crate::types::ProblemSize::new(vars.iter().map(|&v| (v, val)).collect()) -} - -/// Create a `ProblemSize` with all variables set to 1. -fn unit_problem_size(vars: &std::collections::HashSet<&'static str>) -> crate::types::ProblemSize { - make_problem_size(vars, 1) -} - -/// Check dominance numerically by evaluating at two scales. -/// -/// Returns true if `big/small` ratio is > 1 and increasing between the two -/// evaluation points, indicating `big` grows asymptotically faster. -fn numerical_dominance_check( - small: &Expr, - big: &Expr, - vars: &std::collections::HashSet<&'static str>, -) -> bool { - let size1 = make_problem_size(vars, 100); - let size2 = make_problem_size(vars, 10_000); - - let s1 = small.eval(&size1); - let b1 = big.eval(&size1); - let s2 = small.eval(&size2); - let b2 = big.eval(&size2); - - // Both must be finite and positive at both points - if !s1.is_finite() || !b1.is_finite() || !s2.is_finite() || !b2.is_finite() { - return false; - } - if s1 <= 1e-300 || b1 <= 1e-300 || s2 <= 1e-300 || b2 <= 1e-300 { - return false; - } - - let ratio1 = b1 / s1; - let ratio2 = b2 / s2; - - // Dominance: ratio is > 1 at both points and strictly increasing - ratio1 > 1.0 + 1e-10 && ratio2 > ratio1 * (1.0 + 1e-6) + Growth::from_expr(expr) + .to_expr() + .ok_or_else(|| AsymptoticAnalysisError::Unsupported(expr.to_string())) } #[cfg(test)] diff --git a/src/canonical.rs b/src/canonical.rs deleted file mode 100644 index 4f8c73ca7..000000000 --- a/src/canonical.rs +++ /dev/null @@ -1,431 +0,0 @@ -//! Exact symbolic canonicalization for `Expr`. -//! -//! Normalizes expressions into a canonical sum-of-terms form with signed -//! coefficients and deterministic ordering, without losing algebraic precision. - -use std::collections::BTreeMap; - -use crate::expr::{CanonicalizationError, Expr}; - -/// Hard cap on the number of additive terms produced while expanding an -/// expression into canonical sum-of-monomials form. -/// -/// Expanding a nested `(sum)^2 * (sum)^2` structure is exponential in nesting -/// depth: composed-path overheads that traverse quadratic-overhead reductions -/// (e.g. `QuadraticAssignment`) blow up to multi-GB of monomials and OOM/hang. -/// When the intermediate term count would exceed this cap we abandon expansion -/// and report the expression as `Unsupported`; callers (e.g. `big_o_of`) fall -/// back to printing the compact, un-expanded expression. See issue #1069. -/// -/// Legitimate overhead expressions stay far below this bound (the worst -/// non-pathological case is a few hundred terms), so this never affects normal -/// output — it only stops pathological blowups. This is a stopgap guard; the -/// symbolic system is slated for a larger rework. -const MAX_CANONICAL_TERMS: usize = 50_000; - -/// An opaque non-polynomial factor (exp, log, fractional-power base). -/// -/// Stored by its canonical string representation for deterministic ordering. -#[derive(Clone, Debug, PartialEq)] -struct OpaqueFactor { - /// The canonical string form (used for equality and ordering). - key: String, - /// The original `Expr` for reconstruction. - expr: Expr, -} - -impl Eq for OpaqueFactor {} - -impl PartialOrd for OpaqueFactor { - fn partial_cmp(&self, other: &Self) -> Option { - Some(self.cmp(other)) - } -} - -impl Ord for OpaqueFactor { - fn cmp(&self, other: &Self) -> std::cmp::Ordering { - self.key.cmp(&other.key) - } -} - -fn normalized_f64_bits(value: f64) -> u64 { - if value == 0.0 { - 0.0f64.to_bits() - } else { - value.to_bits() - } -} - -/// A single additive term: coefficient × product of canonical factors. -#[derive(Clone, Debug)] -struct CanonicalTerm { - /// Signed numeric coefficient. - coeff: f64, - /// Polynomial variable exponents (variable_name → exponent). - vars: BTreeMap<&'static str, f64>, - /// Non-polynomial opaque factors, sorted by key. - opaque: Vec, -} - -/// Try to merge a new opaque factor into an existing list using transcendental identities. -/// Returns `Some(updated_list)` if a merge happened, `None` if no identity applies. -fn try_merge_opaque(existing: &[OpaqueFactor], new: &OpaqueFactor) -> Option> { - for (i, existing_factor) in existing.iter().enumerate() { - // exp(a) * exp(b) -> exp(a + b) - if let (Expr::Exp(a), Expr::Exp(b)) = (&existing_factor.expr, &new.expr) { - let merged_arg = (**a).clone() + (**b).clone(); - let merged_expr = - Expr::Exp(Box::new(canonical_form(&merged_arg).unwrap_or(merged_arg))); - let mut result = existing.to_vec(); - result[i] = OpaqueFactor { - key: merged_expr.to_string(), - expr: merged_expr, - }; - return Some(result); - } - - // c^a * c^b -> c^(a+b) for matching positive constant base c - if let (Expr::Pow(base1, exp1), Expr::Pow(base2, exp2)) = (&existing_factor.expr, &new.expr) - { - if let (Some(c1), Some(c2)) = (base1.constant_value(), base2.constant_value()) { - if c1 > 0.0 && c2 > 0.0 && (c1 - c2).abs() < 1e-15 { - let merged_exp = (**exp1).clone() + (**exp2).clone(); - let canon_exp = canonical_form(&merged_exp).unwrap_or(merged_exp); - let merged_expr = Expr::Pow(base1.clone(), Box::new(canon_exp)); - let mut result = existing.to_vec(); - result[i] = OpaqueFactor { - key: merged_expr.to_string(), - expr: merged_expr, - }; - return Some(result); - } - } - } - } - None -} - -/// A canonical sum of terms: the exact normal form of an expression. -#[derive(Clone, Debug)] -pub(crate) struct CanonicalSum { - terms: Vec, -} - -impl CanonicalTerm { - fn constant(c: f64) -> Self { - Self { - coeff: c, - vars: BTreeMap::new(), - opaque: Vec::new(), - } - } - - fn variable(name: &'static str) -> Self { - let mut vars = BTreeMap::new(); - vars.insert(name, 1.0); - Self { - coeff: 1.0, - vars, - opaque: Vec::new(), - } - } - - fn opaque_factor(expr: Expr) -> Self { - let key = expr.to_string(); - Self { - coeff: 1.0, - vars: BTreeMap::new(), - opaque: vec![OpaqueFactor { key, expr }], - } - } - - /// Multiply two terms, applying transcendental identities: - /// - `exp(a) * exp(b) -> exp(a + b)` - /// - `c^a * c^b -> c^(a + b)` for matching constant base `c` - fn mul(&self, other: &CanonicalTerm) -> CanonicalTerm { - let coeff = self.coeff * other.coeff; - let mut vars = self.vars.clone(); - for (&v, &e) in &other.vars { - *vars.entry(v).or_insert(0.0) += e; - } - // Remove zero-exponent variables - vars.retain(|_, e| e.abs() > 1e-15); - - // Merge opaque factors with transcendental identities - let mut opaque = self.opaque.clone(); - for other_factor in &other.opaque { - if let Some(merged) = try_merge_opaque(&opaque, other_factor) { - opaque = merged; - } else { - opaque.push(other_factor.clone()); - } - } - opaque.sort(); - CanonicalTerm { - coeff, - vars, - opaque, - } - } - - /// Deterministic sort key for ordering terms in a sum. - fn sort_key(&self) -> (Vec<(&'static str, u64)>, Vec) { - let vars: Vec<_> = self - .vars - .iter() - .map(|(&k, &v)| (k, normalized_f64_bits(v))) - .collect(); - let opaque: Vec<_> = self.opaque.iter().map(|o| o.key.clone()).collect(); - (vars, opaque) - } -} - -impl CanonicalSum { - fn from_term(term: CanonicalTerm) -> Self { - Self { terms: vec![term] } - } - - fn add(mut self, other: CanonicalSum) -> Self { - self.terms.extend(other.terms); - self - } - - fn mul(&self, other: &CanonicalSum) -> CanonicalSum { - let mut terms = Vec::new(); - for a in &self.terms { - for b in &other.terms { - terms.push(a.mul(b)); - } - } - CanonicalSum { terms } - } - - /// Multiply with a guard against pathological expansion (see - /// [`MAX_CANONICAL_TERMS`]). The Cartesian product size is checked *before* - /// it is materialized, so this never allocates the blown-up vector. - fn try_mul(&self, other: &CanonicalSum) -> Result { - let product = self.terms.len().saturating_mul(other.terms.len()); - if product > MAX_CANONICAL_TERMS { - return Err(CanonicalizationError::Unsupported(format!( - "expression too large to canonicalize ({product} terms exceeds cap of {MAX_CANONICAL_TERMS})" - ))); - } - Ok(self.mul(other)) - } - - /// Merge terms with the same signature and drop zero-coefficient terms. - /// Sort the result deterministically. - fn simplify(self) -> Self { - type SortKey = (Vec<(&'static str, u64)>, Vec); - let mut groups: BTreeMap = BTreeMap::new(); - - for term in self.terms { - let key = term.sort_key(); - groups - .entry(key) - .and_modify(|existing| existing.coeff += term.coeff) - .or_insert(term); - } - - let mut terms: Vec<_> = groups - .into_values() - .filter(|t| t.coeff.abs() > 1e-15) - .collect(); - - terms.sort_by(|a, b| a.sort_key().cmp(&b.sort_key())); - - CanonicalSum { terms } - } -} - -/// Normalize an expression into its exact canonical sum-of-terms form. -/// -/// This performs exact symbolic simplification: -/// - Flattens nested Add/Mul -/// - Merges duplicate additive terms by summing coefficients -/// - Merges repeated multiplicative factors into powers -/// - Preserves signed coefficients (supports subtraction) -/// - Preserves transcendental identities: exp(a)*exp(b)=exp(a+b), etc. -/// - Produces deterministic ordering -/// -/// Does NOT drop terms or constant factors — use `big_o_normal_form()` for that. -pub fn canonical_form(expr: &Expr) -> Result { - let sum = expr_to_canonical(expr)?; - let simplified = sum.simplify(); - Ok(canonical_sum_to_expr(&simplified)) -} - -fn expr_to_canonical(expr: &Expr) -> Result { - match expr { - Expr::Const(c) => Ok(CanonicalSum::from_term(CanonicalTerm::constant(*c))), - Expr::Var(name) => Ok(CanonicalSum::from_term(CanonicalTerm::variable(name))), - Expr::Add(a, b) => { - let ca = expr_to_canonical(a)?; - let cb = expr_to_canonical(b)?; - Ok(ca.add(cb)) - } - Expr::Mul(a, b) => { - let ca = expr_to_canonical(a)?; - let cb = expr_to_canonical(b)?; - ca.try_mul(&cb) - } - Expr::Pow(base, exp) => canonicalize_pow(base, exp), - Expr::Exp(arg) => { - // Treat exp(canonicalized_arg) as an opaque factor - let inner = canonical_form(arg)?; - Ok(CanonicalSum::from_term(CanonicalTerm::opaque_factor( - Expr::Exp(Box::new(inner)), - ))) - } - Expr::Log(arg) => { - let inner = canonical_form(arg)?; - Ok(CanonicalSum::from_term(CanonicalTerm::opaque_factor( - Expr::Log(Box::new(inner)), - ))) - } - Expr::Sqrt(arg) => { - // sqrt(x) = x^0.5 — canonicalize as power - canonicalize_pow(arg, &Expr::Const(0.5)) - } - Expr::Factorial(arg) => { - let inner = canonical_form(arg)?; - Ok(CanonicalSum::from_term(CanonicalTerm::opaque_factor( - Expr::Factorial(Box::new(inner)), - ))) - } - } -} - -fn canonicalize_pow(base: &Expr, exp: &Expr) -> Result { - match (base, exp) { - // Constant base, constant exp → numeric constant - (_, _) if base.constant_value().is_some() && exp.constant_value().is_some() => { - let b = base.constant_value().unwrap(); - let e = exp.constant_value().unwrap(); - Ok(CanonicalSum::from_term(CanonicalTerm::constant(b.powf(e)))) - } - // Variable ^ constant exponent → vars map (supports fractional/negative exponents) - (Expr::Var(name), _) if exp.constant_value().is_some() => { - let e = exp.constant_value().unwrap(); - if e.abs() < 1e-15 { - return Ok(CanonicalSum::from_term(CanonicalTerm::constant(1.0))); - } - let mut vars = BTreeMap::new(); - vars.insert(*name, e); - Ok(CanonicalSum::from_term(CanonicalTerm { - coeff: 1.0, - vars, - opaque: Vec::new(), - })) - } - // Polynomial base ^ constant integer exponent → expand - (_, _) if exp.constant_value().is_some() => { - let e = exp.constant_value().unwrap(); - if e >= 0.0 && (e - e.round()).abs() < 1e-10 { - let n = e.round() as usize; - let base_sum = expr_to_canonical(base)?; - if n == 0 { - return Ok(CanonicalSum::from_term(CanonicalTerm::constant(1.0))); - } - let mut result = base_sum.clone(); - for _ in 1..n { - result = result.try_mul(&base_sum)?; - } - Ok(result) - } else { - // Fractional exponent with non-variable base → opaque - let canon_base = canonical_form(base)?; - Ok(CanonicalSum::from_term(CanonicalTerm::opaque_factor( - Expr::Pow(Box::new(canon_base), Box::new(Expr::Const(e))), - ))) - } - } - // Constant base ^ variable exponent → opaque (exponential growth) - (_, _) if base.constant_value().is_some() => { - let c = base.constant_value().unwrap(); - if (c - 1.0).abs() < 1e-15 { - return Ok(CanonicalSum::from_term(CanonicalTerm::constant(1.0))); - } - if c <= 0.0 { - return Err(CanonicalizationError::Unsupported(format!( - "{}^{}", - base, exp - ))); - } - let canon_exp = canonical_form(exp)?; - Ok(CanonicalSum::from_term(CanonicalTerm::opaque_factor( - Expr::Pow(Box::new(base.clone()), Box::new(canon_exp)), - ))) - } - // Variable base ^ variable exponent → unsupported - _ => Err(CanonicalizationError::Unsupported(format!( - "{}^{}", - base, exp - ))), - } -} - -fn canonical_sum_to_expr(sum: &CanonicalSum) -> Expr { - if sum.terms.is_empty() { - return Expr::Const(0.0); - } - - let term_exprs: Vec = sum.terms.iter().map(canonical_term_to_expr).collect(); - - let mut result = term_exprs[0].clone(); - for term in &term_exprs[1..] { - result = result + term.clone(); - } - result -} - -fn canonical_term_to_expr(term: &CanonicalTerm) -> Expr { - let mut factors: Vec = Vec::new(); - - // Add coefficient if not 1.0 (or -1.0, handled specially) - let (coeff_factor, sign) = if term.coeff < 0.0 { - (term.coeff.abs(), true) - } else { - (term.coeff, false) - }; - - let has_other_factors = !term.vars.is_empty() || !term.opaque.is_empty(); - - if (coeff_factor - 1.0).abs() > 1e-15 || !has_other_factors { - factors.push(Expr::Const(coeff_factor)); - } - - // Add variable powers - for (&var, &exp) in &term.vars { - if (exp - 1.0).abs() < 1e-15 { - factors.push(Expr::Var(var)); - } else { - factors.push(Expr::pow(Expr::Var(var), Expr::Const(exp))); - } - } - - // Add opaque factors - for opaque in &term.opaque { - factors.push(opaque.expr.clone()); - } - - let mut result = if factors.is_empty() { - Expr::Const(1.0) - } else { - let mut r = factors[0].clone(); - for f in &factors[1..] { - r = r * f.clone(); - } - r - }; - - if sign { - result = -result; - } - - result -} - -#[cfg(test)] -#[path = "unit_tests/canonical.rs"] -mod tests; diff --git a/src/expr.rs b/src/expr.rs index a880b6a09..fccbab06c 100644 --- a/src/expr.rs +++ b/src/expr.rs @@ -312,32 +312,6 @@ impl fmt::Display for AsymptoticAnalysisError { impl std::error::Error for AsymptoticAnalysisError {} -/// Error returned when exact canonicalization fails. -#[derive(Clone, Debug, PartialEq, Eq)] -pub enum CanonicalizationError { - /// Expression cannot be canonicalized (e.g., variable in both base and exponent). - Unsupported(String), -} - -impl fmt::Display for CanonicalizationError { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - match self { - Self::Unsupported(expr) => { - write!(f, "unsupported expression for canonicalization: {expr}") - } - } - } -} - -impl std::error::Error for CanonicalizationError {} - -/// Return a normalized `Expr` representing the asymptotic behavior of `expr`. -/// -/// This is now a compatibility wrapper for `big_o_normal_form()`. -pub fn asymptotic_normal_form(expr: &Expr) -> Result { - crate::big_o::big_o_normal_form(expr) -} - /// Compute factorial for non-negative values. /// /// For non-negative integers, returns the exact integer factorial. diff --git a/src/growth.rs b/src/growth.rs new file mode 100644 index 000000000..73cba8430 --- /dev/null +++ b/src/growth.rs @@ -0,0 +1,1044 @@ +//! Symbolic growth domain: a dedicated asymptotic normal form for reduction +//! overhead expressions. +//! +//! Where [`crate::canonical`] answers Big-O questions by fully expanding an +//! [`Expr`] to monomial normal form (exponential in nesting depth — the root +//! cause of issue #1069), the growth domain computes an asymptotic upper bound +//! *bottom-up* in a single pass, linear in the tree size, without ever expanding +//! nested sums. +//! +//! # Representation +//! +//! A [`GrowthTerm`] is one growth monomial +//! +//! ```text +//! ∏_v ∏_f base[f]^(coefficient[f] · v) +//! · ∏_v v^(poly[v]) · ∏_v (log v)^(logs[v]) +//! ``` +//! +//! and a [`Growth`] is an *antichain* of pairwise-incomparable dominant terms +//! (each summand of an asymptotic sum), or the absorbing [`Growth::Unknown`] +//! sentinel for content we cannot bound symbolically. +//! +//! # Semantic foundation (the trust contract) +//! +//! Every expression admitted to the domain is assumed **nonnegative** and +//! **weakly monotone** (nondecreasing in each variable) on `vars ≥ 2`. Under +//! these axioms Howell's multivariate-O inconsistencies vanish and +//! `f + g ≍ max(f, g)` up to a constant factor, which licenses +//! `add = antichain union + prune`. All bounds produced are **upper** bounds. +//! +//! Widening (always toward a valid upper bound): +//! - Subtraction `a − b ⇝ a + b`: `a - b` is stored as `Add(a, Mul(-1, b))`; +//! the constant `-1` is dropped by [`Growth::from_expr`], so `from_expr` of a +//! subtraction is exactly the union of the two operands. This also covers the +//! `sqrt((a − b)^2)` absolute-value idiom (`|a − b| ≤ a + b`). +//! - Constants and constant multipliers/divisors are dropped on entry. +//! - Exponentials with a **linear** exponent (`c^x`, `c^(r·x)`, `exp(x)`) are +//! first-class via symbolic base/coefficient factors. The original base is +//! authoritative: it is never normalized through a floating-point logarithm +//! and never reconstructed by rounding. Nonlinear exponents (`2^(n·k)`, +//! `2^sqrt(n)`), `factorial(·)`, and negative polynomial exponents widen to +//! [`Growth::Unknown`], which absorbs through every operation. +//! - [`Expr::Log`] evaluates numerically as the natural logarithm, but all fixed +//! logarithm bases greater than one have the same asymptotic class and are +//! intentionally represented by the single `log(v)` factor. +//! +//! # `Pow` note +//! +//! `Pow(base, k)` for a nonnegative constant `k` raises **each** antichain term +//! of `base` to the power `k` (scaling its exponents). This is the tight +//! asymptotic answer — `(n + m)^2 ≍ max(n, m)^2 = max(n^2, m^2)` by AM-GM, so no +//! binomial cross term is introduced — and it is what makes the widening chain +//! `sqrt((n − m)^2) ≍ n + m` hold exactly. + +use crate::expr::Expr; +use std::cmp::Ordering; +use std::collections::{BTreeMap, BTreeSet}; + +/// Maximum number of terms kept in an antichain. On overflow the antichain is +/// widened to a proven componentwise upper bound when one is representable; +/// otherwise it becomes [`Growth::Unknown`]. It is never truncated by order. +const ANTICHAIN_CAP: usize = 32; + +/// A base retained exactly as it appeared in the input expression. +#[derive(Clone, Debug, PartialEq, serde::Serialize)] +enum ExpBase { + /// A positive, finite constant expression used as the base of `Pow`. + Constant(Expr), + /// The distinguished base of the `exp(...)` AST constructor. + Natural, +} + +#[derive(serde::Deserialize)] +enum OwnedExpr { + Const(f64), + Var(String), + Add(Box, Box), + Mul(Box, Box), + Pow(Box, Box), + Exp(Box), + Log(Box), + Sqrt(Box), + Factorial(Box), +} + +impl OwnedExpr { + fn into_constant_expr(self) -> Option { + match self { + OwnedExpr::Const(value) => Some(Expr::Const(value)), + OwnedExpr::Var(name) => { + drop(name); + None + } + OwnedExpr::Add(a, b) => Some(a.into_constant_expr()? + b.into_constant_expr()?), + OwnedExpr::Mul(a, b) => Some(a.into_constant_expr()? * b.into_constant_expr()?), + OwnedExpr::Pow(base, exponent) => Some(Expr::pow( + base.into_constant_expr()?, + exponent.into_constant_expr()?, + )), + OwnedExpr::Exp(value) => Some(Expr::Exp(Box::new(value.into_constant_expr()?))), + OwnedExpr::Log(value) => Some(Expr::Log(Box::new(value.into_constant_expr()?))), + OwnedExpr::Sqrt(value) => Some(Expr::Sqrt(Box::new(value.into_constant_expr()?))), + OwnedExpr::Factorial(value) => { + Some(Expr::Factorial(Box::new(value.into_constant_expr()?))) + } + } + } +} + +impl<'de> serde::Deserialize<'de> for ExpBase { + fn deserialize(deserializer: D) -> Result + where + D: serde::Deserializer<'de>, + { + #[derive(serde::Deserialize)] + enum Repr { + Constant(OwnedExpr), + Natural, + } + + match Repr::deserialize(deserializer)? { + Repr::Natural => Ok(ExpBase::Natural), + Repr::Constant(base) => { + let base = base.into_constant_expr(); + if let Some(base) = + base.filter(|base| base.constant_value().is_some_and(|value| value.is_finite())) + { + Ok(ExpBase::Constant(base)) + } else { + Err(serde::de::Error::custom( + "symbolic exponential base must be a finite constant", + )) + } + } + } + } +} + +impl ExpBase { + fn structural_key(&self) -> String { + match self { + ExpBase::Constant(base) => format!("C{base:?}"), + ExpBase::Natural => "N".to_string(), + } + } + + /// Directly comparable base values. `Natural` uses the same `E` constant as + /// `Expr::Exp`; arbitrary constant subtrees remain structural-only. + fn directly_comparable_value(&self) -> Option { + match self { + ExpBase::Constant(Expr::Const(value)) => Some(*value), + ExpBase::Natural => Some(std::f64::consts::E), + ExpBase::Constant(_) => None, + } + } + + fn value(&self) -> f64 { + match self { + ExpBase::Constant(base) => base + .constant_value() + .expect("ExpBase::Constant must remain constant"), + ExpBase::Natural => std::f64::consts::E, + } + } + + fn coefficient_cmp(&self, a: f64, b: f64) -> Option { + let order = a.partial_cmp(&b)?; + if self.value() > 1.0 { + Some(order) + } else { + Some(order.reverse()) + } + } +} + +/// One symbolic exponential factor `base^(coefficient * variable)`. +#[derive(Clone, Debug, PartialEq, serde::Serialize, serde::Deserialize)] +struct ExpFactor { + base: ExpBase, + coefficient: f64, +} + +/// Canonical product of growing exponential factors for one variable. +#[derive(Clone, Debug, PartialEq, serde::Serialize, serde::Deserialize)] +struct ExpProduct { + factors: Vec, +} + +impl ExpProduct { + fn empty() -> Self { + ExpProduct { + factors: Vec::new(), + } + } + + fn single(base: ExpBase, coefficient: f64) -> Self { + Self::new(vec![ExpFactor { base, coefficient }]) + } + + /// Canonicalize without translating bases through a common logarithm. + fn new(factors: Vec) -> Self { + let mut combined: Vec = Vec::new(); + for factor in factors { + if factor.coefficient == 0.0 { + continue; + } + if let Some(existing) = combined.iter_mut().find(|f| f.base == factor.base) { + existing.coefficient += factor.coefficient; + } else { + combined.push(factor); + } + } + combined.retain(|factor| factor.coefficient != 0.0); + combined.sort_by_cached_key(|factor| factor.base.structural_key()); + ExpProduct { factors: combined } + } + + fn mul(&self, other: &Self) -> Self { + let mut factors = self.factors.clone(); + factors.extend(other.factors.iter().cloned()); + Self::new(factors) + } + + fn powf(&self, power: f64) -> Self { + let factors = self + .factors + .iter() + .filter_map(|factor| { + let coefficient = factor.coefficient * power; + (coefficient != 0.0).then(|| ExpFactor { + base: factor.base.clone(), + coefficient, + }) + }) + .collect(); + ExpProduct { factors } + } + + fn is_empty(&self) -> bool { + self.factors.is_empty() + } + + fn is_valid(&self) -> bool { + self.factors.iter().all(|factor| { + let base = factor.base.value(); + base.is_finite() + && base > 0.0 + && base != 1.0 + && factor.coefficient.is_finite() + && ((base > 1.0 && factor.coefficient > 0.0) + || (base < 1.0 && factor.coefficient < 0.0)) + }) + } + + /// Prove an ordering using only structural cancellation and direct constant + /// comparisons. `None` means "not proved", never "equal". + fn cmp_proven(&self, other: &Self) -> Option { + if self == other { + return Some(Ordering::Equal); + } + + let mut left_count = 0; + let mut right_count = 0; + let mut left_single: Option<(&ExpBase, f64)> = None; + let mut right_single: Option<(&ExpBase, f64)> = None; + + for a in &self.factors { + if let Some(b) = other.factors.iter().find(|b| a.base == b.base) { + match a.base.coefficient_cmp(a.coefficient, b.coefficient)? { + Ordering::Equal => {} + Ordering::Greater => { + left_count += 1; + left_single = Some((&a.base, a.coefficient - b.coefficient)); + } + Ordering::Less => { + right_count += 1; + right_single = Some((&a.base, b.coefficient - a.coefficient)); + } + } + } else { + left_count += 1; + left_single = Some((&a.base, a.coefficient)); + } + } + + for b in &other.factors { + if !self.factors.iter().any(|a| a.base == b.base) { + right_count += 1; + right_single = Some((&b.base, b.coefficient)); + } + } + + match (left_count, right_count) { + (0, 0) => Some(Ordering::Equal), + (0, _) => Some(Ordering::Less), + (_, 0) => Some(Ordering::Greater), + (1, 1) => { + let (a_base, a_coefficient) = left_single?; + let (b_base, b_coefficient) = right_single?; + Self::cmp_single_factor(a_base, a_coefficient, b_base, b_coefficient) + } + _ => None, + } + } + + fn cmp_single_factor( + a_base: &ExpBase, + a_coefficient: f64, + b_base: &ExpBase, + b_coefficient: f64, + ) -> Option { + if a_base == b_base { + return a_base.coefficient_cmp(a_coefficient, b_coefficient); + } + + let (a_base, b_base) = ( + a_base.directly_comparable_value()?, + b_base.directly_comparable_value()?, + ); + if a_coefficient == b_coefficient { + let base_order = a_base.partial_cmp(&b_base)?; + return if a_coefficient > 0.0 { + Some(base_order) + } else { + Some(base_order.reverse()) + }; + } + + if a_base > 1.0 && b_base > 1.0 { + match ( + a_base.partial_cmp(&b_base)?, + a_coefficient.partial_cmp(&b_coefficient)?, + ) { + (Ordering::Greater | Ordering::Equal, Ordering::Greater | Ordering::Equal) => { + Some(Ordering::Greater) + } + (Ordering::Less | Ordering::Equal, Ordering::Less | Ordering::Equal) => { + Some(Ordering::Less) + } + _ => None, + } + } else if a_base < 1.0 && b_base < 1.0 { + match ( + a_base.partial_cmp(&b_base)?, + a_coefficient.partial_cmp(&b_coefficient)?, + ) { + (Ordering::Less | Ordering::Equal, Ordering::Less | Ordering::Equal) => { + Some(Ordering::Greater) + } + (Ordering::Greater | Ordering::Equal, Ordering::Greater | Ordering::Equal) => { + Some(Ordering::Less) + } + _ => None, + } + } else { + None + } + } + + /// Approximate common-base rate used only to order search work. It is not + /// stored and never participates in equality, dominance, pruning, widening, + /// serialization, or rendering. + fn log2_estimate(&self) -> f64 { + self.factors + .iter() + .map(|factor| factor.coefficient * factor.base.value().log2()) + .sum() + } + + fn sort_key(&self) -> String { + self.factors + .iter() + .map(|factor| format!("{}={:?}", factor.base.structural_key(), factor.coefficient)) + .collect::>() + .join(",") + } +} + +/// One growth monomial, e.g. `2^(3k) · n^2 · m · log(n)`. +/// +/// Empty maps represent `O(1)`. +#[derive(Clone, Debug, PartialEq, serde::Serialize)] +pub struct GrowthTerm { + /// Variable → canonical product of symbolic exponential factors. + exp: BTreeMap<&'static str, ExpProduct>, + /// variable → polynomial degree (`0.5` covers `sqrt`). + poly: BTreeMap<&'static str, f64>, + /// variable → log power. + logs: BTreeMap<&'static str, u32>, +} + +/// The asymptotic growth class of an [`Expr`]. +#[derive(Clone, Debug, PartialEq, serde::Serialize, serde::Deserialize)] +pub enum Growth { + /// Antichain of pairwise-incomparable dominant terms, sorted by a + /// deterministic total order for platform-stable output/serialization. + Terms(Vec), + /// Absorbing sentinel: exp/factorial/negative exponents, or cap overflow + /// that even widening cannot represent. Absorbs through all operations. + Unknown, +} + +impl GrowthTerm { + /// The `O(1)` term (all maps empty). + fn one() -> Self { + GrowthTerm { + exp: BTreeMap::new(), + poly: BTreeMap::new(), + logs: BTreeMap::new(), + } + } + + /// A deterministic, platform-stable total-order key. `{v:?}` renders an + /// `f64` at full precision and is stable across platforms. + fn sort_key(&self) -> String { + let mut s = String::new(); + for (k, v) in &self.exp { + s.push('E'); + s.push_str(k); + s.push('='); + s.push_str(&v.sort_key()); + s.push(';'); + } + s.push('|'); + for (k, v) in &self.poly { + s.push('P'); + s.push_str(k); + s.push('='); + s.push_str(&format!("{v:?}")); + s.push(';'); + } + s.push('|'); + for (k, v) in &self.logs { + s.push('L'); + s.push_str(k); + s.push('='); + s.push_str(&v.to_string()); + s.push(';'); + } + s + } + + /// Raise this term to a nonnegative real power `k` (scale every exponent). + /// Log powers are `u32`; a fractional result is rounded **up** (a valid + /// upper bound, since `(log v)^p ≤ (log v)^⌈p⌉` for `v ≥ 2`). + fn powf(&self, k: f64) -> GrowthTerm { + let mut r = GrowthTerm::one(); + for (v, product) in &self.exp { + let product = product.powf(k); + if !product.is_empty() { + r.exp.insert(v, product); + } + } + for (v, deg) in &self.poly { + r.poly.insert(v, deg * k); + } + for (v, p) in &self.logs { + r.logs.insert(v, ((*p as f64) * k).ceil() as u32); + } + r + } + + /// Multiply two monomials (add matching exponents). + fn mul(&self, other: &GrowthTerm) -> GrowthTerm { + let mut t = self.clone(); + for (k, product) in &other.exp { + let combined = t + .exp + .get(k) + .map_or_else(|| product.clone(), |current| current.mul(product)); + if combined.is_empty() { + t.exp.remove(k); + } else { + t.exp.insert(k, combined); + } + } + for (k, v) in &other.poly { + *t.poly.entry(k).or_insert(0.0) += *v; + } + for (k, v) in &other.logs { + *t.logs.entry(k).or_insert(0) += *v; + } + t + } + + /// Partial order on terms: `Some(Greater)` iff `self` dominates `other` + /// (`≥` on every variable and `>` on at least one). Per variable, + /// exponential products are compared only when a symbolic proof succeeds; + /// polynomial degree and log power then break proven exponential ties. + /// Returns `None` for incomparable or unproved terms. + fn cmp(&self, other: &GrowthTerm) -> Option { + let mut vars: BTreeSet<&'static str> = BTreeSet::new(); + for m in [&self.exp, &other.exp] { + vars.extend(m.keys().copied()); + } + for m in [&self.poly, &other.poly] { + vars.extend(m.keys().copied()); + } + for m in [&self.logs, &other.logs] { + vars.extend(m.keys().copied()); + } + + let mut saw_gt = false; + let mut saw_lt = false; + let empty_exp = ExpProduct::empty(); + for v in &vars { + let exp_a = self.exp.get(v).unwrap_or(&empty_exp); + let exp_b = other.exp.get(v).unwrap_or(&empty_exp); + let exp_order = exp_a.cmp_proven(exp_b)?; + let order = if exp_order == Ordering::Equal { + self.poly + .get(v) + .copied() + .unwrap_or(0.0) + .partial_cmp(&other.poly.get(v).copied().unwrap_or(0.0))? + .then( + self.logs + .get(v) + .copied() + .unwrap_or(0) + .cmp(&other.logs.get(v).copied().unwrap_or(0)), + ) + } else { + exp_order + }; + match order { + Ordering::Greater => saw_gt = true, + Ordering::Less => saw_lt = true, + Ordering::Equal => {} + } + } + match (saw_gt, saw_lt) { + (true, true) => None, + (true, false) => Some(Ordering::Greater), + (false, true) => Some(Ordering::Less), + (false, false) => Some(Ordering::Equal), + } + } + + /// `true` iff `self` dominates `other` (grows at least as fast, and strictly + /// faster on at least one variable). + fn dominates(&self, other: &GrowthTerm) -> bool { + matches!(self.cmp(other), Some(Ordering::Greater)) + } + + /// `true` iff `self` dominates `other` or is asymptotically equal to it. + fn dominates_or_eq(&self, other: &GrowthTerm) -> bool { + matches!( + self.cmp(other), + Some(Ordering::Greater) | Some(Ordering::Equal) + ) + } + + /// A monotone scalar summary of this monomial's growth rate. Exponential rate + /// dominates polynomial degree, which dominates log power. Bigger ⇒ grows + /// faster. Used only as a search-ordering / branch-and-bound heuristic, never + /// for asymptotic dominance decisions (those go through [`GrowthTerm::cmp`]). + fn magnitude(&self) -> f64 { + let e: f64 = self.exp.values().map(ExpProduct::log2_estimate).sum(); + let p: f64 = self.poly.values().sum(); + let l: f64 = self.logs.values().map(|&x| x as f64).sum(); + 1e6 * e + p + 1e-3 * l + } +} + +impl Growth { + /// Compute the growth class of an expression in a single bottom-up pass. + pub fn from_expr(expr: &Expr) -> Growth { + // Any wholly constant subexpression is O(1). Handling it up front keeps + // constant idioms (`n / 2` = `n * 2^(-1)`, `factorial(3)`, `2^3`) out of + // the negative-exponent / factorial `Unknown` bails below. + if expr.constant_value().is_some() { + return Growth::Terms(vec![GrowthTerm::one()]); + } + match expr { + // A pure constant is O(1) — the empty term (also caught above). + Expr::Const(_) => Growth::Terms(vec![GrowthTerm::one()]), + Expr::Var(v) => { + let mut t = GrowthTerm::one(); + t.poly.insert(*v, 1.0); + Growth::Terms(vec![t]) + } + Expr::Add(a, b) => add(Growth::from_expr(a), Growth::from_expr(b)), + Expr::Mul(a, b) => mul(Growth::from_expr(a), Growth::from_expr(b)), + Expr::Pow(base, exp) => pow_expr(base, exp), + Expr::Exp(a) => exponential(ExpBase::Natural, a), + Expr::Log(a) => log_growth(Growth::from_expr(a)), + Expr::Sqrt(a) => pow_const(Growth::from_expr(a), 0.5), + Expr::Factorial(_) => Growth::Unknown, + } + } + + /// Partial order: `true` iff `self` grows at least as fast as `other`. + /// + /// Per the growth-rate reading, [`Growth::Unknown`] is the top element (it + /// may be arbitrarily large, e.g. a factorial), so it dominates everything + /// and nothing known dominates it. For two term antichains, `self` + /// dominates `other` iff every term of `other` is dominated-or-equal by + /// some term of `self` — the standard antichain (Pareto) comparison. + pub fn dominates(&self, other: &Growth) -> bool { + match (self, other) { + (Growth::Unknown, _) => true, + (Growth::Terms(_), Growth::Unknown) => false, + (Growth::Terms(a), Growth::Terms(b)) => { + b.iter().all(|tb| a.iter().any(|ta| ta.dominates_or_eq(tb))) + } + } + } + + /// A deterministic, monotone scalar summary of this growth class (the maximum + /// over its antichain terms). Exponential rate ≫ polynomial degree ≫ log + /// power; [`Growth::Unknown`] maps to a very large finite value so undecidable + /// growth sorts last. This is a *search-ordering* heuristic only (frontier + /// order, branch-and-bound bound); asymptotic dominance is decided exactly by + /// [`Growth::dominates`], never by this scalar. + pub fn magnitude(&self) -> f64 { + match self { + // Large but finite (and well below f64::MAX so sums stay finite). + Growth::Unknown => 1e18, + Growth::Terms(terms) => terms.iter().map(GrowthTerm::magnitude).fold(0.0, f64::max), + } + } + + /// Render this growth class back to a display [`Expr`] (a sum of monomials), + /// or `None` for [`Growth::Unknown`]. Terms are already in the deterministic + /// sort order, so the rendered expression is platform-stable. + /// + /// Exponential factors are rendered directly from their authoritative + /// symbolic bases and coefficients; no base reconstruction is performed. + pub fn to_expr(&self) -> Option { + match self { + Growth::Unknown => None, + Growth::Terms(terms) => { + if terms.is_empty() { + return Some(Expr::Const(1.0)); + } + let mut it = terms.iter().map(term_to_expr); + let mut acc = it.next().unwrap(); + for e in it { + acc = acc + e; + } + Some(acc) + } + } + } + + /// Canonical Big-O string for this growth class: `O()` for a bounded + /// class, or `O(?)` for [`Growth::Unknown`] (no honest asymptotic bound — + /// nonlinear exponent or factorial). This is the single source of truth for + /// how a growth is displayed as Big-O; presentation layers must call it rather + /// than re-deriving the mapping (and the `Unknown` spelling) themselves. + pub fn to_big_o(&self) -> String { + match self.to_expr() { + Some(e) => format!("O({e})"), + None => "O(?)".to_string(), + } + } +} + +/// Render one monomial as a product of its factors (or `Const(1)` when empty). +fn term_to_expr(t: &GrowthTerm) -> Expr { + let mut factors: Vec = Vec::new(); + for (v, product) in &t.exp { + factors.extend(product.factors.iter().map(|factor| exp_factor(v, factor))); + } + for (v, deg) in &t.poly { + factors.push(poly_factor(v, *deg)); + } + for (v, power) in &t.logs { + factors.push(log_factor(v, *power)); + } + let mut it = factors.into_iter(); + match it.next() { + None => Expr::Const(1.0), + Some(first) => it.fold(first, |acc, f| acc * f), + } +} + +/// Render a stored exponential factor without changing its base or coefficient. +fn exp_factor(v: &'static str, factor: &ExpFactor) -> Expr { + let exponent = if factor.coefficient == 1.0 { + Expr::Var(v) + } else { + Expr::Const(factor.coefficient) * Expr::Var(v) + }; + match &factor.base { + ExpBase::Constant(base) => Expr::pow(base.clone(), exponent), + ExpBase::Natural => Expr::Exp(Box::new(exponent)), + } +} + +/// Render `v^degree` (`Display` turns degree `0.5` into `sqrt(v)`). +fn poly_factor(v: &'static str, degree: f64) -> Expr { + if degree == 1.0 { + Expr::Var(v) + } else { + Expr::pow(Expr::Var(v), Expr::Const(degree)) + } +} + +/// Render `(log v)^power`. +fn log_factor(v: &'static str, power: u32) -> Expr { + let log = Expr::Log(Box::new(Expr::Var(v))); + if power == 1 { + log + } else { + Expr::pow(log, Expr::Const(power as f64)) + } +} + +/// Prune a bag of terms to its maximal antichain: drop any term dominated by +/// another and collapse exact duplicates. The resulting *set* is independent of +/// input order. +fn prune(mut terms: Vec) -> Vec { + // Proven-equal terms can retain different symbolic spellings (for example, + // `exp(n)` and a literal-e base). Sort first so the representative does not + // depend on operand order. + terms.sort_by_cached_key(GrowthTerm::sort_key); + let mut result: Vec = Vec::new(); + for t in terms { + if result.iter().any(|r| r.dominates_or_eq(&t)) { + continue; + } + result.retain(|r| !t.dominates(r)); + result.push(t); + } + result +} + +/// Construct a componentwise upper bound when every exponential component has +/// a symbolically proven maximal product. +fn componentwise_max(terms: &[GrowthTerm]) -> Option { + let mut m = GrowthTerm::one(); + let mut vars = BTreeSet::new(); + for term in terms { + vars.extend(term.exp.keys().copied()); + vars.extend(term.poly.keys().copied()); + vars.extend(term.logs.keys().copied()); + } + + for var in vars { + let empty_exp = ExpProduct::empty(); + let mut maximum = &empty_exp; + for product in terms + .iter() + .map(|term| term.exp.get(var).unwrap_or(&empty_exp)) + { + if matches!(product.cmp_proven(maximum), Some(Ordering::Greater)) { + maximum = product; + } + } + if !terms.iter().all(|term| { + matches!( + maximum.cmp_proven(term.exp.get(var).unwrap_or(&empty_exp)), + Some(Ordering::Greater | Ordering::Equal) + ) + }) { + return None; + } + if !maximum.is_empty() { + m.exp.insert(var, maximum.clone()); + } + + let mut max_poly = 0.0_f64; + let mut max_logs = 0_u32; + for term in terms { + let degree = term.poly.get(var).copied().unwrap_or(0.0); + if !degree.is_finite() { + return None; + } + max_poly = max_poly.max(degree); + max_logs = max_logs.max(term.logs.get(var).copied().unwrap_or(0)); + } + if max_poly > 0.0 { + m.poly.insert(var, max_poly); + } + if max_logs > 0 { + m.logs.insert(var, max_logs); + } + } + Some(m) +} + +fn growth_term_is_valid(term: &GrowthTerm) -> bool { + term.exp + .values() + .all(|product| !product.is_empty() && product.is_valid()) + && term + .poly + .values() + .all(|degree| degree.is_finite() && *degree >= 0.0) +} + +/// Prune, apply the antichain cap (widening upward on overflow), and sort into +/// the deterministic total order. +fn make_growth(terms: Vec) -> Growth { + if !terms.iter().all(growth_term_is_valid) { + return Growth::Unknown; + } + let mut pruned = prune(terms); + if pruned.len() > ANTICHAIN_CAP { + let Some(widened) = componentwise_max(&pruned) else { + return Growth::Unknown; + }; + if !pruned.iter().all(|term| widened.dominates_or_eq(term)) { + return Growth::Unknown; + } + pruned = vec![widened]; + } + debug_assert!(pruned.iter().all(growth_term_is_valid)); + Growth::Terms(pruned) +} + +/// Antichain union (asymptotic `+ ≍ max`). +fn add(a: Growth, b: Growth) -> Growth { + match (a, b) { + (Growth::Unknown, _) | (_, Growth::Unknown) => Growth::Unknown, + (Growth::Terms(mut x), Growth::Terms(y)) => { + x.extend(y); + make_growth(x) + } + } +} + +/// Pairwise product of two antichains. +fn mul(a: Growth, b: Growth) -> Growth { + match (a, b) { + (Growth::Unknown, _) | (_, Growth::Unknown) => Growth::Unknown, + (Growth::Terms(x), Growth::Terms(y)) => { + let mut prod = Vec::with_capacity(x.len() * y.len()); + for tx in &x { + for ty in &y { + prod.push(tx.mul(ty)); + } + } + make_growth(prod) + } + } +} + +/// Raise a whole antichain to a nonnegative real power `k` (raise each term). +fn pow_const(g: Growth, k: f64) -> Growth { + match g { + Growth::Unknown => Growth::Unknown, + Growth::Terms(terms) => make_growth(terms.iter().map(|t| t.powf(k)).collect()), + } +} + +/// Transfer function for `Pow(base, exp)`. +fn pow_expr(base: &Expr, exp: &Expr) -> Growth { + if let Some(k) = exp.constant_value() { + // Constant exponent → polynomial power. + if k < 0.0 { + return Growth::Unknown; // negative exponent + } + if k == 0.0 { + return Growth::Terms(vec![GrowthTerm::one()]); // x^0 = O(1) + } + pow_const(Growth::from_expr(base), k) + } else if let Some(c) = base.constant_value() { + // Constant base, variable exponent → exponential. + if c.is_finite() { + exponential(ExpBase::Constant(base.clone()), exp) + } else { + Growth::Unknown + } + } else { + // Variable base and variable exponent (e.g. n^m) → not representable. + Growth::Unknown + } +} + +/// Transfer function for a symbolic fixed-base exponential. The base's numeric +/// value is used only for domain and monotonic-direction checks. +fn exponential(base: ExpBase, exp: &Expr) -> Growth { + let c = base.value(); + if !c.is_finite() || c <= 0.0 { + return Growth::Unknown; + } + if c == 1.0 { + // 1^x = 1 for every x: bounded by O(1). + return Growth::Terms(vec![GrowthTerm::one()]); + } + match linear_form(exp) { + None => Growth::Unknown, // nonlinear exponent + Some(coeffs) => { + let mut term = GrowthTerm::one(); + for (v, coeff) in coeffs { + if !coeff.is_finite() { + return Growth::Unknown; + } + // Drop decaying directions as an upward widening. A fractional + // base grows only along negative exponent coefficients. + if (c > 1.0 && coeff > 0.0) || (c < 1.0 && coeff < 0.0) { + term.exp.insert(v, ExpProduct::single(base.clone(), coeff)); + } + } + make_growth(vec![term]) + } + } +} + +/// Extract the linear coefficients of an expression (variable → coefficient), +/// or `None` if the expression is not linear in its variables. The additive +/// constant term is ignored (dropped). Pure constants map to the empty form. +fn linear_form(expr: &Expr) -> Option> { + if expr.constant_value().is_some() { + return Some(BTreeMap::new()); + } + match expr { + Expr::Var(v) => { + let mut m = BTreeMap::new(); + m.insert(*v, 1.0); + Some(m) + } + Expr::Add(a, b) => { + let mut m = linear_form(a)?; + for (k, v) in linear_form(b)? { + *m.entry(k).or_insert(0.0) += v; + } + Some(m) + } + Expr::Mul(a, b) => { + // A linear term times a variable is nonlinear, so one side must be + // a constant scalar. + if let Some(c) = a.constant_value() { + Some( + linear_form(b)? + .into_iter() + .map(|(k, v)| (k, v * c)) + .collect(), + ) + } else if let Some(c) = b.constant_value() { + Some( + linear_form(a)? + .into_iter() + .map(|(k, v)| (k, v * c)) + .collect(), + ) + } else { + None + } + } + // Pow / Exp / Log / Sqrt / Factorial of variables are nonlinear. + _ => None, + } +} + +/// Transfer function for `Log(a)`: `log` of an antichain is `log` of its +/// dominant term(s), unioned. Uses `log(n^a · m^b) ≍ log n + log m` and +/// `log(2^(r·n)) ≍ n`. +fn log_growth(g: Growth) -> Growth { + match g { + Growth::Unknown => Growth::Unknown, + Growth::Terms(terms) => { + let mut out = Vec::new(); + for t in &terms { + out.extend(log_term(t)); + } + if out.is_empty() { + out.push(GrowthTerm::one()); // log(O(1)) = O(1) + } + make_growth(out) + } + } +} + +/// `log` of a single monomial, returned as its own (small) antichain of +/// summands. `log(∏ baseᵢ^(rᵢ·vᵢ) · ∏vⱼ^aⱼ · ∏(log vₖ)^sₖ)` distributes over the +/// product into a *sum* of the log of each factor, so every factor class of the +/// monomial contributes its own summand — none may be dropped (e.g. `log(2^n·m)` +/// is `n + log m`, not `n`). `make_growth`/`prune` then collapse any dominated +/// summands (so `log(2^n·n^2)` reduces back to `n`). +fn log_term(t: &GrowthTerm) -> Vec { + let mut out = Vec::new(); + // Every stored exponential product grows, so its logarithm is linear. + for v in t.exp.keys().copied() { + let mut g = GrowthTerm::one(); + g.poly.insert(v, 1.0); + out.push(g); + } + // log(v^a) ≍ log v: each positive-degree polynomial factor becomes a log. + for v in t.poly.iter().filter(|(_, d)| **d > 0.0).map(|(k, _)| *k) { + let mut g = GrowthTerm::one(); + g.logs.insert(v, 1); + out.push(g); + } + // log((log v)^s) = log log v, upper-bounded by log v (log log v ≤ log v for + // v ≥ 2): each log factor stays a single log. + for v in t.logs.keys().copied() { + let mut g = GrowthTerm::one(); + g.logs.insert(v, 1); + out.push(g); + } + // Empty term: log(O(1)) = O(1). + if out.is_empty() { + out.push(GrowthTerm::one()); + } + out +} + +// --- serde --- +// +// `GrowthTerm` uses `&'static str` keys (to align with `Expr::Var`), which serde +// cannot deserialize directly. `Deserialize` reads owned `String` keys and leaks +// them to `&'static str`, matching the convention of `Expr`'s runtime parser. +// Each unique key leaks a small allocation that is never freed; acceptable for +// the CLI's one-shot serialization, not for hot loops with adversarial input. + +impl<'de> serde::Deserialize<'de> for GrowthTerm { + fn deserialize(deserializer: D) -> Result + where + D: serde::Deserializer<'de>, + { + #[derive(serde::Deserialize)] + struct Repr { + exp: BTreeMap, + poly: BTreeMap, + logs: BTreeMap, + } + fn leak(s: String) -> &'static str { + Box::leak(s.into_boxed_str()) + } + let r = Repr::deserialize(deserializer)?; + let term = GrowthTerm { + exp: r + .exp + .into_iter() + .map(|(k, product)| (leak(k), ExpProduct::new(product.factors))) + .collect(), + poly: r.poly.into_iter().map(|(k, v)| (leak(k), v)).collect(), + logs: r.logs.into_iter().map(|(k, v)| (leak(k), v)).collect(), + }; + if growth_term_is_valid(&term) { + Ok(term) + } else { + Err(serde::de::Error::custom("invalid symbolic growth term")) + } + } +} + +#[cfg(test)] +#[path = "unit_tests/growth.rs"] +mod tests; diff --git a/src/lib.rs b/src/lib.rs index 2083070c1..4cad72e55 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -20,13 +20,16 @@ extern crate self as problemreductions; pub(crate) mod big_o; -pub(crate) mod canonical; pub mod config; pub mod error; #[cfg(feature = "example-db")] pub mod example_db; pub mod export; pub(crate) mod expr; +// The growth domain backs `big_o_normal_form` (M2) and the asymptotic Pareto path +// search (`GrowthLabel`, M3/F3a). `Growth` is re-exported for CLI/MCP consumers that +// render or serialize the asymptotic front. +pub mod growth; pub mod io; pub mod models; pub mod registry; @@ -110,9 +113,9 @@ pub mod prelude { // Re-export commonly used items at crate root pub use big_o::big_o_normal_form; -pub use canonical::canonical_form; pub use error::{ProblemError, Result}; -pub use expr::{asymptotic_normal_form, AsymptoticAnalysisError, CanonicalizationError, Expr}; +pub use expr::{AsymptoticAnalysisError, Expr}; +pub use growth::Growth; pub use registry::{ComplexityClass, ProblemInfo}; pub use solvers::{BruteForce, Solver}; pub use traits::Problem; diff --git a/src/rules/analysis.rs b/src/rules/analysis.rs index 6d616877d..a54d2dc5f 100644 --- a/src/rules/analysis.rs +++ b/src/rules/analysis.rs @@ -1,14 +1,15 @@ //! Analysis utilities for the reduction graph. //! //! Detects primitive reduction rules that are dominated by composite paths, -//! using asymptotic normalization plus monomial-dominance comparison. +//! comparing overhead expressions through the shared symbolic growth domain +//! ([`crate::growth::Growth`]). //! //! This analysis is **sound but incomplete**: it reports `Dominated` only when -//! the symbolic comparison is trustworthy, and `Unknown` when metadata is too -//! weak to compare safely. +//! the growth comparison is trustworthy, and `Unknown` when a field's growth is +//! [`Growth::Unknown`] (nonlinear exponent, factorial, …). -use crate::canonical::canonical_form; use crate::expr::Expr; +use crate::growth::Growth; use crate::rules::graph::{ReductionGraph, ReductionPath}; use crate::rules::registry::ReductionOverhead; use std::collections::{BTreeMap, BTreeSet}; @@ -94,184 +95,22 @@ pub fn format_problem_variant(name: &str, variant: &BTreeMap) -> format!("{name} {{{vars}}}") } -// ────────── Polynomial normalization ────────── - -/// A monomial: coefficient × ∏(variable ^ exponent). -#[derive(Debug, Clone)] -struct Monomial { - coeff: f64, - /// Variable name → exponent. Only non-zero exponents stored. - vars: BTreeMap<&'static str, f64>, -} - -impl Monomial { - fn constant(c: f64) -> Self { - Self { - coeff: c, - vars: BTreeMap::new(), - } - } - - fn variable(name: &'static str) -> Self { - let mut vars = BTreeMap::new(); - vars.insert(name, 1.0); - Self { coeff: 1.0, vars } - } - - /// Multiply two monomials. - fn mul(&self, other: &Monomial) -> Monomial { - let coeff = self.coeff * other.coeff; - let mut vars = self.vars.clone(); - for (&v, &e) in &other.vars { - *vars.entry(v).or_insert(0.0) += e; - } - Monomial { coeff, vars } - } -} - -/// A polynomial (sum of monomials) in normal form. -#[derive(Debug, Clone)] -struct NormalizedPoly { - terms: Vec, -} - -impl NormalizedPoly { - fn add(mut self, other: NormalizedPoly) -> NormalizedPoly { - self.terms.extend(other.terms); - self - } - - fn mul(&self, other: &NormalizedPoly) -> NormalizedPoly { - let mut terms = Vec::new(); - for a in &self.terms { - for b in &other.terms { - terms.push(a.mul(b)); - } - } - NormalizedPoly { terms } - } - - /// True if any monomial has a negative coefficient. - fn has_negative_coefficients(&self) -> bool { - self.terms.iter().any(|m| m.coeff < -1e-15) - } -} - -/// Normalize an expression into a sum of monomials. -/// -/// Supports: constants, variables, addition, multiplication, -/// and powers with non-negative constant exponents. -/// Returns `Err` for exp, log, sqrt, division, and negative exponents. -fn normalize_polynomial(expr: &Expr) -> Result { - match expr { - Expr::Const(c) => Ok(NormalizedPoly { - terms: vec![Monomial::constant(*c)], - }), - Expr::Var(v) => Ok(NormalizedPoly { - terms: vec![Monomial::variable(v)], - }), - Expr::Add(a, b) => { - let pa = normalize_polynomial(a)?; - let pb = normalize_polynomial(b)?; - Ok(pa.add(pb)) - } - Expr::Mul(a, b) => { - let pa = normalize_polynomial(a)?; - let pb = normalize_polynomial(b)?; - Ok(pa.mul(&pb)) - } - Expr::Pow(base, exp) => { - if let Expr::Const(c) = exp.as_ref() { - if *c < 0.0 { - return Err(format!("negative exponent: {c}")); - } - let pb = normalize_polynomial(base)?; - // Single monomial: multiply exponents - if pb.terms.len() == 1 { - let m = &pb.terms[0]; - let coeff = m.coeff.powf(*c); - let vars: BTreeMap<_, _> = m.vars.iter().map(|(&v, &e)| (v, e * c)).collect(); - return Ok(NormalizedPoly { - terms: vec![Monomial { coeff, vars }], - }); - } - // Multi-term polynomial raised to non-negative integer power - let n = *c as usize; - if c.fract().abs() < 1e-10 { - if n == 0 { - return Ok(NormalizedPoly { - terms: vec![Monomial::constant(1.0)], - }); - } - let mut result = pb.clone(); - for _ in 1..n { - result = result.mul(&pb); - } - return Ok(result); - } - Err(format!( - "non-integer power of multi-term polynomial: ({base})^{c}" - )) - } else { - Err(format!("variable exponent: ({base})^({exp})")) - } - } - Expr::Exp(_) => Err("exp() not supported".into()), - Expr::Log(_) => Err("log() not supported".into()), - Expr::Sqrt(_) => Err("sqrt() not supported".into()), - Expr::Factorial(_) => Err("factorial() not supported".into()), - } -} - -fn prepare_expr_for_comparison(expr: &Expr) -> Expr { - canonical_form(expr).unwrap_or_else(|_| expr.clone()) -} - -// ────────── Monomial-dominance comparison ────────── - -/// Check if monomial `small` is asymptotically dominated by monomial `big`. -/// -/// True iff for every variable in `small`, `big` has at least as large an exponent. -/// This means `small` grows no faster than `big` as all variables → ∞. -fn monomial_dominated_by(small: &Monomial, big: &Monomial) -> bool { - for (&var, &exp_small) in &small.vars { - let exp_big = big.vars.get(var).copied().unwrap_or(0.0); - if exp_small > exp_big + 1e-10 { - return false; - } - } - true -} - -/// Check if polynomial `a` is asymptotically ≤ polynomial `b`. -/// -/// True iff every positive-coefficient monomial in `a` is dominated by -/// some positive-coefficient monomial in `b`. -fn poly_leq(a: &NormalizedPoly, b: &NormalizedPoly) -> bool { - let b_positive: Vec<&Monomial> = b.terms.iter().filter(|m| m.coeff > 1e-15).collect(); - - for a_term in &a.terms { - if a_term.coeff <= 1e-15 { - continue; // zero or negative — can only make `a` smaller - } - let dominated = b_positive - .iter() - .any(|b_term| monomial_dominated_by(a_term, b_term)); - if !dominated { - return false; - } - } - true -} - // ────────── Overhead comparison ────────── -/// Compare two overheads across all common fields. +/// Compare two overheads across all common fields, using the shared symbolic +/// growth domain ([`Growth`]) as the single dominance order. /// -/// Returns `Dominated` if composite ≤ primitive on all common fields. -/// Returns `NotDominated` if composite is worse on any common field. -/// Returns `Unknown` if any common field's expressions cannot be normalized -/// into a comparable polynomial form or contain negative coefficients. +/// Fields present in only one overhead are skipped (common-field semantics). +/// For each common field with primitive growth `pg` and composite growth `cg`: +/// - if either is [`Growth::Unknown`] the whole comparison is `Unknown`; +/// - otherwise the field is fine iff the composite is dominated-or-equal by the +/// primitive (`pg` grows ≥ `cg`, i.e. `pg.dominates(&cg)` — reflexive, so an +/// equal field counts as fine); +/// - otherwise (composite strictly worse, or the two growths incomparable) the +/// comparison is `NotDominated`. +/// +/// Returns `Dominated` when every common field is fine and at least one common +/// field exists; `NotDominated` when there is no common field. pub fn compare_overhead( primitive: &ReductionOverhead, composite: &ReductionOverhead, @@ -290,30 +129,20 @@ pub fn compare_overhead( }; any_common = true; - let primitive_prepared = prepare_expr_for_comparison(prim_expr); - let composite_prepared = prepare_expr_for_comparison(comp_expr); - - if primitive_prepared == composite_prepared { - continue; - } - - let primitive_poly = match normalize_polynomial(&primitive_prepared) { - Ok(p) => p, - Err(_) => return ComparisonStatus::Unknown, - }; - let composite_poly = match normalize_polynomial(&composite_prepared) { - Ok(p) => p, - Err(_) => return ComparisonStatus::Unknown, - }; + let pg = Growth::from_expr(prim_expr); + let cg = Growth::from_expr(comp_expr); - // Reject expressions with negative coefficients - if primitive_poly.has_negative_coefficients() || composite_poly.has_negative_coefficients() - { + // A field whose growth we cannot bound symbolically makes the whole + // comparison undecidable. + if matches!(pg, Growth::Unknown) || matches!(cg, Growth::Unknown) { return ComparisonStatus::Unknown; } - // Check: composite ≤ primitive on this field - if !poly_leq(&composite_poly, &primitive_poly) { + // `pg.dominates(&cg)` means the primitive grows at least as fast as the + // composite on this field (composite ≤ primitive). `dominates` is + // reflexive, so asymptotically-equal fields pass here. Anything else — + // composite strictly worse, or the two growths incomparable — fails. + if !pg.dominates(&cg) { return ComparisonStatus::NotDominated; } } diff --git a/src/rules/cost.rs b/src/rules/cost.rs index 7678d4d87..c44846efe 100644 --- a/src/rules/cost.rs +++ b/src/rules/cost.rs @@ -6,6 +6,12 @@ use crate::types::ProblemSize; /// User-defined cost function for path optimization. pub trait PathCostFn { /// Compute cost of taking an edge given current problem size. + /// + /// Implementations **must** be monotone in `current_size` (a componentwise-larger + /// size never yields a smaller edge cost). The Pareto search prunes by `(cost, size)` + /// dominance, and this monotonicity is what gives the isotonicity that makes such + /// pruning sound. (A nonnegative cost is also expected — all shipped implementations + /// return one — though the kernel no longer branch-and-bounds on it.) fn edge_cost(&self, overhead: &ReductionOverhead, current_size: &ProblemSize) -> f64; } diff --git a/src/rules/graph.rs b/src/rules/graph.rs index ef8a27ff2..32e34cfdd 100644 --- a/src/rules/graph.rs +++ b/src/rules/graph.rs @@ -13,6 +13,9 @@ //! - JSON export for documentation and visualization use crate::rules::cost::PathCostFn; +use crate::rules::pareto::{ + CostLabel, GrowthLabel, MeasuredLabel, PathLabel, ReductionEdge, BAG_CAP, HOP_CAP, +}; use crate::rules::registry::{ AggregateReduceFn, EdgeCapabilities, ReduceFn, ReductionEntry, ReductionOverhead, }; @@ -26,6 +29,7 @@ use serde::Serialize; use std::any::Any; use std::cmp::Reverse; use std::collections::{BTreeMap, BinaryHeap, HashMap, HashSet}; +use std::rc::Rc; /// A source/target pair from the reduction graph, returned by /// [`ReductionGraph::outgoing_reductions`] and [`ReductionGraph::incoming_reductions`]. @@ -464,6 +468,11 @@ impl ReductionGraph { /// Find the cheapest path between two specific problem variants while /// requiring a specific edge capability. + /// + /// Runs the generic [Pareto label-setting search](Self::pareto_search) with a + /// scalar [`CostLabel`], reproducing Dijkstra's single-objective behavior for the + /// given [`PathCostFn`]. Returns the front's best element under the deterministic + /// tie-break (smallest cost, then fewest hops, then lexicographic node names). #[allow(clippy::too_many_arguments)] pub fn find_cheapest_path_mode( &self, @@ -477,71 +486,299 @@ impl ReductionGraph { ) -> Option { let src = self.lookup_node(source, source_variant)?; let dst = self.lookup_node(target, target_variant)?; - let node_path = self.dijkstra(src, dst, mode, input_size, cost_fn)?; - Some(self.node_path_to_reduction_path(&node_path)) + let initial = CostLabel::new(input_size.clone(), cost_fn); + let mut front = self.pareto_search(src, dst, mode, initial, false); + self.pick_best_front(&mut front).map(|(path, _)| path) } - /// Core Dijkstra search on node indices. - fn dijkstra( + /// Generic Pareto label-setting search from `src` to `dst`. + /// + /// Maintains a per-node **bag** (an antichain of non-dominated labels); a label is + /// discarded only when another label at the same node [dominates](PathLabel::dominates) + /// it. Each surviving label carries a predecessor pointer for path reconstruction. + /// Pruning is by dominance alone — always sound for any label domain, unlike a + /// branch-and-bound bound, which would require a monotone scalar `cost` that the + /// measured domain does not have. The frontier is explored in ascending + /// [`cost`](PathLabel::cost) order (a heuristic that finds good paths early). + /// Deterministic safety caps apply: [`HOP_CAP`] bounds path length, and [`BAG_CAP`] + /// bounds each bag with a deterministic tie-break (never iteration-order truncation). + /// Edges are visited in a deterministic (target-name, target-variant) order. + /// + /// When `exhaustive` is `true`, the componentwise dominance guard is disabled (bags + /// retain all labels up to the cap); the sound guards inside [`PathLabel::extend`] + /// still apply. + /// + /// Returns the Pareto front at `dst`: `(path, label)` pairs, deterministically + /// ordered by (cost, hops, node-name path). + pub(crate) fn pareto_search( &self, src: NodeIndex, dst: NodeIndex, mode: ReductionMode, - input_size: &ProblemSize, - cost_fn: &C, - ) -> Option> { - let mut costs: HashMap = HashMap::new(); - let mut sizes: HashMap = HashMap::new(); - let mut prev: HashMap = HashMap::new(); - let mut heap = BinaryHeap::new(); + initial: L, + exhaustive: bool, + ) -> Vec<(ReductionPath, L)> { + // `label` is `Option` so an evicted entry (dominated or cap-truncated) can free its + // label immediately via `take()`. Invariant: any arena index that is a current + // member of some bag has `label == Some`; only non-members may be `None`. + struct Entry { + node: NodeIndex, + label: Option, + pred: Option, + hops: usize, + } - costs.insert(src, 0.0); - sizes.insert(src, input_size.clone()); - heap.push(Reverse((OrderedFloat(0.0), src))); + let mut arena: Vec> = Vec::new(); + let mut bags: HashMap> = HashMap::new(); + let mut frontier: BinaryHeap, usize)>> = BinaryHeap::new(); + let mut adjacency: HashMap> = HashMap::new(); - while let Some(Reverse((cost, node))) = heap.pop() { - if node == dst { - let mut path = vec![dst]; - let mut current = dst; - while current != src { - let &prev_node = prev.get(¤t)?; - path.push(prev_node); - current = prev_node; - } - path.reverse(); - return Some(path); + arena.push(Entry { + node: src, + label: Some(initial.clone()), + pred: None, + hops: 0, + }); + bags.entry(src).or_default().push(0); + frontier.push(Reverse((OrderedFloat(initial.cost()), 0))); + + // Reconstruct the node-name path for an arena entry (used for deterministic + // tie-breaks). Returns the sequence of node names from source to `idx`. + let name_path = |arena: &Vec>, idx: usize| -> Vec<&'static str> { + let mut names = Vec::new(); + let mut cur = Some(idx); + while let Some(i) = cur { + names.push(self.nodes[self.graph[arena[i].node]].name); + cur = arena[i].pred; } + names.reverse(); + names + }; - if cost.0 > *costs.get(&node).unwrap_or(&f64::INFINITY) { + while let Some(Reverse((_cost, idx))) = frontier.pop() { + let node = arena[idx].node; + // Skip stale entries (removed from their bag because dominated / capped out). + if !bags.get(&node).is_some_and(|b| b.contains(&idx)) { continue; } - - let current_size = match sizes.get(&node) { - Some(s) => s.clone(), - None => continue, + // Clone the current label ONCE, up front. A live bag member always has + // `Some` (invariant above), so the `else` is unreachable. Using this local for + // every extend below means we never read `arena[idx].label` inside the edge + // loop — which also removes the self-edge hazard where extending a target == + // `node` edge could `take()` this entry's label mid-loop. + let Some(cur_label) = arena[idx].label.clone() else { + continue; }; + // The destination is terminal: keep it in the front, never expand it. + if node == dst { + continue; + } + if arena[idx].hops >= HOP_CAP { + continue; + } - for edge_ref in self.graph.edges(node) { - if !Self::edge_supports_mode(edge_ref.weight(), mode) { + // Deterministic edge order, cached because many labels can visit one node. + let edges = adjacency + .entry(node) + .or_insert_with(|| { + let mut edges: Vec<(NodeIndex, EdgeIndex)> = self + .graph + .edges(node) + .filter(|e| Self::edge_supports_mode(e.weight(), mode)) + .map(|e| (e.target(), e.id())) + .collect(); + edges.sort_by(|a, b| { + let na = &self.nodes[self.graph[a.0]]; + let nb = &self.nodes[self.graph[b.0]]; + (na.name, &na.variant).cmp(&(nb.name, &nb.variant)) + }); + edges + }) + .clone(); + + let hops = arena[idx].hops; + for (target, edge_idx) in edges { + let weight = &self.graph[edge_idx]; + let target_node = &self.nodes[self.graph[target]]; + let redge = ReductionEdge { + overhead: &weight.overhead, + reduce_fn: weight.reduce_fn, + capabilities: weight.capabilities, + target_name: target_node.name, + target_variant: &target_node.variant, + }; + let Some(new_label) = cur_label.extend(&redge) else { continue; + }; + let new_cost = new_label.cost(); + // Componentwise dominance against the target's bag. + if !exhaustive { + let bag = bags.entry(target).or_default(); + // Dominated by an existing bag member? (Bag members are always `Some`.) + if bag.iter().any(|&j| { + arena[j] + .label + .as_ref() + .is_some_and(|l| l.dominates(&new_label)) + }) { + continue; + } + // Evict every bag member the new label dominates. `Vec::retain` does + // not surface the removed elements, so collect their indices, drop them + // from the bag, then free their labels (`take()`) so nothing dominated + // lingers in the arena. + let mut evicted: Vec = Vec::new(); + bag.retain(|&j| { + let dominated = arena[j] + .label + .as_ref() + .is_some_and(|l| new_label.dominates(l)); + if dominated { + evicted.push(j); + } + !dominated + }); + for j in evicted { + arena[j].label = None; + } } - let overhead = &edge_ref.weight().overhead; - let next = edge_ref.target(); - - let edge_cost = cost_fn.edge_cost(overhead, ¤t_size); - let new_cost = cost.0 + edge_cost; - let new_size = overhead.evaluate_output_size(¤t_size); - - if new_cost < *costs.get(&next).unwrap_or(&f64::INFINITY) { - costs.insert(next, new_cost); - sizes.insert(next, new_size); - prev.insert(next, node); - heap.push(Reverse((OrderedFloat(new_cost), next))); + let nidx = arena.len(); + arena.push(Entry { + node: target, + label: Some(new_label), + pred: Some(idx), + hops: hops + 1, + }); + bags.entry(target).or_default().push(nidx); + frontier.push(Reverse((OrderedFloat(new_cost), nidx))); + + // Enforce the per-node bag cap with a deterministic tie-break. + if bags[&target].len() > BAG_CAP { + let mut entries = bags[&target].clone(); + // Bag members are always `Some`; the `unwrap_or(INFINITY)` is defensive. + let entry_cost = |i: usize| { + arena[i] + .label + .as_ref() + .map(|l| l.cost()) + .unwrap_or(f64::INFINITY) + }; + entries.sort_by(|&a, &b| { + entry_cost(a) + .partial_cmp(&entry_cost(b)) + .unwrap_or(std::cmp::Ordering::Equal) + .then_with(|| arena[a].hops.cmp(&arena[b].hops)) + .then_with(|| name_path(&arena, a).cmp(&name_path(&arena, b))) + }); + // Free the labels of the truncated tail before dropping their indices. + for &j in &entries[BAG_CAP..] { + arena[j].label = None; + } + entries.truncate(BAG_CAP); + bags.insert(target, entries); } } } - None + // The front is the (live) bag at the destination. + let mut front: Vec<(ReductionPath, L)> = bags + .get(&dst) + .map(|b| b.as_slice()) + .unwrap_or(&[]) + .iter() + .map(|&idx| { + let mut node_path = Vec::new(); + let mut cur = Some(idx); + while let Some(i) = cur { + node_path.push(arena[i].node); + cur = arena[i].pred; + } + node_path.reverse(); + ( + self.node_path_to_reduction_path(&node_path), + // Live dst bag members are always `Some` (bag-member invariant). + arena[idx] + .label + .clone() + .expect("live dst bag member has a label"), + ) + }) + .collect(); + + // Deterministic ordering of the front. + front.sort_by(|a, b| { + a.1.cost() + .partial_cmp(&b.1.cost()) + .unwrap_or(std::cmp::Ordering::Equal) + .then_with(|| a.0.len().cmp(&b.0.len())) + .then_with(|| a.0.type_names().cmp(&b.0.type_names())) + }); + front + } + + /// Name-keyed entry to [`pareto_search`](Self::pareto_search): resolves the source + /// and target variant nodes, then runs the generic search. Returns an empty vector + /// if either endpoint is not registered. Test-only: drives the generic kernel with a + /// custom label on a hand-built graph. + #[cfg(test)] + #[allow(clippy::too_many_arguments)] + pub(crate) fn pareto_search_by_name( + &self, + source: &str, + source_variant: &BTreeMap, + target: &str, + target_variant: &BTreeMap, + mode: ReductionMode, + initial: L, + exhaustive: bool, + ) -> Vec<(ReductionPath, L)> { + let (Some(src), Some(dst)) = ( + self.lookup_node(source, source_variant), + self.lookup_node(target, target_variant), + ) else { + return vec![]; + }; + self.pareto_search(src, dst, mode, initial, exhaustive) + } + + /// Pick the best element of a Pareto front under the deterministic tie-break + /// (smallest cost, then fewest hops, then lexicographic node names). The front is + /// already sorted by [`pareto_search`](Self::pareto_search), so this returns the + /// first element. + fn pick_best_front( + &self, + front: &mut Vec<(ReductionPath, L)>, + ) -> Option<(ReductionPath, L)> { + if front.is_empty() { + None + } else { + Some(front.remove(0)) + } + } + + /// Deterministic total-order key for a node-index path. + /// + /// Reproduces the `Name/val1/val2` slash signature the CLI historically used + /// as an ordering tiebreak, but computed purely from library node data so the + /// ordering lives in exactly one place. Within a fixed path length the length + /// contributes nothing, so sorting a same-length level by this key yields a + /// reproducible, build-independent order (BTreeMap variant iteration is + /// deterministic). Distinct simple paths produce distinct keys because each + /// node is a unique `(name, variant)` pair. + fn path_order_key(&self, node_path: &[NodeIndex]) -> String { + let mut key = String::new(); + for (i, &idx) in node_path.iter().enumerate() { + if i > 0 { + key.push('>'); + } + let node = &self.nodes[self.graph[idx]]; + key.push_str(node.name); + for v in node.variant.values() { + key.push('/'); + key.push_str(v); + } + } + key } /// Convert a node index path to a `ReductionPath`. @@ -559,6 +796,90 @@ impl ReductionGraph { ReductionPath { steps } } + /// Enumerate every witness-capable simple path from `src` to `dst`, executing each + /// reduction as it is reached and retaining the measured-smallest completed target. + /// + /// This is deliberately separate from [`pareto_search`](Self::pareto_search): no + /// dominance relation, hop cap, bag cap, or scalar branch-and-bound is valid for a + /// structure-dependent concrete instance. Repeated nodes are excluded because this + /// API searches graph paths (not unbounded walks); that is the sole structural + /// termination condition. + fn measured_best_simple_path<'a>( + &self, + src: NodeIndex, + dst: NodeIndex, + mode: ReductionMode, + initial: MeasuredLabel<'a>, + ) -> Option<(ReductionPath, MeasuredLabel<'a>)> { + let mut stack = vec![(src, vec![src], initial)]; + let mut adjacency: HashMap> = HashMap::new(); + let mut best: Option<(Vec, MeasuredLabel<'a>)> = None; + + while let Some((node, node_path, label)) = stack.pop() { + if node == dst { + let candidate_key = ( + label.measured_size().total(), + node_path.len(), + self.path_order_key(&node_path), + ); + let is_better = best.as_ref().is_none_or(|(best_path, best_label)| { + let best_key = ( + best_label.measured_size().total(), + best_path.len(), + self.path_order_key(best_path), + ); + candidate_key < best_key + }); + if is_better { + best = Some((node_path, label)); + } + continue; + } + + let edges = adjacency + .entry(node) + .or_insert_with(|| { + let mut edges: Vec<(NodeIndex, EdgeIndex)> = self + .graph + .edges(node) + .filter(|e| Self::edge_supports_mode(e.weight(), mode)) + .map(|e| (e.target(), e.id())) + .collect(); + edges.sort_by(|a, b| { + let na = &self.nodes[self.graph[a.0]]; + let nb = &self.nodes[self.graph[b.0]]; + (na.name, &na.variant).cmp(&(nb.name, &nb.variant)) + }); + edges + }) + .clone(); + + // Reverse push order so DFS visits the deterministic ascending edge order. + for (target, edge_idx) in edges.into_iter().rev() { + if node_path.contains(&target) { + continue; + } + let weight = &self.graph[edge_idx]; + let target_node = &self.nodes[self.graph[target]]; + let edge = ReductionEdge { + overhead: &weight.overhead, + reduce_fn: weight.reduce_fn, + capabilities: weight.capabilities, + target_name: target_node.name, + target_variant: &target_node.variant, + }; + let Some(next_label) = label.extend(&edge) else { + continue; + }; + let mut next_path = node_path.clone(); + next_path.push(target); + stack.push((target, next_path, next_label)); + } + } + + best.map(|(path, label)| (self.node_path_to_reduction_path(&path), label)) + } + /// Find all simple paths between two specific problem variants. /// /// Uses `all_simple_paths` on the variant-level graph from the exact @@ -679,18 +1000,41 @@ impl ReductionGraph { None => return vec![], }; - let paths: Vec> = all_simple_paths::< - Vec, - _, - std::hash::RandomState, - >(&self.graph, src, dst, 0, max_intermediate_nodes) - .take(limit) - .collect(); + // Enumerate every simple path in a single DFS pass and keep only the `limit` + // that sort smallest under the deterministic total order: fewest nodes first + // (shortest routes), then by `path_order_key`. Taking `limit` in petgraph's raw + // DFS discovery order (the previous approach) could drop a short route + // discovered late while returning a long route discovered early. A single + // bounded max-heap keyed by `(node count, order key)` retains exactly those + // `limit` paths — push each candidate, and once over capacity pop the current + // largest — so ordering and the truncated subset are reproducible and + // build-independent with O(limit) memory, however many paths the graph holds. + // (`limit == 0` falls out naturally: every push is immediately popped.) + let max_intermediate = + max_intermediate_nodes.unwrap_or_else(|| self.graph.node_count().saturating_sub(2)); + + let mut heap: BinaryHeap<(usize, String, Vec)> = BinaryHeap::new(); + for p in all_simple_paths::, _, std::hash::RandomState>( + &self.graph, + src, + dst, + 0, + Some(max_intermediate), + ) { + if !self.node_path_supports_mode(&p, mode) { + continue; + } + let key = self.path_order_key(&p); + heap.push((p.len(), key, p)); + if heap.len() > limit { + heap.pop(); + } + } - paths - .iter() - .filter(|p| self.node_path_supports_mode(p, mode)) - .map(|p| self.node_path_to_reduction_path(p)) + // `into_sorted_vec` yields ascending `(node count, order key)` order. + heap.into_sorted_vec() + .into_iter() + .map(|(_, _, p)| self.node_path_to_reduction_path(&p)) .collect() } @@ -968,10 +1312,14 @@ impl ReductionGraph { for entry in inventory::iter:: { if entry.source_name == name { - let result = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| { - (entry.source_size_fn)(instance) - })); - if let Ok(size) = result { + // A reduction's `source_size_fn` downcasts `instance` to its own + // source variant and panics on a mismatch; iterating every + // same-name entry means the non-matching variants panic-and-recover. + // Route through the silencer so these expected, caught panics do not + // spam stderr (the plain `catch_unwind` here did). + let result = + crate::rules::pareto::catch_reduction(|| (entry.source_size_fn)(instance)); + if let Some(size) = result { for (k, v) in size.components { if seen.insert(k.clone()) { merged.push((k, v)); @@ -1561,10 +1909,254 @@ impl ReductionGraph { } } +/// A concrete reduction path selected by the measured Pareto search. +/// +/// Holds the winning [`ReductionPath`], its **measured** final target +/// [`ProblemSize`], and the already-constructed reduction chain so downstream +/// solve/witness extraction reuses it without re-executing the reductions. +pub struct MeasuredPath { + /// The variant-level path. + pub path: ReductionPath, + /// Measured size of the final target problem. + pub size: ProblemSize, + /// The executed reduction steps (one per hop), shared via `Rc`. + steps: Vec>, +} + +impl MeasuredPath { + /// Get the final target problem as a type-erased reference. + pub fn target_problem_any(&self) -> &dyn Any { + self.steps + .last() + .expect("MeasuredPath has no steps") + .target_problem_any() + } + + /// Extract a solution from target space back to source space. + pub fn extract_solution(&self, target_solution: &[usize]) -> Vec { + self.steps + .iter() + .rev() + .fold(target_solution.to_vec(), |sol, step| { + step.extract_solution_dyn(&sol) + }) + } +} + +impl ReductionGraph { + /// Find the reduction path with the smallest **measured** final target size. + /// + /// Unlike [`find_cheapest_path_mode`](Self::find_cheapest_path_mode), which ranks + /// paths by overhead *formulas* (scaling upper bounds that can be arbitrarily loose + /// on structure-dependent constructions), this runs the [`MeasuredLabel`] domain: + /// it *actually executes* each reduction on `source_instance` and measures the real + /// constructed target size. Asymptotic overhead formulas are not treated as concrete + /// bounds and do not prune candidates. See design doc M3/F3b. + /// + /// `budget` is the hard total-size limit (sum of `ProblemSize` components); use + /// [`DEFAULT_SIZE_BUDGET`](crate::rules::DEFAULT_SIZE_BUDGET) for the default. + /// The search exhaustively enumerates witness-capable simple paths. It does not use + /// dominance pruning, branch-and-bound, or the generic Pareto kernel's bag/hop caps: + /// neither size vectors nor serialized state equality discard a route. The + /// post-construction measured-budget guard still applies. + /// Because the target must be built before it can be measured, the budget is not an + /// anti-OOM guarantee. + /// + /// Returns `None` if no in-budget witness-capable path exists (or `source == target`). + #[allow(clippy::too_many_arguments)] + pub fn find_measured_best_path( + &self, + source: &str, + source_variant: &BTreeMap, + target: &str, + target_variant: &BTreeMap, + mode: ReductionMode, + source_instance: &dyn Any, + budget: usize, + ) -> Option { + let src = self.lookup_node(source, source_variant)?; + let dst = self.lookup_node(target, target_variant)?; + if src == dst { + return None; + } + let source_size = Self::compute_source_size(source, source_instance); + let initial = MeasuredLabel::new(source_instance, source_size, budget); + let (path, label) = self.measured_best_simple_path(src, dst, mode, initial)?; + let steps: Vec> = label.chain().to_vec(); + if steps.is_empty() { + return None; + } + Some(MeasuredPath { + path, + size: label.measured_size().clone(), + steps, + }) + } + + /// Compute the **asymptotic Pareto front** of reduction paths from `source` to + /// `target` — the instance-free path search (design doc M3/F3a). + /// + /// Runs the generic [Pareto label-setting search](Self::pareto_search) with the + /// [`GrowthLabel`] domain: no concrete instance is needed, and each returned path + /// carries its composed Big-O per target size field (in the source problem's size + /// variables), read off the returned label. Because asymptotic growth over several + /// size variables is a *partial* order, the answer is a front: possibly several + /// mutually incomparable optimal paths (one better in one size field, another in a + /// different one). Paths whose composed growth is [`Growth::Unknown`] (nonlinear + /// exponent, factorial) are still returned, with those fields marked `Unknown` — + /// never a fabricated bound. + /// + /// The front reports **one representative path per distinct growth vector**: the + /// asymptotic front is a Pareto set over *growth vectors*, not routes. Many + /// syntactically different reduction chains compose to the exact same Big-O per size + /// field (e.g. dozens of `MinimumVertexCover → … → ILP` routes all yield + /// `num_constraints = O(num_edges), num_vars = O(num_vertices)`); reporting each + /// route would drown the ~1–3 genuinely distinct trade-offs the user cares about. + /// So equal-growth paths are deduplicated ([`GrowthLabel`] derives `PartialEq`), + /// keeping the deterministic best per group: fewest hops, then lexicographic + /// node-name path. Deduplication is purely by the growth vector, so two paths that + /// reach *different* target variants (e.g. `ILP/bool` vs `ILP/i32`) with the same + /// composed Big-O collapse to a single representative — the endpoint variant is not + /// part of the asymptotic identity. + /// + /// The front is ordered deterministically by (hops, lexicographic node names), so + /// the output is byte-identical across runs and platforms. Returns an empty vector + /// if either endpoint is unregistered or no path exists. + pub fn asymptotic_front( + &self, + source: &str, + source_variant: &BTreeMap, + target: &str, + target_variant: &BTreeMap, + mode: ReductionMode, + ) -> Vec<(ReductionPath, GrowthLabel)> { + let (Some(src), Some(dst)) = ( + self.lookup_node(source, source_variant), + self.lookup_node(target, target_variant), + ) else { + return vec![]; + }; + let source_fields = self.size_field_names(source); + let initial = GrowthLabel::source(&source_fields); + let mut front = self.pareto_search(src, dst, mode, initial, false); + // Order per the issue's contract: (hops, lexicographic node names). The kernel's + // own ordering leads with `cost()`, which is only a search heuristic. Sorting + // first also puts the deterministic best route of each equal-growth group ahead + // of its duplicates, so the dedup below keeps the right representative. + front.sort_by(|a, b| { + a.0.len() + .cmp(&b.0.len()) + .then_with(|| a.0.type_names().cmp(&b.0.type_names())) + }); + // Collapse to one representative per distinct growth vector. `GrowthLabel`'s + // `PartialEq` compares the field → growth map, i.e. the composed Big-O per size + // field; genuinely incomparable vectors are never equal, so they all survive. + // O(n^2), but a front is a handful of entries. + let mut deduped: Vec<(ReductionPath, GrowthLabel)> = Vec::new(); + for entry in front { + if !deduped.iter().any(|(_, label)| *label == entry.1) { + deduped.push(entry); + } + } + deduped + } + + /// Find the measured-smallest path from `source` to **any** variant of the target + /// problem name `target`. + /// + /// Runs [`find_measured_best_path`](Self::find_measured_best_path) once per target + /// variant and returns the overall measured-smallest result, with a deterministic + /// tie-break by (measured total size, hops, node-name path). + pub fn find_measured_best_path_to_name( + &self, + source: &str, + source_variant: &BTreeMap, + target: &str, + mode: ReductionMode, + source_instance: &dyn Any, + budget: usize, + ) -> Option { + let mut best: Option = None; + for tv in self.variants_for(target) { + let Some(candidate) = self.find_measured_best_path( + source, + source_variant, + target, + &tv, + mode, + source_instance, + budget, + ) else { + continue; + }; + let better = match &best { + None => true, + Some(cur) => { + let c = (candidate.size.total(), candidate.path.len()); + let b = (cur.size.total(), cur.path.len()); + c < b || (c == b && candidate.path.type_names() < cur.path.type_names()) + } + }; + if better { + best = Some(candidate); + } + } + best + } +} + +#[cfg(test)] +impl ReductionGraph { + /// Build a bare reduction graph from an explicit node/edge list (test-only). + /// + /// Nodes carry the empty variant and empty complexity; each edge carries a + /// [`ReductionEdgeData`]. This lets tests exercise the generic Pareto search on a + /// hand-built topology (e.g. the negative-control diamond) without depending on the + /// registered inventory. + pub(crate) fn from_test_edges( + node_names: &[&'static str], + edges: &[(&'static str, &'static str, ReductionEdgeData)], + ) -> Self { + let mut graph: DiGraph = DiGraph::new(); + let mut nodes: Vec = Vec::new(); + let mut name_to_nodes: HashMap<&'static str, Vec> = HashMap::new(); + let mut index_of: HashMap<&'static str, NodeIndex> = HashMap::new(); + + for &name in node_names { + let node_id = nodes.len(); + nodes.push(VariantNode { + name, + variant: BTreeMap::new(), + complexity: "", + }); + let idx = graph.add_node(node_id); + index_of.insert(name, idx); + name_to_nodes.entry(name).or_default().push(idx); + } + + for (src, dst, data) in edges { + let s = index_of[src]; + let d = index_of[dst]; + graph.add_edge(s, d, data.clone()); + } + + Self { + graph, + nodes, + name_to_nodes, + default_variants: HashMap::new(), + } + } +} + #[cfg(test)] #[path = "../unit_tests/rules/graph.rs"] mod tests; +#[cfg(test)] +#[path = "../unit_tests/rules/pareto.rs"] +mod pareto_tests; + #[cfg(test)] #[path = "../unit_tests/rules/reduction_path_parity.rs"] mod reduction_path_parity_tests; diff --git a/src/rules/ilp_i32_ilp_bool.rs b/src/rules/ilp_i32_ilp_bool.rs index 6577cb5e5..98460be44 100644 --- a/src/rules/ilp_i32_ilp_bool.rs +++ b/src/rules/ilp_i32_ilp_bool.rs @@ -264,7 +264,7 @@ impl ReductionResult for ReductionIntILPToBinaryILP { } #[reduction(overhead = { - num_vars = "31 * num_variables", + num_vars = "31 * num_vars", num_constraints = "num_constraints", })] impl ReduceTo> for ILP { diff --git a/src/rules/mod.rs b/src/rules/mod.rs index e648997a4..d66a636ee 100644 --- a/src/rules/mod.rs +++ b/src/rules/mod.rs @@ -2,6 +2,7 @@ pub mod analysis; pub mod cost; +pub mod pareto; pub mod registry; pub use cost::{ CustomCost, Minimize, MinimizeOutputSize, MinimizeSteps, MinimizeStepsThenOverhead, PathCostFn, @@ -402,9 +403,15 @@ pub(crate) mod undirectedflowlowerbounds_ilp; #[cfg(feature = "ilp-solver")] pub(crate) mod undirectedtwocommodityintegralflow_ilp; +#[cfg(test)] +pub(crate) use graph::ReductionEdgeData; pub use graph::{ - AggregateReductionChain, NeighborInfo, NeighborTree, ReductionChain, ReductionEdgeInfo, - ReductionGraph, ReductionMode, ReductionPath, ReductionStep, TraversalFlow, + AggregateReductionChain, MeasuredPath, NeighborInfo, NeighborTree, ReductionChain, + ReductionEdgeInfo, ReductionGraph, ReductionMode, ReductionPath, ReductionStep, TraversalFlow, +}; +pub use pareto::{ + CostLabel, GrowthLabel, MeasuredLabel, PathLabel, ReductionEdge, BAG_CAP, DEFAULT_SIZE_BUDGET, + HOP_CAP, }; pub use traits::{ AggregateReductionResult, ReduceTo, ReduceToAggregate, ReductionAutoCast, ReductionResult, diff --git a/src/rules/pareto.rs b/src/rules/pareto.rs new file mode 100644 index 000000000..d85a3aff8 --- /dev/null +++ b/src/rules/pareto.rs @@ -0,0 +1,446 @@ +//! Pareto label-setting search over the reduction graph. +//! +//! This module replaces the old scalar Dijkstra (`ReductionGraph::dijkstra`) with a +//! generic multi-label search. The core motivation (issue #788, design doc +//! `docs/design/symbolic-growth-domain.md`, section M3/F3b) is that edge costs are +//! **path-dependent**: the cost of a reduction depends on the size of the problem +//! accumulated along the path so far. Scalar Dijkstra keeps only the cheapest-so-far +//! label per node, so a cheaper-but-larger intermediate state can poison downstream +//! choices — it can miss the path whose *final* target is smallest. +//! +//! The fix is the standard algorithm for partial-order path costs — **multi-label +//! Pareto search** (Martins 1984; McRAPTOR-style per-node label bags). Each node keeps +//! an antichain of non-dominated labels (a "bag"); a label is only pruned when another +//! label at the same node dominates it. See [`ReductionGraph::pareto_search`]. +//! +//! Two search domains are provided: +//! - [`CostLabel`]: a scalar formula label that reproduces Dijkstra's behavior for the +//! existing `PathCostFn` cost functions (used by `find_cheapest_path*`). It carries the +//! accumulated `ProblemSize` (from overhead formulas) and an additive scalar cost. +//! - [`MeasuredLabel`]: concrete-instance state used by a separate exhaustive simple-path +//! search. It *actually executes* each reduction and measures the real constructed target +//! size. Asymptotic overhead formulas are not used as concrete budget bounds. + +use crate::expr::Expr; +use crate::growth::Growth; +use crate::rules::cost::PathCostFn; +use crate::rules::registry::{EdgeCapabilities, ReduceFn, ReductionOverhead}; +use crate::rules::traits::DynReductionResult; +use crate::types::ProblemSize; +use std::any::Any; +use std::cell::Cell; +use std::collections::{BTreeMap, BTreeSet, HashMap}; +use std::panic; +use std::rc::Rc; +use std::sync::Once; + +thread_local! { + /// When set, the installed panic hook suppresses output on the current thread. + static SILENCE_PANIC: Cell = const { Cell::new(false) }; +} + +static HOOK_INIT: Once = Once::new(); + +/// Run `f`, catching any panic and returning `None`, without printing the panic to +/// stderr on this thread. +/// +/// During the measured search we deliberately execute candidate reductions to measure +/// their real output size. A reduction whose preconditions the current instance violates +/// panics (its macro-generated dispatch downcasts and unwraps); such an edge is simply +/// not a viable path, so we treat the panic as "edge infeasible" and prune it — the +/// design's guarantee that path selection never crashes. The thread-local silencer keeps +/// this expected, recovered panic from spamming stderr while leaving genuine panics on +/// other threads untouched. +pub(crate) fn catch_reduction(f: impl FnOnce() -> R) -> Option { + HOOK_INIT.call_once(|| { + let prev = panic::take_hook(); + panic::set_hook(Box::new(move |info| { + if SILENCE_PANIC.with(|s| s.get()) { + return; + } + prev(info); + })); + }); + SILENCE_PANIC.with(|s| s.set(true)); + let result = panic::catch_unwind(panic::AssertUnwindSafe(f)); + SILENCE_PANIC.with(|s| s.set(false)); + result.ok() +} + +/// Default post-construction total-size budget for the measured search (in "size units", +/// i.e. the sum of all `ProblemSize` components). +/// +/// A reduction's target must exist before it can be measured, so this limits which +/// constructed instances remain eligible for further search; it cannot prevent the +/// construction itself from exhausting memory. +pub const DEFAULT_SIZE_BUDGET: usize = 10_000_000; + +/// Maximum number of reduction steps (hops) explored along any path. +pub const HOP_CAP: usize = 16; + +/// Maximum number of non-dominated labels retained per node. On overflow, the bag is +/// truncated by a deterministic tie-break (never by iteration order). +pub const BAG_CAP: usize = 32; + +/// A borrowed view of one reduction edge, handed to [`PathLabel::extend`]. +/// +/// It exposes exactly what a label needs to advance: the overhead formula (for symbolic +/// and formula-based labels), the executable reduction function (for measured execution), +/// the edge capabilities, and the target node's identity (for measuring the constructed +/// target's size by name). +pub struct ReductionEdge<'g> { + /// Overhead expressions mapping source size fields to target size fields. + pub overhead: &'g ReductionOverhead, + /// Type-erased witness reduction executor, if this edge supports witness/config mode. + pub reduce_fn: Option, + /// Capability metadata for the edge. + pub capabilities: EdgeCapabilities, + /// Target problem name (e.g. "ILP"). + pub target_name: &'static str, + /// Target problem variant. + pub target_variant: &'g BTreeMap, +} + +/// A path cost that composes along reduction edges under a partial order. +/// +/// **Isotonicity invariant (correctness condition for dominance pruning):** if label +/// `A` dominates label `B`, then for any edge `e`, `A.extend(e)` dominates `B.extend(e)` +/// (when both are `Some`). This follows from the monotonicity of overhead / reduction +/// size in the source size. The Pareto search relies on it to safely discard dominated +/// labels. +/// +/// The kernel prunes by [`dominates`](PathLabel::dominates) alone — it does **not** +/// branch-and-bound on [`cost`](PathLabel::cost). Dominance is exact for every label +/// domain, whereas a scalar B&B bound would only be sound for a monotone `cost`; a label's +/// scalar summary may shrink across an edge or summarize an incomparable growth vector. +/// `cost` is used only for frontier ordering and the deterministic final tie-break. +pub trait PathLabel: Clone { + /// Advance this label across `edge`. Returns `None` when a label-domain guard rejects + /// the edge. A `None` must be *isotone*: + /// if `A` dominates `B` and `A.extend(e)` is `None`, that is fine, but a guard must + /// never prune a dominating label while keeping a dominated one. + fn extend(&self, edge: &ReductionEdge) -> Option; + + /// Partial order used to keep each node's bag an antichain. Implementations must + /// satisfy the isotonicity invariant above. + fn dominates(&self, other: &Self) -> bool; + + /// Scalar summary used only for frontier ordering and the deterministic final + /// tie-break — never for pruning (the kernel prunes by [`dominates`] alone). Smaller + /// is better. It need not be monotone along `extend`. + /// + /// [`dominates`]: PathLabel::dominates + fn cost(&self) -> f64; +} + +/// Formula-based label for a [`PathCostFn`]. +/// +/// Carries the accumulated `ProblemSize` (advanced through overhead formulas) and the +/// additive scalar cost. Because a future edge's [`edge_cost`](PathCostFn::edge_cost) +/// depends on the carried size, dominance is **componentwise Pareto over `(cost, size)`**, +/// not scalar: a cheaper-but-larger prefix must not evict a costlier-but-smaller one whose +/// continuation is globally cheapest. Each node therefore keeps the antichain of +/// non-dominated `(cost, size)` labels rather than a single minimum-cost representative. +pub struct CostLabel<'c, C: PathCostFn> { + size: ProblemSize, + cost: f64, + cost_fn: &'c C, +} + +// Manual `Clone` (the derive would wrongly require `C: Clone`; `cost_fn` is a reference). +impl Clone for CostLabel<'_, C> { + fn clone(&self) -> Self { + Self { + size: self.size.clone(), + cost: self.cost, + cost_fn: self.cost_fn, + } + } +} + +impl<'c, C: PathCostFn> CostLabel<'c, C> { + /// Create the initial label at the source node. + pub fn new(input_size: ProblemSize, cost_fn: &'c C) -> Self { + Self { + size: input_size, + cost: 0.0, + cost_fn, + } + } +} + +impl PathLabel for CostLabel<'_, C> { + fn extend(&self, edge: &ReductionEdge) -> Option { + let increment = self.cost_fn.edge_cost(edge.overhead, &self.size); + let new_size = edge.overhead.evaluate_output_size(&self.size); + Some(Self { + size: new_size, + cost: self.cost + increment, + cost_fn: self.cost_fn, + }) + } + + fn dominates(&self, other: &Self) -> bool { + // Path-dependent costs: a future edge's `edge_cost` depends on the carried size, + // so `self` may only evict `other` when it is componentwise no worse in BOTH the + // accumulated cost and the carried size. Scalar `cost <= other.cost` alone would + // let a cheap-but-large prefix evict the globally optimal small one. + self.cost <= other.cost && size_le(&self.size, &other.size) + } + + fn cost(&self) -> f64 { + self.cost + } +} + +/// The current constructed position of a [`MeasuredLabel`]. +#[derive(Clone)] +enum MeasuredPos<'a> { + /// At the source node: the original, un-reduced source instance. + Source(&'a dyn Any), + /// At a reduced node: the last reduction step's result. The current problem instance + /// is `result.target_problem_any()`. + Reduced(Rc), +} + +/// The concrete-instance measured label (design doc M3/F3b). +/// +/// For a concrete source instance, the **measured** target size is authoritative. +/// Asymptotic overhead formulas are deliberately not consulted: evaluating a Big-O +/// expression at one input does not produce a certified concrete upper bound. +/// `extend` runs this stack, in order: +/// +/// 1. **Execute + measure:** run `reduce_to()`, measure the real target size; over budget +/// → `None`. +/// 2. **No comparative pruning:** measured states are enumerated by a separate exhaustive +/// simple-path search. Neither size vectors nor serialized representations discard a +/// constructed route before its downstream reductions are measured, and Pareto bag/hop +/// caps do not apply. +/// +/// **Memory.** The budget is checked only after a reduction has constructed its target, +/// so it cannot prevent a reduction itself from exhausting memory. It limits which +/// constructed instances remain eligible for further search. Exhaustive simple-path +/// enumeration can take exponential time and retain large constructed chains. +#[derive(Clone)] +pub struct MeasuredLabel<'a> { + /// Measured size of the problem instance at the current node. + size: ProblemSize, + /// The reduction steps executed so far (empty at the source). Shared via `Rc` so + /// cloning a label is cheap and never re-executes a reduction. + chain: Vec>, + /// Current constructed position. + pos: MeasuredPos<'a>, + /// Hard total-size budget. + budget: usize, +} + +impl<'a> MeasuredLabel<'a> { + /// Create the initial measured label at the source node. + /// + /// `source_size` is the measured size of `source` (typically + /// `ReductionGraph::compute_source_size`). + pub fn new(source: &'a dyn Any, source_size: ProblemSize, budget: usize) -> Self { + Self { + size: source_size, + chain: Vec::new(), + pos: MeasuredPos::Source(source), + budget, + } + } + + /// The reduction chain executed to reach this label (one entry per hop). + pub(crate) fn chain(&self) -> &[Rc] { + &self.chain + } + + /// The measured problem size at this label's node. + pub(crate) fn measured_size(&self) -> &ProblemSize { + &self.size + } + + /// Execute one reduction and retain the state only when its measured target is + /// within the post-construction budget. + pub(crate) fn extend(&self, edge: &ReductionEdge) -> Option { + // Execute the reduction and measure the real target size. Executing a + // reduction whose preconditions the current instance violates panics; such an + // edge is not a viable path, so a caught panic prunes it (returns `None`). The + // measurement (`compute_source_size`) probes every same-name size function, so + // mismatched-variant probes panic internally too — both are wrapped in one + // silenced `catch_reduction`. + let reduce_fn = edge.reduce_fn?; + let current: &dyn Any = match &self.pos { + MeasuredPos::Source(s) => *s, + MeasuredPos::Reduced(r) => r.target_problem_any(), + }; + let target_name = edge.target_name; + let (result, measured) = catch_reduction(|| { + let result: Rc = Rc::from(reduce_fn(current)); + let measured = crate::rules::ReductionGraph::compute_source_size( + target_name, + result.target_problem_any(), + ); + (result, measured) + })?; + if measured.total() > self.budget { + return None; + } + + let mut chain = self.chain.clone(); + chain.push(result.clone()); + Some(Self { + size: measured, + chain, + pos: MeasuredPos::Reduced(result), + budget: self.budget, + }) + } +} + +/// Componentwise "less-or-equal in every field" test between two sizes. +/// Missing fields are treated as `0`. +fn size_le(a: &ProblemSize, b: &ProblemSize) -> bool { + a.components + .iter() + .all(|(name, av)| *av <= b.get(name).unwrap_or(0)) +} + +/// Asymptotic, **instance-free** label domain (design doc M3/F3a). +/// +/// Each entry maps one size field of the **current** node to its +/// [`Growth`](crate::growth::Growth) expressed in the **source problem's** size +/// variables. The initial label at source `S` maps every one of `S`'s size fields +/// `f` to `Growth::from_expr(Var(f))` — "field `f` grows like itself". +/// +/// [`extend`](PathLabel::extend) composes an edge's overhead into the label: each +/// target size-field's overhead `Expr` is written over the *current* node's field +/// names, so we substitute each current field's rendered growth +/// ([`Growth::to_expr`](crate::growth::Growth::to_expr)) into it and run +/// [`Growth::from_expr`](crate::growth::Growth::from_expr) on the result. This reuses +/// the whole M1+M2 growth pipeline and needs no new growth-domain primitive. A field +/// whose growth is [`Growth::Unknown`](crate::growth::Growth::Unknown) (nonlinear +/// exponent, factorial) has no `Expr`; any target field depending on it becomes +/// `Unknown` too — the bound is never fabricated. +/// +/// [`dominates`](PathLabel::dominates) is componentwise in the **search** sense +/// (smaller growth = better): `self` dominates `other` iff for *every* field `self` +/// grows no faster than `other`, and strictly slower on at least one. Because +/// `Unknown` is the top of the growth order, a label with an `Unknown` field is +/// dominated by any fully-known label — undecidable paths rank last, the honest +/// ranking. +/// +/// **Isotonicity** (the correctness condition for the kernel's dominance pruning) +/// follows from the growth domain's monotonicity axiom: `from_expr` composed with +/// substitution into weakly-monotone overhead expressions preserves the growth +/// order, so `A ⪰ B ⇒ extend(A,e) ⪰ extend(B,e)`. +#[derive(Clone, Debug, PartialEq)] +pub struct GrowthLabel { + /// Current node's size fields → growth in the source problem's variables. + fields: BTreeMap<&'static str, Growth>, +} + +impl GrowthLabel { + /// The initial label at a source node: each size field grows like itself. + /// + /// `source_fields` is the source problem's list of size-field names (e.g. from + /// [`ReductionGraph::size_field_names`](crate::rules::ReductionGraph::size_field_names)). + pub fn source(source_fields: &[&'static str]) -> Self { + let fields = source_fields + .iter() + .map(|&f| (f, Growth::from_expr(&Expr::Var(f)))) + .collect(); + GrowthLabel { fields } + } + + /// Construct directly from a field → growth map (test/introspection helper). + pub fn from_fields(fields: BTreeMap<&'static str, Growth>) -> Self { + GrowthLabel { fields } + } + + /// The current node's size fields mapped to their growth in source variables. + pub fn fields(&self) -> &BTreeMap<&'static str, Growth> { + &self.fields + } +} + +impl PathLabel for GrowthLabel { + fn extend(&self, edge: &ReductionEdge) -> Option { + // Render each current field's growth back to a display `Expr` in the source + // variables. `Unknown` growth has no `Expr` (`None`) and taints any target + // field that references it. + let rendered: BTreeMap<&'static str, Option> = + self.fields.iter().map(|(k, g)| (*k, g.to_expr())).collect(); + + // Substitution map from current field name to its rendered growth `Expr` (in + // source variables). Depends only on `rendered`, so build it once for all edges' + // output fields rather than per target field. Only present-and-known fields are + // mapped. Unlike `ReductionOverhead::compose`, an overhead variable ABSENT from + // this map is NOT a passthrough source variable: in the asymptotic label it is an + // intermediate-only field with no source-variable growth, so any target field that + // references it must be tainted (see below) rather than leaked verbatim. + let mapping: HashMap<&str, &Expr> = rendered + .iter() + .filter_map(|(k, opt)| opt.as_ref().map(|e| (*k, e))) + .collect(); + + let mut new_fields: BTreeMap<&'static str, Growth> = BTreeMap::new(); + for (target_field, expr) in &edge.overhead.output_size { + // Taint the target field if this overhead references any variable we cannot + // express in the source's variables: either a present-but-`Unknown` current + // field, or a variable absent from the label entirely (an intermediate-only + // field that would otherwise leak through `substitute` as a fake source + // variable). Both cases are exactly "not in `mapping`". + let taints = expr.variables().iter().any(|v| !mapping.contains_key(v)); + if taints { + new_fields.insert(target_field, Growth::Unknown); + continue; + } + // Substitute rendered growths into the overhead, then reduce in the growth + // domain. + let substituted = expr.substitute(&mapping); + new_fields.insert(target_field, Growth::from_expr(&substituted)); + } + // Asymptotic mode has no budget, so `extend` never prunes. + Some(GrowthLabel { fields: new_fields }) + } + + fn dominates(&self, other: &Self) -> bool { + // Search-sense componentwise dominance over the union of fields (labels + // compared are at the same node, so their field sets coincide; the union is + // defensive). `self` dominates `other` iff `self` grows no faster on every + // field and strictly slower on at least one. + // + // `Growth::dominates(a, b)` means "a grows ≥ b", with `Unknown` as top. So: + // self ≤ other on field f ⟺ other_f.dominates(self_f) + // and self is strictly better on f iff additionally NOT self_f.dominates(other_f). + let o1 = Growth::Terms(Vec::new()); // O(1): the bottom, for absent fields. + let keys: BTreeSet<&'static str> = self + .fields + .keys() + .chain(other.fields.keys()) + .copied() + .collect(); + let mut strict = false; + for k in keys { + let s = self.fields.get(k).unwrap_or(&o1); + let o = other.fields.get(k).unwrap_or(&o1); + if !o.dominates(s) { + // self grows strictly faster than other here → self does not dominate. + return false; + } + if !s.dominates(o) { + // other ≥ self but self ⋡ other ⇒ self strictly slower on this field. + strict = true; + } + } + strict + } + + fn cost(&self) -> f64 { + // Heuristic scalar summary for frontier ordering and the deterministic final + // tie-break ONLY — never for pruning (dominance is the exact partial order, and + // asymptotic growth is incomparable so no scalar bound could separate front + // members). Summed field magnitudes; `Unknown` fields dominate the sum, ranking + // undecidable paths last. + self.fields.values().map(|g| g.magnitude()).sum() + } +} diff --git a/src/solvers/ilp/solver.rs b/src/solvers/ilp/solver.rs index 51b2a0df2..08ebbacb6 100644 --- a/src/solvers/ilp/solver.rs +++ b/src/solvers/ilp/solver.rs @@ -240,48 +240,48 @@ impl ILPSolver { any.is::>() || any.is::>() || any.is::() } - /// Two-level path selection: - /// 1. Dijkstra finds the cheapest path to each ILP variant using - /// `MinimizeStepsThenOverhead` (additive edge costs: step count + log overhead). - /// 2. Across ILP variants, we pick the path whose composed final output size - /// is smallest — this is the actual ILP problem size the solver will face. - fn best_path_to_ilp( + /// Execute the first constructible preferred witness path to an ILP variant. + /// + /// Solving only requires a valid formulation; it does not require proving which of + /// every possible multi-hop formulation is concretely smallest. One shortest path is + /// considered per ILP variant, ordered deterministically by hops and node names. + fn preferred_chain_to_ilp( &self, graph: &crate::rules::ReductionGraph, name: &str, variant: &std::collections::BTreeMap, - mode: ReductionMode, instance: &dyn std::any::Any, - ) -> Option { - let ilp_variants = graph.variants_for("ILP"); + ) -> Option { let input_size = crate::rules::ReductionGraph::compute_source_size(name, instance); - let mut best_path: Option = None; - let mut best_cost = f64::INFINITY; - - for dv in &ilp_variants { - if let Some(path) = graph.find_cheapest_path_mode( - name, - variant, - "ILP", - dv, - mode, - &input_size, - &crate::rules::MinimizeStepsThenOverhead, - ) { - // Use composed final output size for cross-variant comparison, - // since this determines the actual ILP problem size. - let final_size = graph - .evaluate_path_overhead(&path, &input_size) - .unwrap_or_default(); - let cost = final_size.total() as f64; - if cost < best_cost { - best_cost = cost; - best_path = Some(path); - } + let mut candidates: Vec<_> = graph + .variants_for("ILP") + .into_iter() + .filter_map(|target_variant| { + graph.find_cheapest_path_mode( + name, + variant, + "ILP", + &target_variant, + ReductionMode::Witness, + &input_size, + &crate::rules::MinimizeSteps, + ) + }) + .collect(); + candidates.sort_by(|a, b| { + a.len() + .cmp(&b.len()) + .then_with(|| a.type_names().cmp(&b.type_names())) + }); + for path in candidates { + if let Some(chain) = + crate::rules::pareto::catch_reduction(|| graph.reduce_along_path(&path, instance)) + .flatten() + { + return Some(chain); } } - - best_path + None } pub fn try_solve_via_reduction( @@ -300,41 +300,76 @@ impl ILPSolver { let graph = crate::rules::ReductionGraph::new(); - let Some(path) = - self.best_path_to_ilp(&graph, name, variant, ReductionMode::Witness, instance) - else { - if self - .best_path_to_ilp(&graph, name, variant, ReductionMode::Aggregate, instance) - .is_some() - { - return Err(SolveViaReductionError::WitnessPathRequired { + if let Some(chain) = self.preferred_chain_to_ilp(&graph, name, variant, instance) { + let ilp_solution = self.solve_dyn(chain.target_problem_any()).ok_or_else(|| { + SolveViaReductionError::NoSolution { name: name.to_string(), - }); - } + } + })?; + return Ok(chain.extract_solution(&ilp_solution)); + } + + // A preferred shortest path can be instance-infeasible even when another route + // works. Fall back to the uncapped, execution-aware measured enumeration before + // reporting that no witness path exists. + if let Some(measured) = graph.find_measured_best_path_to_name( + name, + variant, + "ILP", + ReductionMode::Witness, + instance, + crate::rules::DEFAULT_SIZE_BUDGET, + ) { + let ilp_solution = self + .solve_dyn(measured.target_problem_any()) + .ok_or_else(|| SolveViaReductionError::NoSolution { + name: name.to_string(), + })?; + return Ok(measured.extract_solution(&ilp_solution)); + } - return Err(SolveViaReductionError::NoReductionPath { + if self.has_aggregate_path_to_ilp(&graph, name, variant) { + return Err(SolveViaReductionError::WitnessPathRequired { name: name.to_string(), }); - }; + } - let chain = graph.reduce_along_path(&path, instance).ok_or_else(|| { - SolveViaReductionError::WitnessPathRequired { - name: name.to_string(), - } - })?; - let ilp_solution = self.solve_dyn(chain.target_problem_any()).ok_or_else(|| { - SolveViaReductionError::NoSolution { - name: name.to_string(), - } - })?; - Ok(chain.extract_solution(&ilp_solution)) + Err(SolveViaReductionError::NoReductionPath { + name: name.to_string(), + }) + } + + /// Whether an aggregate-capable (but possibly not witness-capable) reduction path to + /// some ILP variant exists. Used only to distinguish "no path at all" from "a path + /// exists but cannot recover a witness" for error reporting. + fn has_aggregate_path_to_ilp( + &self, + graph: &crate::rules::ReductionGraph, + name: &str, + variant: &std::collections::BTreeMap, + ) -> bool { + let input_size = crate::types::ProblemSize::new(vec![]); + graph.variants_for("ILP").iter().any(|dv| { + graph + .find_cheapest_path_mode( + name, + variant, + "ILP", + dv, + ReductionMode::Aggregate, + &input_size, + &crate::rules::MinimizeSteps, + ) + .is_some() + }) } /// Solve a type-erased problem by finding a reduction path to ILP. /// - /// Tries all ILP variants, picks the cheapest path, reduces, solves, - /// and extracts the solution back. Falls back to direct ILP solve if - /// the problem is already an ILP type. + /// Prefers a shortest witness path to an ILP variant, reduces, solves, and extracts + /// the solution back. If the preferred constructions are instance-infeasible, it + /// falls back to exhaustive measured simple-path search. Problems already represented + /// as ILP are solved directly. /// /// Returns `None` if no path to ILP exists or the solver finds no solution. pub fn solve_via_reduction( diff --git a/src/unit_tests/big_o.rs b/src/unit_tests/big_o.rs index 6dab26625..9ed1cefc8 100644 --- a/src/unit_tests/big_o.rs +++ b/src/unit_tests/big_o.rs @@ -109,9 +109,12 @@ fn test_big_o_rejects_division() { } #[test] -fn test_big_o_rejects_negative_dominant_term() { +fn test_big_o_drops_negative_constant_factor() { + // The growth domain drops constant multipliers, sign included, so `-1 * n` + // widens to `n` (an upper bound on its magnitude) instead of being rejected. let e = Expr::Const(-1.0) * Expr::Var("n"); - assert!(big_o_normal_form(&e).is_err()); + let result = big_o_normal_form(&e).unwrap(); + assert_eq!(result.to_string(), "n"); } #[test] @@ -219,11 +222,18 @@ fn test_big_o_multivar_exp_dominates_poly() { } #[test] -fn test_big_o_pathological_nesting_errors_instead_of_hanging() { - // Regression for issue #1069: a deeply-nested power that expands - // exponentially must return an error promptly (so callers like `big_o_of` - // fall back to the un-expanded expression) rather than OOM/hang. +fn test_big_o_pathological_nesting_returns_bound_instantly() { + // Regression for issue #1069: a deeply-nested power that the old expansion + // pipeline could not normalize (it OOM'd, then refused via the term cap). + // The growth domain answers it bottom-up: `((a+b+c+d)^4)^4` raises each + // variable term to degree 16, so it returns a real bound, instantly. let sum = Expr::Var("a") + Expr::Var("b") + Expr::Var("c") + Expr::Var("d"); let e = Expr::pow(Expr::pow(sum, Expr::Const(4.0)), Expr::Const(4.0)); - assert!(big_o_normal_form(&e).is_err()); + let start = std::time::Instant::now(); + let result = big_o_normal_form(&e).unwrap(); + assert!(start.elapsed().as_millis() < 50, "should be instant"); + let s = result.to_string(); + for v in ["a^16", "b^16", "c^16", "d^16"] { + assert!(s.contains(v), "expected {v} in {s}"); + } } diff --git a/src/unit_tests/canonical.rs b/src/unit_tests/canonical.rs deleted file mode 100644 index dcf3f8fd0..000000000 --- a/src/unit_tests/canonical.rs +++ /dev/null @@ -1,165 +0,0 @@ -use super::*; -use crate::expr::Expr; - -#[test] -fn test_canonical_identity() { - let e = Expr::Var("n"); - let c = canonical_form(&e).unwrap(); - assert_eq!(c.to_string(), "n"); -} - -#[test] -fn test_canonical_add_like_terms() { - // n + n → 2 * n - let e = Expr::Var("n") + Expr::Var("n"); - let c = canonical_form(&e).unwrap(); - assert_eq!(c.to_string(), "2 * n"); -} - -#[test] -fn test_canonical_subtract_to_zero() { - // n - n → 0 - let e = Expr::Var("n") - Expr::Var("n"); - let c = canonical_form(&e).unwrap(); - assert_eq!(c.to_string(), "0"); -} - -#[test] -fn test_canonical_mixed_addition() { - // n + n - m + 2*m → 2*n + m - let e = Expr::Var("n") + Expr::Var("n") - Expr::Var("m") + Expr::Const(2.0) * Expr::Var("m"); - let c = canonical_form(&e).unwrap(); - assert_eq!(c.to_string(), "m + 2 * n"); -} - -#[test] -fn test_canonical_exp_product_identity() { - // exp(n) * exp(m) -> exp(m + n) (transcendental identity, alphabetical order) - let e = Expr::Exp(Box::new(Expr::Var("n"))) * Expr::Exp(Box::new(Expr::Var("m"))); - let c = canonical_form(&e).unwrap(); - // Verify numerical equivalence - let size = crate::types::ProblemSize::new(vec![("n", 2), ("m", 3)]); - assert!((c.eval(&size) - (2.0_f64.exp() * 3.0_f64.exp())).abs() < 1e-6); -} - -#[test] -fn test_canonical_constant_base_exp_identity() { - // 2^n * 2^m -> 2^(m + n) - let e = - Expr::pow(Expr::Const(2.0), Expr::Var("n")) * Expr::pow(Expr::Const(2.0), Expr::Var("m")); - let c = canonical_form(&e).unwrap(); - let size = crate::types::ProblemSize::new(vec![("n", 3), ("m", 4)]); - assert!((c.eval(&size) - 2.0_f64.powf(7.0)).abs() < 1e-6); -} - -#[test] -fn test_canonical_polynomial_expansion() { - // (n + m)^2 = n^2 + 2*n*m + m^2 - let e = Expr::pow(Expr::Var("n") + Expr::Var("m"), Expr::Const(2.0)); - let c = canonical_form(&e).unwrap(); - let size = crate::types::ProblemSize::new(vec![("n", 3), ("m", 4)]); - assert_eq!(c.eval(&size), 49.0); // (3+4)^2 = 49 -} - -#[test] -fn test_canonical_signed_polynomial() { - // n^3 - n^2 + 2*n + 4*n*m — should remain exact - let e = Expr::pow(Expr::Var("n"), Expr::Const(3.0)) - - Expr::pow(Expr::Var("n"), Expr::Const(2.0)) - + Expr::Const(2.0) * Expr::Var("n") - + Expr::Const(4.0) * Expr::Var("n") * Expr::Var("m"); - let c = canonical_form(&e).unwrap(); - let size = crate::types::ProblemSize::new(vec![("n", 3), ("m", 2)]); - // 27 - 9 + 6 + 24 = 48 - assert_eq!(c.eval(&size), 48.0); -} - -#[test] -fn test_canonical_division_becomes_negative_exponent() { - // n / m should canonicalize; the division is represented as m^(-1) - // which becomes an opaque factor (negative exponent) - let e = Expr::Var("n") / Expr::Var("m"); - let c = canonical_form(&e).unwrap(); - let size = crate::types::ProblemSize::new(vec![("n", 6), ("m", 3)]); - assert!((c.eval(&size) - 2.0).abs() < 1e-10); -} - -#[test] -fn test_canonical_distinct_fractional_exponents_do_not_merge() { - let e = Expr::pow(Expr::Var("n"), Expr::Const(1.0004)) - Expr::Var("n"); - let c = canonical_form(&e).unwrap(); - assert_ne!(c.to_string(), "0"); - let size = crate::types::ProblemSize::new(vec![("n", 2)]); - assert_ne!(c.eval(&size), 0.0); -} - -#[test] -fn test_canonical_constant_base_one_folds_to_constant() { - let e = Expr::pow(Expr::Const(1.0), Expr::Var("n")); - let c = canonical_form(&e).unwrap(); - assert_eq!(c.to_string(), "1"); -} - -#[test] -fn test_canonical_negative_constant_base_with_symbolic_exponent_is_rejected() { - let e = Expr::pow(Expr::Const(-2.0), Expr::Var("n")); - let err = canonical_form(&e).unwrap_err(); - assert!(matches!(err, CanonicalizationError::Unsupported(_))); -} - -#[test] -fn test_canonical_zero_constant_base_with_symbolic_exponent_is_rejected() { - let e = Expr::pow(Expr::Const(0.0), Expr::Var("n")); - let err = canonical_form(&e).unwrap_err(); - assert!(matches!(err, CanonicalizationError::Unsupported(_))); -} - -#[test] -fn test_canonical_deterministic_order() { - // m + n and n + m should produce the same canonical form - let a = canonical_form(&(Expr::Var("m") + Expr::Var("n"))).unwrap(); - let b = canonical_form(&(Expr::Var("n") + Expr::Var("m"))).unwrap(); - assert_eq!(a.to_string(), b.to_string()); -} - -#[test] -fn test_canonical_constant_folding() { - // 2 + 3 → 5 - let e = Expr::Const(2.0) + Expr::Const(3.0); - let c = canonical_form(&e).unwrap(); - assert_eq!(c.to_string(), "5"); -} - -#[test] -fn test_canonical_sqrt_as_power() { - // sqrt(n) should canonicalize the same as n^0.5 - let a = canonical_form(&Expr::Sqrt(Box::new(Expr::Var("n")))).unwrap(); - let b = canonical_form(&Expr::pow(Expr::Var("n"), Expr::Const(0.5))).unwrap(); - assert_eq!(a.to_string(), b.to_string()); -} - -#[test] -fn test_canonical_nested_power_blowup_is_capped() { - // Regression for issue #1069: a "square of a square of a sum" structure — - // the shape composed-path overheads take when they traverse - // quadratic-overhead reductions — expands exponentially. Before the cap - // this OOM'd / hung indefinitely; now it must fail fast with Unsupported - // rather than try to materialize the blown-up monomial expansion. - let sum = Expr::Var("a") + Expr::Var("b") + Expr::Var("c") + Expr::Var("d"); - // ((a+b+c+d)^4)^4 expands to >50_000 intermediate terms. - let e = Expr::pow(Expr::pow(sum, Expr::Const(4.0)), Expr::Const(4.0)); - let err = canonical_form(&e).unwrap_err(); - assert!(matches!(err, CanonicalizationError::Unsupported(_))); -} - -#[test] -fn test_canonical_moderate_power_still_expands() { - // The cap must not perturb legitimate, modestly-sized expressions: - // (a+b)^3 stays well under the cap and expands normally. - let e = Expr::pow(Expr::Var("a") + Expr::Var("b"), Expr::Const(3.0)); - let c = canonical_form(&e).unwrap(); - // a^3 + 3 a^2 b + 3 a b^2 + b^3 — compare against the same expansion - // written out flat (both go through canonical_form for identical ordering). - let expected = canonical_form(&Expr::parse("a^3 + 3*a^2*b + 3*a*b^2 + b^3")).unwrap(); - assert_eq!(c.to_string(), expected.to_string()); -} diff --git a/src/unit_tests/expr.rs b/src/unit_tests/expr.rs index 037f39c8c..fd1e217aa 100644 --- a/src/unit_tests/expr.rs +++ b/src/unit_tests/expr.rs @@ -168,95 +168,6 @@ fn test_expr_display_pow_with_complex_exponent() { assert_eq!(format!("{expr}"), "2^(m + n)"); } -#[test] -fn test_asymptotic_normal_form_drops_constant_factors() { - let expr = Expr::parse("3 * num_variables^2"); - let normalized = asymptotic_normal_form(&expr).unwrap(); - assert_eq!(normalized.to_string(), "num_variables^2"); -} - -#[test] -fn test_asymptotic_normal_form_drops_additive_constants() { - let expr = Expr::parse("num_variables + 1"); - let normalized = asymptotic_normal_form(&expr).unwrap(); - assert_eq!(normalized.to_string(), "num_variables"); -} - -#[test] -fn test_asymptotic_normal_form_canonicalizes_commutative_sum() { - let a = asymptotic_normal_form(&Expr::parse("n + m")).unwrap(); - let b = asymptotic_normal_form(&Expr::parse("m + n")).unwrap(); - assert_eq!(a, b); - assert_eq!(a.to_string(), "m + n"); -} - -#[test] -fn test_asymptotic_normal_form_canonicalizes_commutative_product() { - let a = asymptotic_normal_form(&Expr::parse("n * m")).unwrap(); - let b = asymptotic_normal_form(&Expr::parse("m * n")).unwrap(); - assert_eq!(a, b); - assert_eq!(a.to_string(), "m * n"); -} - -#[test] -fn test_asymptotic_normal_form_combines_repeated_factors() { - let normalized = asymptotic_normal_form(&Expr::parse("n * n^(1/2)")).unwrap(); - assert_eq!(normalized.to_string(), "n^1.5"); -} - -#[test] -fn test_asymptotic_normal_form_canonicalizes_exponential_product() { - let a = asymptotic_normal_form(&Expr::parse("exp(n) * exp(m)")).unwrap(); - let b = asymptotic_normal_form(&Expr::parse("exp(n + m)")).unwrap(); - assert_eq!(a, b); - assert_eq!(a.to_string(), "exp(m + n)"); -} - -#[test] -fn test_asymptotic_normal_form_canonicalizes_constant_base_exponential_product() { - let a = asymptotic_normal_form(&Expr::parse("2^n * 2^m")).unwrap(); - let b = asymptotic_normal_form(&Expr::parse("2^(n + m)")).unwrap(); - assert_eq!(a, b); - assert_eq!(a.to_string(), "2^(m + n)"); -} - -#[test] -fn test_asymptotic_normal_form_sqrt_matches_fractional_power() { - let a = asymptotic_normal_form(&Expr::parse("sqrt(n * m)")).unwrap(); - let b = asymptotic_normal_form(&Expr::parse("(n * m)^(1/2)")).unwrap(); - assert_eq!(a, b); -} - -#[test] -fn test_asymptotic_normal_form_log_of_power() { - // log(n^2) = 2*log(n) — the new engine keeps log(n^2) which is O(log(n)) - let normalized = asymptotic_normal_form(&Expr::parse("log(n^2)")).unwrap(); - // Both log(n^2) and log(n) are asymptotically equivalent - let s = normalized.to_string(); - assert!(s.contains("log"), "expected log in result, got: {s}"); - assert!(s.contains("n"), "expected n in result, got: {s}"); -} - -#[test] -fn test_asymptotic_normal_form_substitution_is_closed() { - let notation = asymptotic_normal_form(&Expr::parse("n * m")).unwrap(); - let k = Expr::parse("k"); - let k_squared = Expr::parse("k^2"); - let mapping = HashMap::from([("n", &k), ("m", &k_squared)]); - let substituted = asymptotic_normal_form(¬ation.substitute(&mapping)).unwrap(); - assert_eq!(substituted.to_string(), "k^3"); -} - -#[test] -fn test_asymptotic_normal_form_handles_subtraction() { - // n - m: the -m term survives as a negative dominant term → unsupported - assert!(asymptotic_normal_form(&Expr::parse("n - m")).is_err()); - - // n^2 - n: -n is dominated by n^2 and eliminated → works - let result = asymptotic_normal_form(&Expr::parse("n^2 - n")).unwrap(); - assert_eq!(result.to_string(), "n^2"); -} - #[test] fn test_expr_display_fractional_constant() { assert_eq!(format!("{}", Expr::Const(2.75)), "2.75"); diff --git a/src/unit_tests/growth.rs b/src/unit_tests/growth.rs new file mode 100644 index 000000000..f1d7b99dd --- /dev/null +++ b/src/unit_tests/growth.rs @@ -0,0 +1,1104 @@ +//! Unit tests for the symbolic growth domain (`src/growth.rs`). + +use super::{ + add, componentwise_max, make_growth, mul, ExpBase, ExpFactor, ExpProduct, Growth, GrowthTerm, +}; +use crate::expr::Expr; +use std::cmp::Ordering; + +/// Build a term from `(exp, poly, logs)` entry lists. +fn term( + exp: &[(&'static str, f64)], + poly: &[(&'static str, f64)], + logs: &[(&'static str, u32)], +) -> GrowthTerm { + GrowthTerm { + exp: exp + .iter() + .map(|(variable, rate)| { + ( + *variable, + ExpProduct::single(ExpBase::Constant(Expr::Const(2.0)), *rate), + ) + }) + .collect(), + poly: poly.iter().copied().collect(), + logs: logs.iter().copied().collect(), + } +} + +fn terms_of(g: &Growth) -> &[GrowthTerm] { + match g { + Growth::Terms(t) => t, + Growth::Unknown => panic!("expected Terms, got Unknown"), + } +} + +fn g(s: &str) -> Growth { + Growth::from_expr(&Expr::parse(s)) +} + +fn exp_product(factors: &[(f64, f64)]) -> ExpProduct { + ExpProduct::new( + factors + .iter() + .map(|(base, coefficient)| ExpFactor { + base: ExpBase::Constant(Expr::Const(*base)), + coefficient: *coefficient, + }) + .collect(), + ) +} + +// --- The six named verification cases from issue #1075 --- + +/// 1. No-expansion regression: the nested sum-of-squares shape that OOM'd in +/// issue #1069 is handled without expansion, quickly, with few terms. +#[test] +fn test_growth_no_expansion_regression() { + let e = Expr::parse("(12*(n + 3*m) + 5)^2 * (12*(n + 3*m) + 5)^2"); + let start = std::time::Instant::now(); + let result = Growth::from_expr(&e); + let elapsed = start.elapsed(); + + let ts = terms_of(&result); + assert!( + ts.contains(&term(&[], &[("n", 4.0)], &[])), + "expected n^4 in {ts:?}" + ); + assert!( + ts.contains(&term(&[], &[("m", 4.0)], &[])), + "expected m^4 in {ts:?}" + ); + assert!(ts.len() <= 6, "expected <= 6 terms, got {}", ts.len()); + assert!(elapsed.as_millis() < 10, "from_expr took {elapsed:?}"); +} + +/// 2. Dominance beats the old sampling heuristic: `1.001^n` dominates `n^100` +/// (any positive exponential rate outranks any polynomial degree). +#[test] +fn test_growth_exponential_dominates_polynomial() { + let exp = g("1.001^n"); + let poly = g("n^100"); + assert!(exp.dominates(&poly)); + assert!(!poly.dominates(&exp)); +} + +/// 3. Incomparability is honest: neither `n^2` nor `n*m` dominates the other, +/// and both are kept in the sum. +#[test] +fn test_growth_incomparable_terms_both_kept() { + let n2 = g("n^2"); + let nm = g("n*m"); + assert!(!n2.dominates(&nm)); + assert!(!nm.dominates(&n2)); + + let sum = g("n^2 + n*m"); + assert_eq!(terms_of(&sum).len(), 2); +} + +/// 4. Exponent rates are exact: `2^(2n)` dominates `2^n` (not conversely), and +/// `3^n` dominates `2^n` via direct symbolic base comparison. +#[test] +fn test_growth_exponent_rates_exact() { + let two_2n = g("2^(2*n)"); + let two_n = g("2^n"); + assert!(two_2n.dominates(&two_n)); + assert!(!two_n.dominates(&two_2n)); + + let three_n = g("3^n"); + assert!(three_n.dominates(&two_n)); + assert!(!two_n.dominates(&three_n)); + + let exp_2n = g("exp(2*n)"); + let exp_n = g("exp(n)"); + assert!(exp_2n.dominates(&exp_n)); + assert!(!exp_n.dominates(&exp_2n)); + + assert!(g("0.5^(-2*n)").dominates(&g("0.5^(-n)"))); + assert!(g("0.25^(-n)").dominates(&g("0.5^(-n)"))); +} + +/// Multi-base products remain incomparable when the conservative symbolic +/// rules cannot prove an ordering, even when a stronger algebra system could. +#[test] +fn test_growth_unproved_multi_base_comparison_is_retained() { + let left = g("2^(2*n) * 3^n"); + let right = g("2^n * 4^n"); + assert!(!left.dominates(&right)); + assert!(!right.dominates(&left)); + assert_eq!(terms_of(&g("2^(2*n) * 3^n + 2^n * 4^n")).len(), 2); +} + +#[test] +fn test_exponential_product_proof_rules() { + let empty = ExpProduct::empty(); + let two = exp_product(&[(2.0, 1.0)]); + let two_squared = exp_product(&[(2.0, 2.0)]); + let three = exp_product(&[(3.0, 1.0)]); + + assert_eq!(empty.cmp_proven(&empty), Some(Ordering::Equal)); + assert_eq!(empty.cmp_proven(&two), Some(Ordering::Less)); + assert_eq!(two.cmp_proven(&empty), Some(Ordering::Greater)); + assert_eq!(two_squared.cmp_proven(&two), Some(Ordering::Greater)); + assert_eq!(two.cmp_proven(&two_squared), Some(Ordering::Less)); + assert_eq!(three.cmp_proven(&two), Some(Ordering::Greater)); + assert_eq!(two.cmp_proven(&three), Some(Ordering::Less)); + + assert_eq!( + exp_product(&[(3.0, 2.0)]).cmp_proven(&exp_product(&[(2.0, 1.0)])), + Some(Ordering::Greater) + ); + assert_eq!( + exp_product(&[(2.0, 1.0)]).cmp_proven(&exp_product(&[(3.0, 2.0)])), + Some(Ordering::Less) + ); + assert_eq!( + exp_product(&[(2.0, 3.0)]).cmp_proven(&exp_product(&[(3.0, 1.0)])), + None + ); + + assert_eq!( + exp_product(&[(0.25, -1.0)]).cmp_proven(&exp_product(&[(0.5, -1.0)])), + Some(Ordering::Greater) + ); + assert_eq!( + exp_product(&[(0.5, -1.0)]).cmp_proven(&exp_product(&[(0.25, -1.0)])), + Some(Ordering::Less) + ); + assert_eq!( + exp_product(&[(0.25, -2.0)]).cmp_proven(&exp_product(&[(0.5, -1.0)])), + Some(Ordering::Greater) + ); + assert_eq!( + exp_product(&[(0.5, -1.0)]).cmp_proven(&exp_product(&[(0.25, -2.0)])), + Some(Ordering::Less) + ); + assert_eq!( + exp_product(&[(0.25, -1.0)]).cmp_proven(&exp_product(&[(0.5, -2.0)])), + None + ); + assert_eq!(two.cmp_proven(&exp_product(&[(0.5, -1.0)])), None); + + let natural = ExpProduct::single(ExpBase::Natural, 1.0); + assert_eq!(natural.cmp_proven(&two), Some(Ordering::Greater)); + assert_eq!(two.cmp_proven(&natural), Some(Ordering::Less)); + + // Arbitrary constant subtrees are preserved but compared structurally only. + let composite = ExpProduct::single(ExpBase::Constant(Expr::parse("1 + 2")), 1.0); + assert_eq!(composite.cmp_proven(&three), None); + + // Two residual products with no factorwise proof remain incomparable. + assert_eq!( + exp_product(&[(2.0, 2.0), (3.0, 1.0)]).cmp_proven(&exp_product(&[(2.0, 1.0), (4.0, 1.0)])), + None + ); +} + +#[test] +fn test_exponential_product_canonicalization() { + let combined = ExpProduct::new(vec![ + ExpFactor { + base: ExpBase::Constant(Expr::Const(2.0)), + coefficient: 1.0, + }, + ExpFactor { + base: ExpBase::Constant(Expr::Const(2.0)), + coefficient: 2.0, + }, + ExpFactor { + base: ExpBase::Constant(Expr::Const(3.0)), + coefficient: 0.0, + }, + ]); + assert_eq!(combined, exp_product(&[(2.0, 3.0)])); + + let cancelled = ExpProduct::new(vec![ + ExpFactor { + base: ExpBase::Constant(Expr::Const(2.0)), + coefficient: 1.0, + }, + ExpFactor { + base: ExpBase::Constant(Expr::Const(2.0)), + coefficient: -1.0, + }, + ]); + assert!(cancelled.is_empty()); +} + +#[test] +fn test_growth_multi_base_product_is_deterministic() { + let left = g("2^n * 3^n"); + let right = g("3^n * 2^n"); + assert_eq!(left, right); + assert_eq!( + serde_json::to_string(&left).unwrap(), + serde_json::to_string(&right).unwrap() + ); +} + +#[test] +fn test_proven_equal_exponential_spelling_is_deterministic() { + let natural_first = g("exp(n) + 2.718281828459045^n"); + let literal_first = g("2.718281828459045^n + exp(n)"); + assert_eq!(natural_first, literal_first); + assert_eq!(natural_first.to_big_o(), literal_first.to_big_o()); +} + +/// 5. Widening: subtraction widens to addition, including the `sqrt((a-b)^2)` +/// absolute-value idiom. +#[test] +fn test_growth_widening() { + assert_eq!(g("n - m"), g("n + m")); + assert_eq!(g("sqrt((n - m)^2)"), g("n + m")); +} + +/// 6. Determinism: the antichain is canonically sorted, so structurally +/// equivalent inputs are equal regardless of term order. +#[test] +fn test_growth_determinism() { + assert_eq!(g("n*m + m*n"), g("m*n + n*m")); +} + +// --- Negative control --- + +/// Unsupported content widens to `Unknown`, and `Unknown` absorbs through add +/// and mul — unsupported content can never silently produce a fake bound. +#[test] +fn test_growth_unknown_negative_control() { + assert_eq!(g("2^(n*k)"), Growth::Unknown); + assert_eq!(g("factorial(n)"), Growth::Unknown); + + // Absorption through the real `from_expr` add/mul paths. + assert_eq!(g("factorial(n) + n^2"), Growth::Unknown); + assert_eq!(g("n^2 + factorial(n)"), Growth::Unknown); + assert_eq!(g("factorial(n) * n^2"), Growth::Unknown); + assert_eq!(g("n^2 * factorial(n)"), Growth::Unknown); + + // Absorption at the operation level too. + let n2 = g("n^2"); + assert_eq!(add(Growth::Unknown, n2.clone()), Growth::Unknown); + assert_eq!(add(n2.clone(), Growth::Unknown), Growth::Unknown); + assert_eq!(mul(Growth::Unknown, n2.clone()), Growth::Unknown); + assert_eq!(mul(n2, Growth::Unknown), Growth::Unknown); +} + +// --- Additional coverage --- + +/// Pure constants, constant factors, and constant division are all O(1) / dropped. +#[test] +fn test_growth_constants_are_o1() { + let c = g("42"); + assert_eq!(terms_of(&c), [GrowthTerm::one()]); + + // A wholly constant subtree (including `2^3`, `factorial(3)`, `1/2`) is O(1). + assert_eq!(g("2^3"), c); + assert_eq!(g("factorial(3)"), c); + + // Constant multiplier and constant divisor drop out. + assert_eq!(g("3 * n"), g("n")); + assert_eq!(g("n / 2"), g("n")); +} + +/// `x^0` is O(1); a negative exponent on a variable base is not admitted. +#[test] +fn test_growth_pow_special_cases() { + assert_eq!(terms_of(&g("n^0")), [GrowthTerm::one()]); + assert_eq!(g("n^(-1)"), Growth::Unknown); + // Variable base with variable exponent is not representable. + assert_eq!(g("n^m"), Growth::Unknown); +} + +/// Canonical Big-O rendering: bounded classes get `O()`, `Unknown` gets `O(?)`. +#[test] +fn test_growth_to_big_o() { + // The dominated `n` summand is dropped by the antichain, leaving just `n^2`. + assert_eq!(g("n^2 + n").to_big_o(), "O(n^2)"); + assert_eq!(g("2^n").to_big_o(), "O(2^n)"); + assert_eq!(g("5").to_big_o(), "O(1)"); + assert_eq!(Growth::Unknown.to_big_o(), "O(?)"); + // Renders exactly `O()` for bounded classes. + let bounded = g("n * m"); + assert_eq!( + bounded.to_big_o(), + format!("O({})", bounded.to_expr().unwrap()) + ); +} + +/// Exponential bases are authoritative symbolic data, not values reconstructed +/// from a rounded base-2 logarithm. +#[test] +fn test_growth_preserves_exponential_base() { + assert_eq!(g("3^n").to_big_o(), "O(3^n)"); + assert_eq!(g("1.0000000001^n").to_big_o(), "O(1.0000000001^n)"); + assert_eq!(g("2.7182818289^n").to_big_o(), "O(2.7182818289^n)"); + assert_eq!(g("2^(n / 2)").to_big_o(), "O(2^(0.5 * n))"); +} + +#[test] +fn test_growth_exponential_roundtrip_is_exact() { + for source in [ + "3^n", + "2^(n / 2)", + "exp(2 * n)", + "2^n * 3^n", + "3^n * n^2 * log(n)", + ] { + let growth = g(source); + let rendered = growth.to_expr().expect("growth should be representable"); + assert_eq!( + Growth::from_expr(&rendered), + growth, + "exponential growth changed while round-tripping {source} via {rendered}" + ); + } +} + +/// `exp(n)` uses base e; a decaying/unit base is bounded by O(1). +#[test] +fn test_growth_exponential_variants() { + // exp(n) is represented directly as e^n: exponential, dominates any polynomial. + let en = g("exp(n)"); + assert!(en.dominates(&g("n^5"))); + // 2^(n-m) ≤ 2^n after dropping the negative rate. + assert_eq!(g("2^(n - m)"), g("2^n")); + // Unit base is O(1); a decaying base with a growing exponent is O(1) too. + assert_eq!(g("1^n"), g("7")); + assert_eq!(g("0.5^n"), g("7")); + // A fractional base with a negative exponent grows and retains that exact + // symbolic base instead of being translated through a common logarithm. + assert_eq!(g("0.5^(-n)").to_big_o(), "O(0.5^(-1 * n))"); + assert!(g("0.5^(-n)").dominates(&g("n^100"))); +} + +/// `log` lowers each level: log of an exponential is linear, log of a +/// polynomial is a log, and log distributes over products as a sum. +#[test] +fn test_growth_log_levels() { + // log(2^n) ≍ n. + assert_eq!(g("log(2^n)"), g("n")); + assert_eq!(g("log(3^n)"), g("n")); + assert_eq!(g("log(exp(n))"), g("n")); + // log(n) is a single log term. + assert_eq!( + g("log(n)"), + Growth::Terms(vec![term(&[], &[], &[("n", 1)])]) + ); + // log(n*m) ≍ log n + log m (two summands, not a product). + assert_eq!(terms_of(&g("log(n*m)")).len(), 2); + // log of a constant is O(1). + assert_eq!(terms_of(&g("log(5)")), [GrowthTerm::one()]); + + // A mixed monomial's log keeps *every* factor class: log(2^n * m) ≍ n + log m. + // The exponential factor must not swallow the polynomial one. + let mixed = g("log(2^n * m)"); + let expected = make_growth(vec![ + term(&[], &[("n", 1.0)], &[]), + term(&[], &[], &[("m", 1)]), + ]); + assert_eq!(mixed, expected); + assert_eq!(terms_of(&mixed).len(), 2, "expected n + log m: {mixed:?}"); + + // When the classes share a variable the dominated summand is pruned: + // log(2^n * n^2) ≍ n + log n ≍ n (a single summand). + let shared = g("log(2^n * n^2)"); + assert_eq!(shared, g("n")); + assert_eq!(terms_of(&shared), [term(&[], &[("n", 1.0)], &[])]); +} + +/// `Unknown` is the top of the growth order. +#[test] +fn test_growth_unknown_dominance() { + let n2 = g("n^2"); + assert!(Growth::Unknown.dominates(&n2)); + assert!(!n2.dominates(&Growth::Unknown)); + assert!(Growth::Unknown.dominates(&Growth::Unknown)); +} + +/// On antichain-cap overflow the domain widens up to the single componentwise +/// max term (a valid upper bound), never truncating by iteration order. +#[test] +fn test_growth_antichain_cap_widens() { + // 40 distinct single-variable terms are pairwise incomparable. + let vars: Vec<&'static str> = (0..40) + .map(|i| &*Box::leak(format!("v{i}").into_boxed_str())) + .collect(); + let many: Vec = vars.iter().map(|v| term(&[], &[(*v, 1.0)], &[])).collect(); + + let widened = make_growth(many); + let ts = terms_of(&widened); + assert_eq!(ts.len(), 1, "cap overflow should widen to one term"); + // The single term dominates every original (it carries all variables). + for v in &vars { + assert!( + ts[0].dominates(&term(&[], &[(*v, 1.0)], &[])) || ts[0] == term(&[], &[(*v, 1.0)], &[]) + ); + } +} + +#[test] +fn test_growth_componentwise_max_with_symbolic_exponentials() { + let inputs = vec![ + terms_of(&g("2^n * n")).first().unwrap().clone(), + terms_of(&g("3^n * log(n)")).first().unwrap().clone(), + ]; + let upper = componentwise_max(&inputs).expect("3^n is a proven exponential maximum"); + assert!(inputs.iter().all(|term| upper.dominates_or_eq(term))); + assert_eq!(Growth::Terms(vec![upper]).to_big_o(), "O(3^n * n * log(n))"); + + let invalid = GrowthTerm { + exp: BTreeMap::new(), + poly: [("n", f64::NAN)].into_iter().collect(), + logs: BTreeMap::new(), + }; + assert_eq!(componentwise_max(&[invalid]), None); +} + +/// If symbolic exponential products have no provable componentwise maximum, +/// cap overflow widens to Unknown instead of guessing an under-bound. +#[test] +fn test_growth_antichain_cap_with_unproved_exponentials_is_unknown() { + let terms = (1..=33) + .map(|i| GrowthTerm { + exp: [("n", exp_product(&[(2.0, i as f64), (3.0, 1.0 / i as f64)]))] + .into_iter() + .collect(), + poly: BTreeMap::new(), + logs: BTreeMap::new(), + }) + .collect(); + + assert_eq!(make_growth(terms), Growth::Unknown); +} + +/// Structured serde round-trips (with `&'static str` keys leaked on read), and +/// `Unknown` round-trips. +#[test] +fn test_growth_serde_roundtrip() { + let value = g("2^n * m^2 + n * log(k)"); + let json = serde_json::to_string(&value).unwrap(); + let back: Growth = serde_json::from_str(&json).unwrap(); + assert_eq!(value, back); + + let unknown_json = serde_json::to_string(&Growth::Unknown).unwrap(); + assert_eq!( + serde_json::from_str::(&unknown_json).unwrap(), + Growth::Unknown + ); + + // Every constant Expr form admitted as a symbolic base remains lossless. + for source in [ + "(1 + 1)^n", + "(2 * 2)^n", + "(2^2)^n", + "exp(1)^n", + "log(3)^n", + "sqrt(4)^n", + "factorial(3)^n", + "exp(n)", + ] { + let value = g(source); + let json = serde_json::to_string(&value).unwrap(); + assert_eq!(serde_json::from_str::(&json).unwrap(), value); + } + + // The transient base-2-rate representation from the unmerged PR is not + // guessed back into a symbolic base. + let old_rate_only = r#"{"Terms":[{"exp":{"n":1.0},"poly":{},"logs":{}}]}"#; + assert!(serde_json::from_str::(old_rate_only).is_err()); + + let variable_base = r#"{"Constant":{"Var":"n"}}"#; + assert!(serde_json::from_str::(variable_base).is_err()); + + let invalid = Growth::Terms(vec![GrowthTerm { + exp: [("n", ExpProduct::empty())].into_iter().collect(), + poly: BTreeMap::new(), + logs: BTreeMap::new(), + }]); + let invalid_json = serde_json::to_string(&invalid).unwrap(); + assert!(serde_json::from_str::(&invalid_json).is_err()); +} + +// --- Randomized property tests (#1077) --- +// +// These cross-validate the symbolic growth domain against the numeric ground +// truth (`Expr::eval`) over a large, seeded input space, in the spirit of the +// repo's `/verify-reduction` adversarial culture. Three contracts are exercised +// ≥ 5000 times each with a hand-rolled, deterministic RNG (no wall-clock, no +// entropy — CI must be byte-reproducible across platforms): +// +// 1. Upper-bound soundness: `eval(e, s) ≤ C·eval(render(growth(e)), s)` at +// sizes larger than the anchor from which `C` was calibrated. +// 2. Idempotence: `growth(render(growth(e))) == growth(e)`. +// 3. Dominance soundness: when `dominates(b, a)`, the numeric ratio +// `eval(b)/eval(a)` does not shrink and exceeds 1 at the larger size. +// +// A #[test] negative control runs the same upper-bound harness against a +// deliberately broken transfer function and asserts the harness catches it, so +// the property tests are demonstrably capable of failing. +// +// Why the domain exists at all is *why* some numeric checks are unreachable: +// crossovers like `2^n ≻ n^100` lie far beyond f64 range. The harnesses handle +// this honestly — they skip (and count) samples where numerics are +// indeterminate (both sides overflow to `inf`), never by hiding a failing +// assertion. The dominance contract additionally restricts its numeric +// cross-check to single-term, in-band growths, the regime where the crossover +// is reachable; that regime targets exactly the lexicographic per-variable +// comparison (`GrowthTerm::cmp`) at the heart of the order, so the restriction +// is well-aimed, not vacuous. + +use super::{exponential, log_growth, pow_const}; +use crate::types::ProblemSize; +use std::collections::BTreeMap; + +/// Fixed master seed. Every contract derives its own stream by offsetting this, +/// so the whole suite is deterministic and reproducible on any platform. +const MASTER_SEED: u64 = 0xD1CE_2026_1077_ABCD; + +/// SplitMix64 — a tiny, fully specified PRNG. Hand-rolled (rather than +/// `rand::StdRng`) precisely because its output must be identical across crate +/// versions and platforms; the constants below are the published SplitMix64 +/// mixing constants and will never change. +struct SplitMix64 { + state: u64, +} + +impl SplitMix64 { + fn new(seed: u64) -> Self { + SplitMix64 { state: seed } + } + + fn next_u64(&mut self) -> u64 { + self.state = self.state.wrapping_add(0x9E37_79B9_7F4A_7C15); + let mut z = self.state; + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + z ^ (z >> 31) + } + + /// Uniform integer in `[0, n)`. + fn below(&mut self, n: u64) -> u64 { + self.next_u64() % n + } +} + +fn b(e: Expr) -> Box { + Box::new(e) +} + +/// Variable pool — `&'static str` literals so they satisfy `Expr::Var` and match +/// the `ProblemSize` keys built by [`joint_size`]. +const VARS: [&str; 3] = ["n", "m", "k"]; + +fn gen_var(rng: &mut SplitMix64) -> Expr { + Expr::Var(VARS[rng.below(VARS.len() as u64) as usize]) +} + +/// All variables set jointly to `s` (the contracts evaluate on the diagonal). +fn joint_size(s: usize) -> ProblemSize { + ProblemSize::new(vec![("n", s), ("m", s), ("k", s)]) +} + +// --- General expression generator (contracts 1 and 2) --- +// +// Bounded depth, variables {n, m, k}, constructors Const/Var/Add/Mul/Pow(const)/ +// Sqrt/Log plus linear `2^x` and `exp(x)` forms. A small (~1% per node) branch +// emits a nonlinear exponent (`2^(n*m)`, `2^sqrt(n)`) so the `Unknown` widening +// path is genuinely exercised while staying a minority of whole trees. + +const MAX_DEPTH: u32 = 5; + +fn gen_leaf(rng: &mut SplitMix64) -> Expr { + // Bias toward variables; keep constants small and positive. + if rng.below(4) == 0 { + Expr::Const((1 + rng.below(4)) as f64) + } else { + gen_var(rng) + } +} + +/// A linear expression in the variables (so `2^x` stays first-class in the +/// domain): a sum of 1..=3 terms `c·v` with small positive integer coefficients. +fn gen_linear(rng: &mut SplitMix64) -> Expr { + let nterms = 1 + rng.below(3); + let mut e = gen_lin_term(rng); + for _ in 1..nterms { + e = e + gen_lin_term(rng); + } + e +} + +fn gen_lin_term(rng: &mut SplitMix64) -> Expr { + let v = gen_var(rng); + let c = 1 + rng.below(3); + if c == 1 { + v + } else { + Expr::Const(c as f64) * v + } +} + +/// A deliberately nonlinear exponent, driving `2^(·)` to `Growth::Unknown`. +fn gen_nonlinear(rng: &mut SplitMix64) -> Expr { + if rng.below(2) == 0 { + Expr::Mul(b(gen_var(rng)), b(gen_var(rng))) + } else { + Expr::Sqrt(b(gen_var(rng))) + } +} + +const E_BELOW: f64 = std::f64::consts::E - 1e-10; +const E_ABOVE: f64 = std::f64::consts::E + 1e-10; +const STABLE_EXPONENTIAL_BASES: &[f64] = &[2.0, E_BELOW, E_ABOVE, 3.0]; +const ADVERSARIAL_EXPONENTIAL_BASES: &[f64] = &[1.0000000001, 2.0, E_BELOW, E_ABOVE, 3.0]; + +fn gen_exponential_base(rng: &mut SplitMix64, bases: &[f64]) -> Expr { + Expr::Const(bases[rng.below(bases.len() as u64) as usize]) +} + +fn gen_expr(rng: &mut SplitMix64, depth: u32, exponential_bases: &[f64]) -> Expr { + if depth == 0 { + return gen_leaf(rng); + } + match rng.below(100) { + 0..=19 => gen_leaf(rng), + 20..=39 => Expr::Add( + b(gen_expr(rng, depth - 1, exponential_bases)), + b(gen_expr(rng, depth - 1, exponential_bases)), + ), + 40..=54 => Expr::Mul( + b(gen_expr(rng, depth - 1, exponential_bases)), + b(gen_expr(rng, depth - 1, exponential_bases)), + ), + 55..=69 => Expr::pow( + gen_expr(rng, depth - 1, exponential_bases), + Expr::Const((1 + rng.below(3)) as f64), + ), + 70..=79 => Expr::Sqrt(b(gen_expr(rng, depth - 1, exponential_bases))), + 80..=89 => Expr::Log(b(gen_expr(rng, depth - 1, exponential_bases))), + 90..=96 => Expr::pow( + gen_exponential_base(rng, exponential_bases), + gen_linear(rng), + ), + 97..=98 => Expr::Exp(b(gen_var(rng))), + // ~1% per node: a nonlinear exponent → Unknown (a minority of trees). + _ => Expr::pow( + gen_exponential_base(rng, exponential_bases), + gen_nonlinear(rng), + ), + } +} + +// --- Monomial generator (contract 3) --- +// +// A product of single-term factors, so its growth is always a single antichain +// term. This isolates the lexicographic per-variable dominance decision. + +fn gen_factor(rng: &mut SplitMix64) -> Expr { + let v = gen_var(rng); + match rng.below(6) { + 0 => v, + 1 => Expr::pow(v, Expr::Const((1 + rng.below(3)) as f64)), + 2 => Expr::Sqrt(b(v)), + 3 => Expr::Log(b(v)), + // Keep the numeric dominance harness on one common base: different + // fixed bases can have crossovers beyond its finite observation window. + // Multi-base behavior is covered by symbolic proof tests above. + 4 => Expr::pow(Expr::Const(2.0), v), + _ => Expr::pow(Expr::Const(2.0), Expr::Const((1 + rng.below(3)) as f64) * v), + } +} + +fn gen_monomial(rng: &mut SplitMix64) -> Expr { + let nf = 1 + rng.below(4); + let mut e = gen_factor(rng); + for _ in 1..nf { + e = e * gen_factor(rng); + } + e +} + +// --- Contract 1: upper-bound soundness --- + +/// The number of independent `#[test]`-level iterations for the upper-bound and +/// idempotence contracts (each well above the 5000-meaningful-check floor after +/// `Unknown`/overflow skips). +const UB_ITERS: usize = 20_000; + +/// Outcome tallies for the upper-bound harness. `meaningful` counts samples that +/// produced at least one *conclusive* large-size comparison. +#[derive(Default)] +struct UbResult { + meaningful: usize, + unknown: usize, + skipped: usize, + violations: usize, + first_violation: Option, +} + +/// Run the upper-bound harness against an arbitrary transfer function. The real +/// test passes `Growth::from_expr`; the negative control passes +/// `broken_from_expr`. Parameterizing here is what gives the harness teeth: the +/// exact same code must accept the sound transfer and reject the broken one. +fn run_upper_bound(transfer: fn(&Expr) -> Growth, seed: u64, iters: usize) -> UbResult { + // Anchor 2^6; check at 2^8, 2^10, 2^12 — all *larger* than the anchor. + let anchor = 64.0_f64; + let large = [256.0_f64, 1024.0, 4096.0]; + let slack = 16.0_f64; + + let mut rng = SplitMix64::new(seed); + let mut r = UbResult::default(); + + for _ in 0..iters { + let e = gen_expr(&mut rng, MAX_DEPTH, STABLE_EXPONENTIAL_BASES); + let g = transfer(&e); + let gexpr = match g.to_expr() { + Some(x) => x, + None => { + r.unknown += 1; + continue; + } + }; + + // Calibrate C from the observed ratio at the (smaller) anchor. + let sz0 = joint_size(anchor as usize); + let ve0 = e.eval(&sz0); + let vg0 = gexpr.eval(&sz0); + // Nonnegativity is a domain precondition. A negative anchor value means + // the generated expression is outside the domain's contract (e.g. deeply + // nested `log`s that are negative at these sizes) — skip it, don't hold + // the domain to a bound it never promised for such inputs. + if !ve0.is_finite() || !vg0.is_finite() || ve0 <= 0.0 || vg0 <= 0.0 { + r.skipped += 1; + continue; + } + let c = (ve0 / vg0) * slack; + + let mut conclusive = false; + for &s in &large { + let sz = joint_size(s as usize); + let ve = e.eval(&sz); + let vg = gexpr.eval(&sz); + if ve.is_nan() || vg.is_nan() { + continue; + } + if vg.is_infinite() { + // The bound overestimates. Holds trivially unless `e` also blew + // up, in which case the comparison is indeterminate — skip it. + if ve.is_finite() { + conclusive = true; + } + continue; + } + if ve.is_infinite() { + // `eval(e)` can overflow to `inf` at intermediate steps even + // when the true value is finite (e.g. `log(n^2 * exp(n))` blows + // up at the inner `exp` before the outer `log` tames it back to + // `n`). Such a numeric artifact is indeterminate, not a genuine + // violation of a finite bound — skip this size. + continue; + } + if ve <= 0.0 || vg <= 0.0 { + // Out of the nonnegative domain at this size — indeterminate. + continue; + } + // Both finite and positive: a real, decidable comparison. + conclusive = true; + let bound = c * vg; + if ve > bound { + r.violations += 1; + if r.first_violation.is_none() { + r.first_violation = Some(format!( + "e = {e} | g = {gexpr} | s = {s}: eval(e) = {ve} > {c} * {vg} = {bound}" + )); + } + } + } + + if conclusive { + r.meaningful += 1; + } else { + r.skipped += 1; + } + } + r +} + +/// A deliberately broken transfer function: `Add` keeps only its *first* +/// operand's growth, dropping the second. This is an under-approximation — it +/// can miss the dominant summand — so the upper bound must fail somewhere. +/// Every other node mirrors the real `Growth::from_expr` (reusing its private +/// transfer helpers), so the only defect is the seeded `Add` bug. +fn broken_from_expr(e: &Expr) -> Growth { + if e.constant_value().is_some() { + return Growth::Terms(vec![GrowthTerm::one()]); + } + match e { + Expr::Const(_) => Growth::Terms(vec![GrowthTerm::one()]), + Expr::Var(v) => { + let mut t = GrowthTerm::one(); + t.poly.insert(v, 1.0); + Growth::Terms(vec![t]) + } + // The seeded bug: drop the second summand. + Expr::Add(a, _b) => broken_from_expr(a), + Expr::Mul(a, b) => mul(broken_from_expr(a), broken_from_expr(b)), + Expr::Pow(base, exp) => { + if let Some(k) = exp.constant_value() { + if k < 0.0 { + Growth::Unknown + } else if k == 0.0 { + Growth::Terms(vec![GrowthTerm::one()]) + } else { + pow_const(broken_from_expr(base), k) + } + } else if base.constant_value().is_some() { + exponential(ExpBase::Constant(base.as_ref().clone()), exp) + } else { + Growth::Unknown + } + } + Expr::Exp(a) => exponential(ExpBase::Natural, a), + Expr::Log(a) => log_growth(broken_from_expr(a)), + Expr::Sqrt(a) => pow_const(broken_from_expr(a), 0.5), + Expr::Factorial(_) => Growth::Unknown, + } +} + +#[test] +fn test_growth_property_upper_bound_sound() { + let r = run_upper_bound(Growth::from_expr, MASTER_SEED ^ 0x01, UB_ITERS); + + assert_eq!( + r.violations, + 0, + "upper-bound violation ({} total); first: {}", + r.violations, + r.first_violation.as_deref().unwrap_or("") + ); + assert!( + r.meaningful >= 5000, + "need >= 5000 meaningful checks, got {} (unknown {}, skipped {})", + r.meaningful, + r.unknown, + r.skipped + ); + // The generator must actually exercise the domain, not mostly produce Unknown. + let total = r.meaningful + r.unknown + r.skipped; + assert!( + r.unknown * 2 < total, + "Unknown must be a minority: {}/{}", + r.unknown, + total + ); + assert!(r.unknown > 0, "generator never exercised the Unknown path"); +} + +#[test] +fn test_growth_property_upper_bound_negative_control() { + // The SAME harness, run against the broken transfer, must detect a + // violation. If it cannot, the property tests have no teeth and this fails. + let r = run_upper_bound(broken_from_expr, MASTER_SEED ^ 0x01, UB_ITERS); + assert!( + r.violations > 0, + "harness failed to catch the seeded Add bug (meaningful {}, violations {})", + r.meaningful, + r.violations + ); +} + +// --- Contract 2: idempotence --- + +/// Exponential factors round-trip exactly. Polynomial degrees retain the +/// pre-existing tolerance for unrelated floating-point power composition. +fn map_approx_eq(a: &BTreeMap<&'static str, f64>, b: &BTreeMap<&'static str, f64>) -> bool { + a.len() == b.len() + && a.iter() + .all(|(k, v)| b.get(k).is_some_and(|w| (v - w).abs() < 1e-6)) +} + +fn term_approx_eq(x: &GrowthTerm, y: &GrowthTerm) -> bool { + x.exp == y.exp && map_approx_eq(&x.poly, &y.poly) && x.logs == y.logs +} + +fn growth_approx_eq(a: &Growth, b: &Growth) -> bool { + match (a, b) { + (Growth::Unknown, Growth::Unknown) => true, + (Growth::Terms(ta), Growth::Terms(tb)) => { + ta.len() == tb.len() + && ta.iter().all(|t| tb.iter().any(|u| term_approx_eq(t, u))) + && tb.iter().all(|u| ta.iter().any(|t| term_approx_eq(t, u))) + } + _ => false, + } +} + +#[test] +fn test_growth_property_idempotence() { + let mut rng = SplitMix64::new(MASTER_SEED ^ 0x02); + let mut meaningful = 0usize; + let mut unknown = 0usize; + + for _ in 0..UB_ITERS { + // Idempotence is purely symbolic, so it can safely exercise bases near + // one whose numeric crossover lies far beyond the f64 test window. + let e = gen_expr(&mut rng, MAX_DEPTH, ADVERSARIAL_EXPONENTIAL_BASES); + let g = Growth::from_expr(&e); + let rendered = match g.to_expr() { + Some(x) => x, + None => { + unknown += 1; + continue; + } + }; + let g2 = Growth::from_expr(&rendered); + assert!( + growth_approx_eq(&g, &g2), + "growth not idempotent: e = {e} | render = {rendered}\n g = {g:?}\n g2 = {g2:?}" + ); + meaningful += 1; + } + + assert!( + meaningful >= 5000, + "need >= 5000 meaningful checks, got {meaningful} (unknown {unknown})" + ); +} + +// --- Contract 3: dominance soundness --- + +const DOM_ITERS: usize = 120_000; + +/// A single antichain term, or `None` if the growth is `Unknown` or a +/// multi-term antichain. Restricting to single terms keeps the numeric ratio a +/// pure monomial ratio: multi-term dominance can add a *lower-order* summand +/// (`{n^2, m}` dominates `{n^2}`) whose ratio shrinks toward 1 — a real feature +/// of the antichain order, but not what this monomial cross-check targets. The +/// single-term regime isolates the lexicographic per-variable comparison +/// (`GrowthTerm::cmp`) that is the heart of the order. +fn single_term(g: &Growth) -> Option<&GrowthTerm> { + match g { + Growth::Terms(ts) if ts.len() == 1 => Some(&ts[0]), + _ => None, + } +} + +/// `(total exp rate, total poly degree, total log power)` on the joint diagonal. +fn totals(t: &GrowthTerm) -> (f64, f64, f64) { + ( + t.exp.values().map(ExpProduct::log2_estimate).sum(), + t.poly.values().sum(), + t.logs.values().map(|&x| x as f64).sum(), + ) +} + +#[test] +fn test_growth_property_dominance_sound() { + let mut rng = SplitMix64::new(MASTER_SEED ^ 0x03); + let mut meaningful = 0usize; + let mut skipped = 0usize; + let mut unreachable = 0usize; + const LN2: f64 = std::f64::consts::LN_2; + + for _ in 0..DOM_ITERS { + let ga = Growth::from_expr(&gen_monomial(&mut rng)); + let gb = Growth::from_expr(&gen_monomial(&mut rng)); + + let (ta, tb) = match (single_term(&ga), single_term(&gb)) { + (Some(a), Some(b)) => (a.clone(), b.clone()), + _ => { + skipped += 1; + continue; + } + }; + + // Orient to the strict dominator; skip incomparable or asymptotically + // equal pairs (a flat ratio has nothing to assert). + let ab = ga.dominates(&gb); + let ba = gb.dominates(&ga); + let (hi, lo) = if ba && !ab { + (&tb, &ta) + } else if ab && !ba { + (&ta, &tb) + } else { + skipped += 1; + continue; + }; + + // Choose the evaluation window from the *magnitude* of the exponent gap + // — a structural property of the two terms, computed independently of + // which direction `dominates` picked. This places the check in the + // numerically-informative regime (past the ratio's minimum, past the + // crossover, below f64 overflow) so the assertions are meaningful; it + // does NOT peek at the assertion outcome, so a mis-ordering by + // `dominates` still fails the signed check below. + let (eh, ph, lh) = totals(hi); + let (el, pl, ll) = totals(lo); + let (de, dp, dl) = (eh - el, ph - pl, lh - ll); + const EPS: f64 = 1e-9; + let exp_max = eh.max(el); + + let (s1, s2): (usize, usize) = if de.abs() > EPS { + // Exponential gap: crossover is at moderate size; keep exp finite. + (16, 64) + } else if dp.abs() > EPS { + // Polynomial gap under a *common* exponent: the crossover (e.g. + // sqrt(n) vs (log n)^3 at n≈2.4e7) needs large sizes where any + // shared exponential would overflow. Reachable only with no + // exponential — and then poly values stay finite to astronomical + // sizes, so a wide window clears even the fractional-poly-vs-high- + // log-power crossovers our generator can produce (dp≥0.5, |dl|≤4). + if exp_max > EPS { + unreachable += 1; + continue; + } + (8192, 1usize << 42) + } else if dl.abs() > EPS { + // Log-power gap only: manifest at any modest size. + (16, 64) + } else { + // No gap on the diagonal (strict domination on an off-diagonal + // variable that collapses here) — nothing to assert numerically. + skipped += 1; + continue; + }; + + // Overflow guard for the (in-principle reachable) exponential cases. + if exp_max * (s2 as f64) * LN2 > 700.0 { + unreachable += 1; + continue; + } + + let a = Growth::Terms(vec![lo.clone()]).to_expr().unwrap(); + let bx = Growth::Terms(vec![hi.clone()]).to_expr().unwrap(); + let (z1, z2) = (joint_size(s1), joint_size(s2)); + let (a1, a2) = (a.eval(&z1), a.eval(&z2)); + let (b1, b2) = (bx.eval(&z1), bx.eval(&z2)); + if [a1, a2, b1, b2].iter().any(|v| !v.is_finite() || *v <= 0.0) { + skipped += 1; + continue; + } + + let r1 = b1 / a1; + let r2 = b2 / a2; + meaningful += 1; + + // The ratio does not shrink from s1 to s2 (tiny tolerance for float + // noise), and it exceeds 1 at the larger size. A wrong-direction + // dominance decision flips the signed gap and fails both. + assert!( + r2 >= r1 * (1.0 - 1e-9), + "dominance ratio shrank: {bx} over {a}; r({s1}) = {r1}, r({s2}) = {r2}" + ); + assert!( + r2 > 1.0, + "dominator not numerically ahead at s2: {bx} over {a}; r({s2}) = {r2}" + ); + } + + assert!( + meaningful >= 5000, + "need >= 5000 meaningful dominating pairs, got {meaningful} \ + (skipped {skipped}, unreachable {unreachable})" + ); +} diff --git a/src/unit_tests/reduction_graph.rs b/src/unit_tests/reduction_graph.rs index be612b57b..88454a153 100644 --- a/src/unit_tests/reduction_graph.rs +++ b/src/unit_tests/reduction_graph.rs @@ -988,3 +988,59 @@ fn test_find_paths_bounded_limits_depth() { "MIS→QUBO has no direct edge, so bound=0 should return empty" ); } + +#[test] +fn test_find_paths_bounded_returns_shortest_when_truncated() { + use crate::expr::Expr; + use crate::rules::registry::{EdgeCapabilities, ReductionOverhead}; + use crate::rules::ReductionEdgeData; + + fn edge() -> ReductionEdgeData { + ReductionEdgeData { + overhead: ReductionOverhead::new(vec![("n", Expr::Var("n"))]), + reduce_fn: None, + reduce_aggregate_fn: None, + capabilities: EdgeCapabilities::witness_only(), + } + } + + // Topology where DFS discovery order surfaces a LONG route before the SHORT one. + // From S the first outgoing edge (S->A) leads into a long chain A->B->C->T, while a + // later edge S->T is a direct hop. petgraph's DFS explores S->A first, so the + // 4-edge route is discovered before the 1-edge direct route. With a tight limit, + // the old `.take(limit)` in discovery order would keep the long route and drop the + // short one; length-first enumeration must return the short route. + let graph = ReductionGraph::from_test_edges( + &["S", "A", "B", "C", "T"], + &[ + ("S", "A", edge()), + ("A", "B", edge()), + ("B", "C", edge()), + ("C", "T", edge()), + ("S", "T", edge()), + ], + ); + + let empty = BTreeMap::new(); + + // Sanity: both routes exist when unbounded. + let all = graph.find_paths_up_to("S", &empty, "T", &empty, 100); + assert_eq!(all.len(), 2, "expected the direct route and the long chain"); + + // With limit 1, the SHORT (direct) route must be the one returned. + let limited = graph.find_paths_up_to("S", &empty, "T", &empty, 1); + assert_eq!(limited.len(), 1); + assert_eq!( + limited[0].len(), + 1, + "truncated result must keep the shortest (direct) route, not the long chain" + ); + + // Results are length-sorted (non-decreasing edge counts). + let lens: Vec = all.iter().map(|p| p.len()).collect(); + assert!( + lens.windows(2).all(|w| w[0] <= w[1]), + "paths must be returned shortest-first, got lengths {lens:?}" + ); + assert_eq!(lens, vec![1, 4]); +} diff --git a/src/unit_tests/rules/analysis.rs b/src/unit_tests/rules/analysis.rs index a97f6060e..54cff218b 100644 --- a/src/unit_tests/rules/analysis.rs +++ b/src/unit_tests/rules/analysis.rs @@ -71,45 +71,87 @@ fn test_compare_overhead_no_common_fields() { } #[test] -fn test_compare_overhead_unknown_exp() { - // Different exponential-vs-polynomial growth is still not decided by the - // monomial comparison fallback. +fn test_compare_overhead_exp_dominates_poly() { + // primitive exp(n) grows faster than composite n, so composite ≤ primitive + // on the only common field → dominated. (The old polynomial engine rejected + // exp outright and returned Unknown; the growth domain decides it.) let prim = ReductionOverhead::new(vec![("num_vars", Expr::Exp(Box::new(Expr::Var("n"))))]); let comp = ReductionOverhead::new(vec![("num_vars", Expr::Var("n"))]); - assert_eq!(compare_overhead(&prim, &comp), ComparisonStatus::Unknown); + assert_eq!(compare_overhead(&prim, &comp), ComparisonStatus::Dominated); } #[test] -fn test_compare_overhead_unknown_log() { +fn test_compare_overhead_poly_dominates_log() { + // primitive n vs composite log(n): n grows faster than log(n), so the + // composite is dominated. Previously Unknown (the polynomial engine could + // not normalize `log`); now decided by the growth domain. let prim = ReductionOverhead::new(vec![("num_vars", Expr::Var("n"))]); let comp = ReductionOverhead::new(vec![("num_vars", Expr::Log(Box::new(Expr::Var("n"))))]); - assert_eq!(compare_overhead(&prim, &comp), ComparisonStatus::Unknown); + assert_eq!(compare_overhead(&prim, &comp), ComparisonStatus::Dominated); } #[test] -fn test_compare_overhead_exp_identity_after_asymptotic_normalization() { +fn test_compare_overhead_exp_identity_decided() { + // `exp(n + m)` and `exp(n) * exp(m)` are asymptotically equal. The growth + // domain normalizes both to the same exponential term, so the (reflexive) + // dominance holds → dominated. (Was temporarily asserted Unknown while the + // bespoke engine — which could not handle exp — was still in place.) let prim = ReductionOverhead::new(vec![("num_vars", Expr::parse("exp(n + m)"))]); let comp = ReductionOverhead::new(vec![("num_vars", Expr::parse("exp(n) * exp(m)"))]); assert_eq!(compare_overhead(&prim, &comp), ComparisonStatus::Dominated); } #[test] -fn test_compare_overhead_log_identity_after_asymptotic_normalization() { - // log(n) vs log(n^2): the new canonicalization engine keeps log(n^2) as-is - // (it doesn't simplify log(x^k) = k*log(x)), so polynomial comparison - // returns Unknown for non-polynomial log terms. +fn test_compare_overhead_log_identity_decided() { + // log(n) vs log(n^2): the growth domain uses log(n^k) ≍ log(n), so both + // fields collapse to the same growth → dominated. (Was temporarily Unknown + // because the polynomial engine could not normalize `log`.) let prim = ReductionOverhead::new(vec![("num_vars", Expr::parse("log(n)"))]); let comp = ReductionOverhead::new(vec![("num_vars", Expr::parse("log(n^2)"))]); - assert_eq!(compare_overhead(&prim, &comp), ComparisonStatus::Unknown); + assert_eq!(compare_overhead(&prim, &comp), ComparisonStatus::Dominated); } #[test] -fn test_compare_overhead_sqrt_identity_after_asymptotic_normalization() { +fn test_compare_overhead_sqrt_identity_decided() { + // `sqrt(n * m)` and `(n * m)^(1/2)` are equal; the growth domain maps both + // to poly degree 0.5 in n and m → dominated. (Was temporarily Unknown while + // the sqrt-rejecting polynomial engine was in place.) let prim = ReductionOverhead::new(vec![("num_vars", Expr::parse("sqrt(n * m)"))]); let comp = ReductionOverhead::new(vec![("num_vars", Expr::parse("(n * m)^(1/2)"))]); assert_eq!(compare_overhead(&prim, &comp), ComparisonStatus::Dominated); } +#[test] +fn test_compare_overhead_subtraction_now_decided() { + // Subtraction: primitive n^2 vs composite n^2 - n. The growth domain widens + // `a - b ⇝ a + b` so n^2 - n ≍ n^2, asymptotically equal to the primitive → + // dominated. The old polynomial engine rejected negative coefficients and + // returned Unknown. + let prim = ReductionOverhead::new(vec![("num_vars", Expr::parse("n^2"))]); + let comp = ReductionOverhead::new(vec![("num_vars", Expr::parse("n^2 - n"))]); + assert_eq!(compare_overhead(&prim, &comp), ComparisonStatus::Dominated); +} + +#[test] +fn test_compare_overhead_negative_control_cubic_worse() { + // Negative control: primitive num_vertices = n^2 vs composite num_vertices = + // n^3, all other common fields equal. The composite grows strictly faster on + // the differing field, so this MUST be NotDominated — a direction inversion + // or an ignored field would flip it to Dominated. + let prim = ReductionOverhead::new(vec![ + ("num_vertices", Expr::pow(Expr::Var("n"), Expr::Const(2.0))), + ("num_edges", Expr::Var("n")), + ]); + let comp = ReductionOverhead::new(vec![ + ("num_vertices", Expr::pow(Expr::Var("n"), Expr::Const(3.0))), + ("num_edges", Expr::Var("n")), + ]); + assert_eq!( + compare_overhead(&prim, &comp), + ComparisonStatus::NotDominated + ); +} + #[test] fn test_compare_overhead_additive_constant_after_asymptotic_normalization() { let prim = ReductionOverhead::new(vec![("num_vars", Expr::parse("n"))]); @@ -119,10 +161,9 @@ fn test_compare_overhead_additive_constant_after_asymptotic_normalization() { #[test] fn test_compare_overhead_multivariate_product_vs_sum() { - // n * m (degree 2) vs n + m (degree 1): - // monomial n*m has exponents {n:1, m:1} - // monomials n, m each have exponent 1 in one variable - // n*m is NOT dominated by either n or m → composite is worse + // primitive n + m ≍ {n, m} (two incomparable terms) vs composite n * m ≍ + // {n·m}. The single composite term n·m is dominated by neither n nor m, so + // the primitive does not dominate the composite → not dominated. let prim = ReductionOverhead::new(vec![("num_vars", Expr::Var("n") + Expr::Var("m"))]); let comp = ReductionOverhead::new(vec![("num_vars", Expr::Var("n") * Expr::Var("m"))]); assert_eq!( @@ -132,10 +173,10 @@ fn test_compare_overhead_multivariate_product_vs_sum() { } #[test] -fn test_compare_overhead_multivariate_product_vs_square() { - // n * m (has m) vs n^2 (no m): incomparable - // n*m monomial {n:1, m:1} — dominated by n^2 {n:2}? - // exponent_n: 1 <= 2 ✓, exponent_m: 1 <= 0 ✗ → not dominated +fn test_compare_overhead_incomparable_field_not_dominated() { + // Incomparable growths on a field: primitive n^2 vs composite n * m. n^2 has + // degree 2 in n and 0 in m; n·m has degree 1 in each. Neither dominates the + // other (n^2 wins on n, n·m wins on m) → not dominated. let prim = ReductionOverhead::new(vec![( "num_vars", Expr::pow(Expr::Var("n"), Expr::Const(2.0)), @@ -165,11 +206,11 @@ fn test_compare_overhead_constant_factor() { #[test] fn test_compare_overhead_polynomial_expansion() { - // (n + m)^2 = n^2 + 2nm + m^2 (degree 2) vs n^3 (degree 3) - // Each monomial of composite has total degree ≤ 2, primitive has degree 3 - // n^2 dominated by n^3? exponent_n: 2 ≤ 3 ✓ → yes - // 2*n*m dominated by n^3? exponent_n: 1 ≤ 3 ✓, exponent_m: 1 ≤ 0 ✗ → no! - // So composite is NOT dominated — (n+m)^2 can exceed n^3 when m is large + // Composite (n + m)^2 ≍ max(n, m)^2 = {n^2, m^2} in the growth domain (no + // binomial cross term). Primitive n^3 ≍ {n^3}. n^3 dominates n^2, but n^3 + // does not dominate m^2 (it has degree 0 in m), so the primitive does not + // dominate the composite → not dominated — (n+m)^2 can exceed n^3 when m is + // large. let prim = ReductionOverhead::new(vec![( "num_vars", Expr::pow(Expr::Var("n"), Expr::Const(3.0)), @@ -271,6 +312,16 @@ fn test_find_dominated_rules_returns_known_set() { "KSatisfiability {k: \"K3\"}", "MinimumVertexCover {graph: \"SimpleGraph\", weight: \"i32\"}", ), + // Newly decided by the growth-domain rewire (#1081): PartitionIntoPathsOfLength2 + // → BCSF → ILP{i32} → ILP{bool}. The composite's composed num_vars/num_constraints + // carry a `num_vertices / 3` factor (from max_components = V/3); the old polynomial + // engine rejected that constant divisor as a negative-exponent power and returned + // Unknown, while the growth domain drops constant divisors, giving both fields + // growth {V^2, E*V} — asymptotically equal to the direct edge, hence Dominated. + ( + "PartitionIntoPathsOfLength2 {graph: \"SimpleGraph\"}", + "ILP {variable: \"bool\"}", + ), ] .into_iter() .collect(); diff --git a/src/unit_tests/rules/pareto.rs b/src/unit_tests/rules/pareto.rs new file mode 100644 index 000000000..1fc843570 --- /dev/null +++ b/src/unit_tests/rules/pareto.rs @@ -0,0 +1,1292 @@ +//! Tests for the Pareto label-setting search (`src/rules/pareto.rs`) and its two label +//! domains. Covers: +//! - The measured concrete-instance search (issue #788 known-answer and budget semantics). +//! - The generic kernel's correctness on a hand-built diamond (negative control): a +//! scalar-cost path selection commits to the wrong prefix, while the Pareto search +//! returns the path with the strictly-better final measured size. + +use super::*; +use crate::expr::Expr; +use crate::growth::Growth; +use crate::models::algebraic::{ObjectiveSense, ILP}; +use crate::models::formula::{CNFClause, Satisfiability}; +use crate::models::graph::HamiltonianCircuit; +use crate::rules::cost::CustomCost; +use crate::rules::pareto::{GrowthLabel, PathLabel, ReductionEdge}; +use crate::rules::registry::{EdgeCapabilities, ReductionOverhead}; +use crate::rules::traits::DynReductionResult; +use crate::rules::{ReductionAutoCast, ReductionGraph, ReductionMode}; +use crate::topology::SimpleGraph; +use crate::traits::Problem; +use crate::types::{Or, ProblemSize}; +use std::any::Any; +use std::cell::Cell; +use std::collections::BTreeMap; +use std::rc::Rc; + +#[derive(Clone)] +struct MeasuredSource; + +#[derive(Clone)] +struct MeasuredBranchA; + +#[derive(Clone)] +struct MeasuredBranchB; + +macro_rules! impl_measured_test_problem { + ($ty:ty, $name:literal) => { + impl Problem for $ty { + const NAME: &'static str = $name; + type Value = Or; + + fn dims(&self) -> Vec { + vec![] + } + + fn evaluate(&self, _config: &[usize]) -> Or { + Or(true) + } + + fn variant() -> Vec<(&'static str, &'static str)> { + vec![] + } + } + }; +} + +impl_measured_test_problem!(MeasuredSource, "MeasuredSource"); +impl_measured_test_problem!(MeasuredBranchA, "MeasuredBranchA"); +impl_measured_test_problem!(MeasuredBranchB, "MeasuredBranchB"); + +fn measured_source_to_a(any: &dyn Any) -> Box { + any.downcast_ref::() + .expect("expected MeasuredSource"); + Box::new(ReductionAutoCast::::new( + MeasuredBranchA, + )) +} + +fn measured_source_to_b(any: &dyn Any) -> Box { + any.downcast_ref::() + .expect("expected MeasuredSource"); + Box::new(ReductionAutoCast::::new( + MeasuredBranchB, + )) +} + +fn measured_a_to_sat(any: &dyn Any) -> Box { + any.downcast_ref::() + .expect("expected MeasuredBranchA"); + Box::new(ReductionAutoCast::::new( + Satisfiability::new(1, vec![CNFClause::new(vec![1])]), + )) +} + +fn measured_b_to_sat(any: &dyn Any) -> Box { + any.downcast_ref::() + .expect("expected MeasuredBranchB"); + Box::new(ReductionAutoCast::::new( + Satisfiability::new(1, vec![CNFClause::new(vec![-1])]), + )) +} + +fn measured_sat_to_structure_dependent_ilp(any: &dyn Any) -> Box { + let sat = any + .downcast_ref::() + .expect("expected Satisfiability"); + let first_literal = sat.clauses()[0].literals[0]; + let num_vars = if first_literal > 0 { 100 } else { 1 }; + let target = ILP::::new(num_vars, vec![], vec![], ObjectiveSense::Minimize); + Box::new(ReductionAutoCast::>::new(target)) +} + +fn measured_source_to_small_ilp(any: &dyn Any) -> Box { + any.downcast_ref::() + .expect("expected MeasuredSource"); + let target = ILP::::new(1, vec![], vec![], ObjectiveSense::Minimize); + Box::new(ReductionAutoCast::>::new(target)) +} + +fn measured_edge( + reduce_fn: fn(&dyn Any) -> Box, + asymptotic_prediction: f64, +) -> ReductionEdgeData { + ReductionEdgeData { + overhead: ReductionOverhead::new(vec![( + "predicted_total", + Expr::Const(asymptotic_prediction), + )]), + reduce_fn: Some(reduce_fn), + reduce_aggregate_fn: None, + capabilities: EdgeCapabilities::witness_only(), + } +} + +// --------------------------------------------------------------------------- +// Verification 1: issue #788 known-answer check. +// --------------------------------------------------------------------------- + +/// The prism (triangular-prism) graph from issue #788: 6 vertices, 9 edges. +fn prism_hamiltonian_circuit() -> HamiltonianCircuit { + let prism = SimpleGraph::new( + 6, + vec![ + (0, 1), + (1, 2), + (2, 0), + (3, 4), + (4, 5), + (5, 3), + (0, 3), + (1, 4), + (2, 5), + ], + ); + HamiltonianCircuit::new(prism) +} + +/// #788: the measured Pareto search selects the path whose *measured* final ILP size is +/// smallest. +/// +/// The literal reduction chain quoted in issue #788 (HC → HP → ConsecutiveOnesSubmatrix → +/// ILP, total 60) no longer exists on the current reduction graph. The *current* measured +/// optimum is HC → LongestCircuit → ILP with a measured total of 232 +/// (num_constraints=127, num_vars=105); the next candidates are RuralPostman → ILP +/// (366) and TravelingSalesman → ILP (768). This test pins the measured optimum so +/// the selector is proven to rank by *measured* final size, not by step count or formula. +#[test] +fn test_hamiltoniancircuit_to_ilp_measured_optimum_788() { + let hc = prism_hamiltonian_circuit(); + let graph = ReductionGraph::new(); + let variant = ReductionGraph::variant_to_map(&[("graph", "SimpleGraph")]); + + let measured = graph + .find_measured_best_path_to_name( + "HamiltonianCircuit", + &variant, + "ILP", + ReductionMode::Witness, + &hc as &dyn Any, + 1_000, + ) + .expect("a measured witness path from HamiltonianCircuit to ILP"); + + // Measured final ILP size is the current-graph optimum. + assert_eq!( + measured.size.total(), + 232, + "measured optimum should be 232, got {:?}", + measured.size + ); + // Via LongestCircuit, to the bool ILP variant. + assert_eq!( + measured.path.type_names(), + vec!["HamiltonianCircuit", "LongestCircuit", "ILP"], + ); + + // The constructed chain is reusable: the final target is a genuine ILP. + use crate::models::algebraic::ILP; + let ilp = measured + .target_problem_any() + .downcast_ref::>() + .expect("final target is ILP"); + assert_eq!(ilp.num_vars, 105); +} + +// --------------------------------------------------------------------------- +// Verification 2: measured search does not discard equal-size concrete states. +// --------------------------------------------------------------------------- + +#[test] +fn test_measured_search_keeps_equal_size_structure_dependent_instances() { + let graph = ReductionGraph::from_test_edges( + &[ + "MeasuredSource", + "MeasuredBranchA", + "MeasuredBranchB", + "Satisfiability", + "ILP", + ], + &[ + ( + "MeasuredSource", + "MeasuredBranchA", + measured_edge(measured_source_to_a, 0.0), + ), + ( + "MeasuredSource", + "MeasuredBranchB", + measured_edge(measured_source_to_b, 0.0), + ), + ( + "MeasuredBranchA", + "Satisfiability", + measured_edge(measured_a_to_sat, 0.0), + ), + ( + "MeasuredBranchB", + "Satisfiability", + measured_edge(measured_b_to_sat, 0.0), + ), + ( + "Satisfiability", + "ILP", + measured_edge(measured_sat_to_structure_dependent_ilp, 0.0), + ), + ], + ); + let empty = BTreeMap::new(); + let source = MeasuredSource; + + let bad_sat = Satisfiability::new(1, vec![CNFClause::new(vec![1])]); + let good_sat = Satisfiability::new(1, vec![CNFClause::new(vec![-1])]); + assert_eq!( + ReductionGraph::compute_source_size("Satisfiability", &bad_sat), + ReductionGraph::compute_source_size("Satisfiability", &good_sat), + "the two structurally different hub instances must have identical measured sizes", + ); + + let measured = graph + .find_measured_best_path( + "MeasuredSource", + &empty, + "ILP", + &empty, + ReductionMode::Witness, + &source, + 1_000, + ) + .expect("the structure-dependent small continuation must survive"); + + assert_eq!( + measured.path.type_names(), + ["MeasuredSource", "MeasuredBranchB", "Satisfiability", "ILP",], + ); + assert_eq!(measured.size.total(), 1); +} + +#[test] +fn test_asymptotic_overhead_is_not_a_concrete_budget_guard() { + let graph = ReductionGraph::from_test_edges( + &["MeasuredSource", "ILP"], + &[( + "MeasuredSource", + "ILP", + measured_edge(measured_source_to_small_ilp, 1_000_000.0), + )], + ); + let empty = BTreeMap::new(); + let source = MeasuredSource; + + let measured = graph + .find_measured_best_path( + "MeasuredSource", + &empty, + "ILP", + &empty, + ReductionMode::Witness, + &source, + 1, + ) + .expect("a loose asymptotic expression must not prune an actually in-budget target"); + + assert_eq!(measured.size.total(), 1); +} + +// --------------------------------------------------------------------------- +// Verification 4: negative control on a hand-built diamond. +// --------------------------------------------------------------------------- + +/// A test label whose objective is the *final* measured size `s`, while carrying a +/// separate accumulated step cost `c`. Dominance is componentwise Pareto over `(c, s)`, +/// so two labels that trade off `c` against `s` are incomparable and both survive — the +/// exact structure a scalar Dijkstra collapses (keeping only the min-`c` label, and thus +/// its `s`). +#[derive(Clone)] +struct DiamondLabel { + /// Accumulated step cost. + c: f64, + /// Current (path-dependent) measured size. + s: f64, +} + +impl DiamondLabel { + fn ctx(&self) -> ProblemSize { + ProblemSize::new(vec![("s", self.s.round().max(0.0) as usize)]) + } +} + +impl PathLabel for DiamondLabel { + fn extend(&self, edge: &ReductionEdge) -> Option { + let ctx = self.ctx(); + let add_c = edge.overhead.get("c").map(|e| e.eval(&ctx)).unwrap_or(0.0); + let new_s = edge + .overhead + .get("s") + .map(|e| e.eval(&ctx)) + .unwrap_or(self.s); + Some(DiamondLabel { + c: self.c + add_c, + s: new_s, + }) + } + + fn dominates(&self, other: &Self) -> bool { + self.c <= other.c && self.s <= other.s + } + + fn cost(&self) -> f64 { + self.s + } +} + +fn diamond_edge(c: f64, s: Expr) -> ReductionEdgeData { + ReductionEdgeData { + overhead: ReductionOverhead::new(vec![("c", Expr::Const(c)), ("s", s)]), + reduce_fn: None, + reduce_aggregate_fn: None, + capabilities: EdgeCapabilities::witness_only(), + } +} + +/// Negative control: P1 (S→M→T) has the lower first-edge cost but a larger measured +/// intermediate size at M; P2 (S→P→M→T) has a higher first-edge cost but a strictly +/// smaller final measured size. A scalar-cost path selection (`find_cheapest_path` over +/// the additive step cost) commits to P1's prefix at M and returns P1; the measured +/// Pareto search keeps both routes into M (they are incomparable) and returns P2. +#[test] +fn test_negative_control_diamond_pareto_beats_scalar() { + let empty = std::collections::BTreeMap::new(); + let graph = ReductionGraph::from_test_edges( + &["S", "M", "P", "T"], + &[ + // S -> M: cheap first edge (c=1), large intermediate size (s=100). + ("S", "M", diamond_edge(1.0, Expr::Const(100.0))), + // S -> P: pricier first edge (c=2), small size (s=5). + ("S", "P", diamond_edge(2.0, Expr::Const(5.0))), + // P -> M: small size (s=6). + ("P", "M", diamond_edge(1.0, Expr::Const(6.0))), + // M -> T: identity on size (final size = size at M). + ("M", "T", diamond_edge(1.0, Expr::Var("s"))), + ], + ); + + // (a) Scalar-cost selection (minimize additive step cost `c`) commits to P1. + let scalar = graph + .find_cheapest_path( + "S", + &empty, + "T", + &empty, + &ProblemSize::new(vec![]), + &CustomCost(|oh: &ReductionOverhead, sz: &ProblemSize| { + oh.get("c").map(|e| e.eval(sz)).unwrap_or(0.0) + }), + ) + .expect("scalar path S -> T"); + assert_eq!( + scalar.type_names(), + vec!["S", "M", "T"], + "scalar cost selection should commit to the cheap-prefix P1" + ); + + // (b) The measured Pareto search returns P2 (strictly smaller final size). + let initial = DiamondLabel { c: 0.0, s: 0.0 }; + let front = graph.pareto_search_by_name( + "S", + &empty, + "T", + &empty, + ReductionMode::Witness, + initial, + false, + ); + assert!(!front.is_empty(), "front should reach T"); + let (best_path, best_label) = &front[0]; + assert_eq!( + best_path.type_names(), + vec!["S", "P", "M", "T"], + "Pareto search should return the better-final-size P2" + ); + assert_eq!(best_label.cost(), 6.0, "P2's final measured size is 6"); +} + +/// The `exhaustive` flag disables only the heuristic componentwise-dominance guard; the +/// front still contains the true optimum. On the diamond, both routes into M survive +/// regardless, so the answer is unchanged. +#[test] +fn test_diamond_exhaustive_matches_pruned() { + let empty = std::collections::BTreeMap::new(); + let graph = ReductionGraph::from_test_edges( + &["S", "M", "P", "T"], + &[ + ("S", "M", diamond_edge(1.0, Expr::Const(100.0))), + ("S", "P", diamond_edge(2.0, Expr::Const(5.0))), + ("P", "M", diamond_edge(1.0, Expr::Const(6.0))), + ("M", "T", diamond_edge(1.0, Expr::Var("s"))), + ], + ); + let front = graph.pareto_search_by_name( + "S", + &empty, + "T", + &empty, + ReductionMode::Witness, + DiamondLabel { c: 0.0, s: 0.0 }, + true, + ); + assert_eq!(front[0].0.type_names(), vec!["S", "P", "M", "T"]); + assert_eq!(front[0].1.cost(), 6.0); +} + +// --------------------------------------------------------------------------- +// GrowthLabel (asymptotic, instance-free) domain — issue #1080 / design M3/F3a. +// --------------------------------------------------------------------------- + +/// A power `Var(v)^k`. +fn powk(v: &'static str, k: f64) -> Expr { + Expr::pow(Expr::Var(v), Expr::Const(k)) +} + +/// A test edge carrying only a symbolic overhead (target field → Expr over the +/// current node's fields), no executable reduction. +fn growth_edge(fields: Vec<(&'static str, Expr)>) -> ReductionEdgeData { + ReductionEdgeData { + overhead: ReductionOverhead::new(fields), + reduce_fn: None, + reduce_aggregate_fn: None, + capabilities: EdgeCapabilities::witness_only(), + } +} + +/// The rendered Big-O string for one field of a growth label (or `"?"` for +/// `Unknown`), for compact assertions. +fn field_big_o(label: &GrowthLabel, field: &str) -> String { + match label.fields().get(field) { + Some(g) => match g.to_expr() { + Some(e) => e.to_string(), + None => "?".to_string(), + }, + None => "".to_string(), + } +} + +/// `extend` substitutes the current label's growth into an edge's overhead and +/// reduces in the growth domain, yielding the target field's growth in source vars. +#[test] +fn test_growth_label_extend_composes_overhead() { + // Source S has fields n, m; edge maps a = n^2, b = m (in the source's variables). + let edge_data = growth_edge(vec![("a", powk("n", 2.0)), ("b", Expr::Var("m"))]); + let target_variant = BTreeMap::new(); + let redge = ReductionEdge { + overhead: &edge_data.overhead, + reduce_fn: None, + capabilities: EdgeCapabilities::witness_only(), + target_name: "Target", + target_variant: &target_variant, + }; + + let initial = GrowthLabel::source(&["n", "m"]); + let next = initial + .extend(&redge) + .expect("asymptotic extend never prunes"); + assert_eq!(field_big_o(&next, "a"), "n^2"); + assert_eq!(field_big_o(&next, "b"), "m"); + + // A second hop composes: c = a * b substitutes a→n^2, b→m ⇒ n^2 * m. + let edge2 = growth_edge(vec![("c", Expr::Var("a") * Expr::Var("b"))]); + let redge2 = ReductionEdge { + overhead: &edge2.overhead, + reduce_fn: None, + capabilities: EdgeCapabilities::witness_only(), + target_name: "Target2", + target_variant: &target_variant, + }; + let composed = next.extend(&redge2).expect("extend"); + assert_eq!(field_big_o(&composed, "c"), "m * n^2"); +} + +/// An overhead field that depends on an `Unknown`-growth current field stays +/// `Unknown` — the bound is never fabricated. +#[test] +fn test_growth_label_propagates_unknown() { + // Build a label whose field `x` is Unknown (factorial growth). + let mut fields = BTreeMap::new(); + fields.insert( + "x", + Growth::from_expr(&Expr::Factorial(Box::new(Expr::Var("n")))), + ); + fields.insert("y", Growth::from_expr(&Expr::Var("n"))); + let label = GrowthLabel::from_fields(fields); + assert!(matches!(label.fields().get("x"), Some(Growth::Unknown))); + + // out1 uses x (Unknown) → Unknown; out2 uses only y → bounded. + let edge = growth_edge(vec![ + ("out1", Expr::Var("x") * Expr::Var("y")), + ("out2", powk("y", 2.0)), + ]); + let tv = BTreeMap::new(); + let redge = ReductionEdge { + overhead: &edge.overhead, + reduce_fn: None, + capabilities: EdgeCapabilities::witness_only(), + target_name: "T", + target_variant: &tv, + }; + let next = label.extend(&redge).expect("extend"); + assert_eq!(field_big_o(&next, "out1"), "?"); + assert_eq!(field_big_o(&next, "out2"), "n^2"); +} + +/// A label with an `Unknown` field is dominated by any fully-known label, and never +/// dominates one — undecidable paths rank last. +#[test] +fn test_growth_label_unknown_ranks_last() { + let known = GrowthLabel::from_fields({ + let mut m = BTreeMap::new(); + m.insert("a", Growth::from_expr(&powk("n", 2.0))); + m.insert("b", Growth::from_expr(&Expr::Var("m"))); + m + }); + let with_unknown = GrowthLabel::from_fields({ + let mut m = BTreeMap::new(); + m.insert("a", Growth::from_expr(&powk("n", 2.0))); + m.insert("b", Growth::Unknown); + m + }); + // Known is strictly better on field b (n^0? no: bounded vs Unknown) ⇒ known dominates. + assert!(known.dominates(&with_unknown)); + assert!(!with_unknown.dominates(&known)); +} + +/// Componentwise search-sense dominance: `self` dominates `other` iff it grows no +/// faster on every field and strictly slower on at least one. +#[test] +fn test_growth_label_dominance_partial_order() { + let a = GrowthLabel::from_fields({ + let mut m = BTreeMap::new(); + m.insert("v", Growth::from_expr(&Expr::Var("n"))); // n + m.insert("e", Growth::from_expr(&Expr::Var("m"))); // m + m + }); + let b = GrowthLabel::from_fields({ + let mut m = BTreeMap::new(); + m.insert("v", Growth::from_expr(&powk("n", 2.0))); // n^2 + m.insert("e", Growth::from_expr(&Expr::Var("m"))); // m + m + }); + // a (n, m) grows slower in v, equal in e ⇒ a dominates b; b does not dominate a. + assert!(a.dominates(&b)); + assert!(!b.dominates(&a)); + // Reflexivity is *not* strict dominance: equal labels do not dominate each other. + assert!(!a.dominates(&a.clone())); + + // Incomparable pair: one better in v, the other better in e. + let c = GrowthLabel::from_fields({ + let mut m = BTreeMap::new(); + m.insert("v", Growth::from_expr(&powk("n", 2.0))); // n^2 + m.insert("e", Growth::from_expr(&Expr::Var("m"))); // m + m + }); + let d = GrowthLabel::from_fields({ + let mut m = BTreeMap::new(); + m.insert("v", Growth::from_expr(&Expr::Var("n"))); // n + m.insert("e", Growth::from_expr(&powk("m", 2.0))); // m^2 + m + }); + assert!(!c.dominates(&d)); + assert!(!d.dominates(&c)); +} + +/// **Negative control (issue #1080):** two S→T paths whose composed growths are +/// incomparable — path A costs `O(n^2)` in `vertices` / `O(m)` in `edges`, path B +/// costs `O(n)` / `O(m^2)` — must *both* appear in the asymptotic Pareto front. An +/// implementation that scalarizes or keeps a single representative fails this. +#[test] +fn test_growth_negative_control_incomparable_front() { + let empty = BTreeMap::new(); + let graph = ReductionGraph::from_test_edges( + &["S", "A", "B", "T"], + &[ + // Both prefixes just carry the source fields n, m through unchanged. + ( + "S", + "A", + growth_edge(vec![("n", Expr::Var("n")), ("m", Expr::Var("m"))]), + ), + ( + "S", + "B", + growth_edge(vec![("n", Expr::Var("n")), ("m", Expr::Var("m"))]), + ), + // Path A: vertices = n^2, edges = m. + ( + "A", + "T", + growth_edge(vec![ + ("vertices", powk("n", 2.0)), + ("edges", Expr::Var("m")), + ]), + ), + // Path B: vertices = n, edges = m^2. + ( + "B", + "T", + growth_edge(vec![ + ("vertices", Expr::Var("n")), + ("edges", powk("m", 2.0)), + ]), + ), + ], + ); + + let initial = GrowthLabel::source(&["n", "m"]); + let front = graph.pareto_search_by_name( + "S", + &empty, + "T", + &empty, + ReductionMode::Witness, + initial, + false, + ); + + // The front must contain BOTH incomparable paths — not one representative. + assert_eq!( + front.len(), + 2, + "front should keep both incomparable paths, got {:?}", + front + .iter() + .map(|(p, _)| p.type_names()) + .collect::>() + ); + let mut seen: Vec<(String, String)> = front + .iter() + .map(|(p, label)| { + ( + p.type_names().join("→"), + format!( + "v={} e={}", + field_big_o(label, "vertices"), + field_big_o(label, "edges") + ), + ) + }) + .collect(); + seen.sort(); + assert_eq!( + seen, + vec![ + ("S→A→T".to_string(), "v=n^2 e=m".to_string()), + ("S→B→T".to_string(), "v=n e=m^2".to_string()), + ], + ); +} + +// Completeness under ASYMMETRIC magnitudes: the two incomparable paths have +// different scalar `cost` summaries (A: n^2 + m ⇒ magnitude 3; B: n + m^3 ⇒ +// magnitude 4). A scalar branch-and-bound (were the kernel to use one) would let the +// cheaper path A complete first and then prune B (cost 4 ≥ 3), silently dropping a +// Pareto-optimal path. This is the case the equal-magnitude negative control above +// does NOT catch; it passes because the kernel prunes by exact dominance only, never +// by the scalar `cost`. +#[test] +fn test_growth_asymmetric_incomparable_front_complete() { + let empty = BTreeMap::new(); + let graph = ReductionGraph::from_test_edges( + &["S", "A", "B", "T"], + &[ + ( + "S", + "A", + growth_edge(vec![("n", Expr::Var("n")), ("m", Expr::Var("m"))]), + ), + ( + "S", + "B", + growth_edge(vec![("n", Expr::Var("n")), ("m", Expr::Var("m"))]), + ), + // Path A: vertices = n^2, edges = m (magnitude 2 + 1 = 3). + ( + "A", + "T", + growth_edge(vec![ + ("vertices", powk("n", 2.0)), + ("edges", Expr::Var("m")), + ]), + ), + // Path B: vertices = n, edges = m^3 (magnitude 1 + 3 = 4). + ( + "B", + "T", + growth_edge(vec![ + ("vertices", Expr::Var("n")), + ("edges", powk("m", 3.0)), + ]), + ), + ], + ); + + let front = graph.pareto_search_by_name( + "S", + &empty, + "T", + &empty, + ReductionMode::Witness, + GrowthLabel::source(&["n", "m"]), + false, + ); + + let mut seen: Vec<(String, String)> = front + .iter() + .map(|(p, label)| { + ( + p.type_names().join("→"), + format!( + "v={} e={}", + field_big_o(label, "vertices"), + field_big_o(label, "edges") + ), + ) + }) + .collect(); + seen.sort(); + assert_eq!( + seen, + vec![ + ("S→A→T".to_string(), "v=n^2 e=m".to_string()), + ("S→B→T".to_string(), "v=n e=m^3".to_string()), + ], + "both incomparable paths must survive despite different scalar magnitudes", + ); +} + +/// Isotonicity of `extend` (design invariant): if `A` dominates `B`, then +/// `extend(A, e)` dominates `extend(B, e)` for the same edge — the correctness +/// condition for the kernel's dominance pruning. +#[test] +fn test_growth_label_extend_isotone() { + // A = (n, m) dominates B = (n^2, m^2) componentwise. + let a = GrowthLabel::source(&["n", "m"]); + let b = GrowthLabel::from_fields({ + let mut mm = BTreeMap::new(); + mm.insert("n", Growth::from_expr(&powk("n", 2.0))); + mm.insert("m", Growth::from_expr(&powk("m", 2.0))); + mm + }); + assert!(a.dominates(&b)); + + let tv = BTreeMap::new(); + // A monotone overhead in both fields. + for overhead in [ + growth_edge(vec![("x", Expr::Var("n") * Expr::Var("m"))]), + growth_edge(vec![("x", powk("n", 3.0)), ("y", Expr::Var("m"))]), + ] { + let redge = ReductionEdge { + overhead: &overhead.overhead, + reduce_fn: None, + capabilities: EdgeCapabilities::witness_only(), + target_name: "T", + target_variant: &tv, + }; + let ea = a.extend(&redge).unwrap(); + let eb = b.extend(&redge).unwrap(); + // A ⪰ B ⇒ extend(A) ⪰ extend(B) (dominates-or-equal). Equality is possible + // when the overhead collapses the difference, so accept dominate-or-equal. + assert!( + ea.dominates(&eb) || ea == eb, + "isotonicity violated: {ea:?} vs {eb:?}" + ); + } +} + +/// `asymptotic_front` reports **one representative per distinct growth vector**, not +/// one per route. On the real graph, `MinimumVertexCover → ILP` has dozens of +/// syntactically distinct reduction chains that compose to only a handful of Big-O +/// profiles; the front must (a) contain no two entries with identical growth vectors +/// and (b) collapse to that small handful — while the raw kernel front (same search, +/// no dedup) still holds the many redundant routes. +#[test] +fn test_asymptotic_front_dedups_by_growth_vector() { + let graph = ReductionGraph::new(); + let src_v = graph + .default_variant_for("MinimumVertexCover") + .or_else(|| graph.variants_for("MinimumVertexCover").into_iter().next()) + .expect("MinimumVertexCover registered"); + let dst_v = graph + .default_variant_for("ILP") + .or_else(|| graph.variants_for("ILP").into_iter().next()) + .expect("ILP registered"); + + let front = graph.asymptotic_front( + "MinimumVertexCover", + &src_v, + "ILP", + &dst_v, + ReductionMode::Witness, + ); + assert!(!front.is_empty(), "MVC -> ILP must have a path"); + + // (a) No two front entries share a growth vector (GrowthLabel PartialEq). + for i in 0..front.len() { + for j in (i + 1)..front.len() { + assert!( + front[i].1 != front[j].1, + "duplicate growth vector in front:\n {}\n {}", + front[i].0.type_names().join("→"), + front[j].0.type_names().join("→"), + ); + } + } + // (b) A proper Pareto front is a small handful, not the dozens of redundant routes. + assert!( + (1..=4).contains(&front.len()), + "expected 1..=4 distinct growth vectors, got {}", + front.len() + ); + + // The dedup genuinely collapsed routes: the raw kernel front (same search, no + // dedup) is strictly larger and does contain repeated growth vectors. + let src_fields = graph.size_field_names("MinimumVertexCover"); + let raw = graph.pareto_search_by_name( + "MinimumVertexCover", + &src_v, + "ILP", + &dst_v, + ReductionMode::Witness, + GrowthLabel::source(&src_fields), + false, + ); + assert!( + raw.len() > front.len(), + "dedup should collapse redundant routes: raw {} vs deduped {}", + raw.len(), + front.len() + ); +} + +/// A composed front label must express every size field's growth purely in the +/// **source problem's** own size variables — never in a downstream getter alias or an +/// intermediate node's field name. +/// +/// Regression for the `MinimumFeedbackVertexSet → ILP` bug: the `ILP → ILP` +/// binary-encoding cast declared its overhead as `num_vars = "31 * num_variables"`, +/// referencing the getter *alias* `num_variables()` instead of ILP's size-field *name* +/// `num_vars`. Instance mode and raw-overhead rendering both resolve the getter, so the +/// mistake was invisible there — but growth composition threads field *names*, so the +/// alias was unmapped and leaked through as `num_vars = O(num_variables)` instead of +/// the correct `O(num_vertices)`. +#[test] +fn test_asymptotic_front_uses_only_source_variables_mfvs_ilp() { + let graph = ReductionGraph::new(); + let src_v = graph + .default_variant_for("MinimumFeedbackVertexSet") + .or_else(|| { + graph + .variants_for("MinimumFeedbackVertexSet") + .into_iter() + .next() + }) + .expect("MinimumFeedbackVertexSet registered"); + let dst_v = graph + .default_variant_for("ILP") + .or_else(|| graph.variants_for("ILP").into_iter().next()) + .expect("ILP registered"); + + let front = graph.asymptotic_front( + "MinimumFeedbackVertexSet", + &src_v, + "ILP", + &dst_v, + ReductionMode::Witness, + ); + + // The direct route (MFVS → ILP/i32 → ILP/bool; the ILP variants collapse in the + // deduplicated node-name view) is the one exercised by the fixed cast. + let (_, label) = front + .iter() + .find(|(p, _)| p.type_names() == ["MinimumFeedbackVertexSet", "ILP"]) + .expect("direct MinimumFeedbackVertexSet -> ILP path"); + + // The size fields of MinimumFeedbackVertexSet — the only variables any composed + // growth is allowed to mention. + let allowed = ["num_arcs", "num_vertices"]; + for (field, growth) in label.fields() { + let expr = growth + .to_expr() + .unwrap_or_else(|| panic!("field {field} should have a bounded growth")); + for var in expr.variables() { + assert!( + allowed.contains(&var), + "field `{field}` growth O({expr}) references `{var}`, which is not a \ + MinimumFeedbackVertexSet source variable {allowed:?}", + ); + } + } + + // The previously-buggy field, pinned to the correct source-variable Big-O. + let num_vars = label + .fields() + .get("num_vars") + .expect("ILP has a num_vars size field"); + assert_eq!( + num_vars.to_expr().unwrap().to_string(), + "num_vertices", + "ILP num_vars must compose to O(num_vertices), not the getter alias num_variables" + ); +} + +// --------------------------------------------------------------------------- +// Fix A: the kernel prunes by dominance only — never (unsound) branch-and-bound. +// --------------------------------------------------------------------------- + +/// A test label whose `cost` is the label's current absolute value — a value a late edge +/// can *shrink* below an already-completed route's final value. It verifies that the +/// generic kernel does not silently add scalar branch-and-bound. +#[derive(Clone)] +struct ShrinkLabel { + v: f64, +} + +impl PathLabel for ShrinkLabel { + fn extend(&self, edge: &ReductionEdge) -> Option { + // The edge sets a new absolute value (`v`), which may be smaller than the current. + let z = ProblemSize::new(vec![]); + let v = edge.overhead.get("v").map(|e| e.eval(&z)).unwrap_or(self.v); + Some(ShrinkLabel { v }) + } + + fn dominates(&self, other: &Self) -> bool { + self.v <= other.v + } + + fn cost(&self) -> f64 { + self.v + } +} + +/// Kernel regression for Fix A: a route that *shrinks late* (its intermediate cost 100 is +/// higher than a rival route that completes early at 50, but a final edge drops it to 10) +/// must survive to the front. A kernel that applied branch-and-bound would prune the +/// intermediate node (100 ≥ best-so-far 50) and silently drop the true optimum. Because +/// the kernel prunes by dominance only, the shrink-late route reaches the front even under +/// `exhaustive = true` (which disables only the dominance guard). +#[test] +fn test_kernel_keeps_shrink_late_route_dominance_only() { + let empty = std::collections::BTreeMap::new(); + let graph = ReductionGraph::from_test_edges( + &["S", "A", "T"], + &[ + // S -> T: completes early with final value 50. + ("S", "T", growth_edge(vec![("v", Expr::Const(50.0))])), + // S -> A: intermediate value 100 (would trip a B&B bound of 50). + ("S", "A", growth_edge(vec![("v", Expr::Const(100.0))])), + // A -> T: shrinks the value to 10 (globally best). + ("A", "T", growth_edge(vec![("v", Expr::Const(10.0))])), + ], + ); + + let front = graph.pareto_search_by_name( + "S", + &empty, + "T", + &empty, + ReductionMode::Witness, + ShrinkLabel { v: 0.0 }, + true, + ); + + // The shrink-late route S -> A -> T (final value 10) must be present in the front. + let shrink_late = front + .iter() + .find(|(p, _)| p.type_names() == ["S", "A", "T"]) + .expect("shrink-late route S -> A -> T must survive without branch-and-bound"); + assert_eq!( + shrink_late.1.cost(), + 10.0, + "the shrink-late route finishes at the global optimum value 10" + ); + // The kernel's best (lowest cost) front element is that shrink-late route. + assert_eq!(front[0].0.type_names(), ["S", "A", "T"]); + assert_eq!(front[0].1.cost(), 10.0); +} + +// --------------------------------------------------------------------------- +// Fix B: CostLabel dominance is componentwise over (cost, size). +// --------------------------------------------------------------------------- + +/// Fix B regression: an edge cost that DEPENDS on the carried size makes a cheaper-so-far +/// prefix with a *larger* intermediate size a trap — a scalar `cost <= other.cost` +/// dominance would evict the costlier-but-smaller prefix whose continuation is globally +/// cheapest. With componentwise `(cost, size)` dominance both prefixes survive at the hub +/// and `find_cheapest_path` returns the globally optimal route. +#[test] +fn test_cost_label_path_dependent_dominance() { + let empty = std::collections::BTreeMap::new(); + // Edges carry `c` (base edge cost), `wf` (weight on the size-dependent term) and `w` + // (the tracked size field). The cost function is `c + wf * current_w`, so the M -> T + // edge's cost is exactly the size `w` accumulated at M. + let graph = ReductionGraph::from_test_edges( + &["S", "M", "P", "T"], + &[ + // S -> M: cheap prefix (c = 1) but produces a LARGE intermediate size w = 100. + ( + "S", + "M", + growth_edge(vec![ + ("c", Expr::Const(1.0)), + ("wf", Expr::Const(0.0)), + ("w", Expr::Const(100.0)), + ]), + ), + // S -> P: pricier prefix (c = 3) but a SMALL size w = 1. + ( + "S", + "P", + growth_edge(vec![ + ("c", Expr::Const(3.0)), + ("wf", Expr::Const(0.0)), + ("w", Expr::Const(1.0)), + ]), + ), + // P -> M: cheap (c = 1), keeps the small size w = 1. + ( + "P", + "M", + growth_edge(vec![ + ("c", Expr::Const(1.0)), + ("wf", Expr::Const(0.0)), + ("w", Expr::Const(1.0)), + ]), + ), + // M -> T: cost = current w (wf = 1, c = 0); identity on size. + ( + "M", + "T", + growth_edge(vec![ + ("c", Expr::Const(0.0)), + ("wf", Expr::Const(1.0)), + ("w", Expr::Var("w")), + ]), + ), + ], + ); + + // Cost function: c + wf * current_w. Depends on the carried size, so the two prefixes + // into M are incomparable and must both be kept. + let cost_fn = CustomCost(|oh: &ReductionOverhead, sz: &ProblemSize| { + let c = oh.get("c").map(|e| e.eval(sz)).unwrap_or(0.0); + let wf = oh.get("wf").map(|e| e.eval(sz)).unwrap_or(0.0); + c + wf * sz.get("w").unwrap_or(0) as f64 + }); + + let best = graph + .find_cheapest_path( + "S", + &empty, + "T", + &empty, + &ProblemSize::new(vec![("w", 0)]), + &cost_fn, + ) + .expect("cheapest path S -> T"); + + // Globally cheapest: S -> P -> M -> T (total 3 + 1 + 1 = 5), NOT the cheap-prefix trap + // S -> M -> T (total 1 + 100 = 101). A scalar-dominance CostLabel would evict the + // small-w prefix at M and return the S -> M -> T trap. + assert_eq!( + best.type_names(), + vec!["S", "P", "M", "T"], + "componentwise (cost, size) dominance must keep the globally optimal small-w prefix" + ); +} + +// --------------------------------------------------------------------------- +// Fix C: GrowthLabel taints target fields referencing intermediate-only variables. +// --------------------------------------------------------------------------- + +/// Fix C regression: an overhead output expression that references a variable ABSENT from +/// the current label (an intermediate-only field, e.g. `tseitin_*`, `num_encoding_bits`) +/// must taint its target field to `Growth::Unknown` — it must NOT pass through +/// `substitute` verbatim and surface as a fake source variable in the final bound. +#[test] +fn test_growth_label_taints_absent_variable() { + // The label knows only the source field `n`. + let label = GrowthLabel::source(&["n"]); + // Edge output: `bounded` depends only on `n`; `leaky` references `tseitin`, which is + // absent from the label (an intermediate-only construction variable). + let edge = growth_edge(vec![ + ("bounded", Expr::Var("n")), + ("leaky", Expr::Var("n") * Expr::Var("tseitin")), + ]); + let tv = BTreeMap::new(); + let redge = ReductionEdge { + overhead: &edge.overhead, + reduce_fn: None, + capabilities: EdgeCapabilities::witness_only(), + target_name: "T", + target_variant: &tv, + }; + let next = label.extend(&redge).expect("extend"); + + // Depends only on a mapped source variable ⇒ stays bounded. + assert_eq!(field_big_o(&next, "bounded"), "n"); + // References an unmapped, intermediate-only variable ⇒ tainted to Unknown, never + // leaked as `O(n * tseitin)`. + assert!( + matches!(next.fields().get("leaky"), Some(Growth::Unknown)), + "a target field referencing an absent variable must become Unknown, got {:?}", + next.fields().get("leaky") + ); +} + +// --------------------------------------------------------------------------- +// Fix D: the arena frees evicted labels (bag cap bounds retained instance memory). +// --------------------------------------------------------------------------- + +thread_local! { + /// Live token instances on this thread. + static TOK_LIVE: Cell = const { Cell::new(0) }; + /// Peak live token instances observed. + static TOK_PEAK: Cell = const { Cell::new(0) }; + /// Total token instances ever created. + static TOK_CREATED: Cell = const { Cell::new(0) }; +} + +/// A drop-tracking token. Each `new()` is a distinct live instance; `Drop` frees it. Held +/// behind `Rc` inside a label, so cloning a label shares the token. If the arena pinned +/// evicted labels, their tokens would stay live until the search ended, so `TOK_PEAK` +/// would reach `TOK_CREATED`. +struct DropToken; + +impl DropToken { + fn new() -> Self { + let live = TOK_LIVE.with(|c| { + let v = c.get() + 1; + c.set(v); + v + }); + TOK_PEAK.with(|p| { + if live > p.get() { + p.set(live); + } + }); + TOK_CREATED.with(|c| c.set(c.get() + 1)); + DropToken + } +} + +impl Drop for DropToken { + fn drop(&mut self) { + TOK_LIVE.with(|c| c.set(c.get() - 1)); + } +} + +/// A label carrying an `Rc` and a two-component `(c, s)` value. The engineered +/// `(c, s)` pairs are pairwise incomparable, so no label evicts another by dominance and +/// the per-node bag grows until the cap truncates it — exercising the truncation free path. +#[derive(Clone)] +struct TokenLabel { + c: f64, + s: f64, + _tok: Rc, +} + +impl PathLabel for TokenLabel { + fn extend(&self, edge: &ReductionEdge) -> Option { + let z = ProblemSize::new(vec![]); + let c = edge.overhead.get("c").map(|e| e.eval(&z)).unwrap_or(self.c); + let s = edge.overhead.get("s").map(|e| e.eval(&z)).unwrap_or(self.s); + Some(TokenLabel { + c, + s, + _tok: Rc::new(DropToken::new()), + }) + } + + fn dominates(&self, other: &Self) -> bool { + self.c <= other.c && self.s <= other.s + } + + fn cost(&self) -> f64 { + self.c + } +} + +/// Fix D regression: drive the kernel on a graph that generates far more labels at one hub +/// than `BAG_CAP`, all incomparable so the bag truncates repeatedly. Because evicted / +/// truncated arena entries free their labels immediately, the *peak* number of live +/// `DropToken` instances stays well below the *total* ever created. If the arena pinned +/// evicted labels (the bug), peak would equal total. +#[test] +fn test_arena_frees_evicted_labels_bounds_live_memory() { + TOK_LIVE.with(|c| c.set(0)); + TOK_PEAK.with(|c| c.set(0)); + TOK_CREATED.with(|c| c.set(0)); + + // One hub M fed by N ≫ BAG_CAP parallel S -> M edges with pairwise-incomparable + // (c = i+1, s = N-i) labels, then M -> T (identity). The M bag truncates repeatedly. + let n: usize = 200; + let mut edges: Vec<(&'static str, &'static str, ReductionEdgeData)> = Vec::new(); + // Leak small &'static str-free constants via Expr::Const (no string needed for values). + for i in 0..n { + edges.push(( + "S", + "M", + growth_edge(vec![ + ("c", Expr::Const((i + 1) as f64)), + ("s", Expr::Const((n - i) as f64)), + ]), + )); + } + edges.push(( + "M", + "T", + growth_edge(vec![("c", Expr::Var("c")), ("s", Expr::Var("s"))]), + )); + let graph = ReductionGraph::from_test_edges(&["S", "M", "T"], &edges); + + let empty = std::collections::BTreeMap::new(); + let initial = TokenLabel { + c: 0.0, + s: 0.0, + _tok: Rc::new(DropToken::new()), + }; + let front = graph.pareto_search_by_name( + "S", + &empty, + "T", + &empty, + ReductionMode::Witness, + initial, + false, + ); + // Sanity: the search reached T. + assert!(!front.is_empty(), "front should reach T"); + + let created = TOK_CREATED.with(|c| c.get()); + let peak = TOK_PEAK.with(|c| c.get()); + // Many labels were created (≥ the N hub edges). + assert!( + created >= n as i64, + "expected many token instances created, got {created}" + ); + // Eviction frees labels: peak live is strictly below total created. With the bug + // (arena pins evicted labels) peak would equal created; the margin here is large + // (peak is bounded by ~BAG_CAP per live node, created scales with N) so this is not + // flaky. + assert!( + peak < created, + "arena must free evicted labels: peak {peak} should be < created {created}" + ); + + // The retained tokens are bounded by the live bag entries, not by N. Concretely, far + // fewer than the total are still live once the search completes. + drop(front); + let live_after = TOK_LIVE.with(|c| c.get()); + assert!( + live_after < created, + "retained tokens {live_after} must be bounded well below total {created}" + ); +}