Conversation
- Use the new interface of AbstractOperators.jl v0.4 - Add new parent package (OperatorCore) and subpackages (FFTWOperators, DSPOperators) of AbstractOperators - Use TermSet instead of tuple of Terms - Implement parsing for LeastSquaresTerm - Add name field for Variable - rename back L1,2-norm from mixednorm to norm - separate Project.toml for test
Adds ArrayPartition dispatch for SeparableSum prox/gradient, fixes slicing-mask helper calls to go through AbstractOperators, marks Aqua persistent_tasks as broken (Julia 1.12 HPC false positive), and updates test/doc Project.toml dependency bounds. Excludes AGENTS.md and Manifest.toml, tracked separately. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01478URHkYh8YPDHBLznsAR7
…perators.jl ProximalOperators.jl gained a native RecursiveArrayToolsExt extension for SeparableSum + ArrayPartition, which now conflicts (method overwriting) with the equivalent dispatch this package defined locally, breaking precompilation. Also un-marks Aqua persistent_tasks as broken now that it passes, and drops the now-unused SeparableSum piracy allowlist entry. Ignore generated LocalCoverage output (*.cov, coverage/, coverage_html/). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01478URHkYh8YPDHBLznsAR7
… on nonconvex problems ProximalOperators.jl correctly marks IndBallL1 as is_proximable=false (no guaranteed exact prox), so tests asserting it should be proximable, and solves relying on that assumption, were wrong. Switch the affected tests to IndBallL2, which is genuinely proximable and exercises the same multi-variable/constrained-solve code paths. Also fixes two masked failures uncovered once the above was corrected: test_build_minimize.jl solved nonconvex problems (nonlinear sigmoid composition, Rosenbrock function) with FastForwardBackward, whose stated assumption requires convexity of f. Swapped for ZeroFPR, which handles nonconvex smooth f. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01478URHkYh8YPDHBLznsAR7
Add coverage visibility and doctest execution in CI, without gating on either. - codecov.yml in informational mode (records the 68% baseline; never fails CI) - ci.yml: run tests with coverage, process to lcov.info, upload via codecov-action; bump actions to v2/v4 and add julia-actions/cache - docs/make.jl: doctest = true so doc code blocks run as tests in CI Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01478URHkYh8YPDHBLznsAR7
Fix ten confirmed defects, each covered by a value-asserting regression test in
test/test_phase1_regressions.jl:
- parse.jl 1.1/1.2/1.3: stop double-counting displacement/λ when a sum of smooth
terms contains a nonlinear composition, and in the OperatorTerm/InfConv TermSet
paths (carry displacement once in the affine operator); return the checked op in
the InfConv func₂ branch.
- parse.jl 1.4 (LeastSquaresTerm): scale the residual by √(term.lambda*f.lambda) so
the data term is weighted correctly relative to the SquaredL2 regularizer; fold in
f.lambda; fix the b sign (b = -displacement, since A*x - b stores displacement -b);
reject LeastSquares functions whose embedded A/b this path cannot read.
- sqrNormL2WithNormalOp.jl 1.5: weighted gradient is Aᴴ·diag(λ)·A·x (weights in the
codomain); is_strongly_convex requires full column rank and its dispatch now matches
the fully-parameterized type.
- build_solve.jl 1.6/1.7: accept AbstractVector{<:IterativeAlgorithm}; guard the Tuple
minimizer in the auto-select path; factor the three solve bodies into _run_solver.
- term.jl 1.8: a * TermSet stays a TermSet; scalar-* on Term preserves repr.
- addition_tricky_part.jl 1.9: UnregularIndex length is prod(max), not sum.
- 1.10: value_and_gradient stays a generic bridge (documented) — it must accept any
foreign smooth function the package composes, so it cannot be narrowed to owned
types; kept in the Aqua treat_as_own allowlist with a rationale comment.
Also replaces a few term.A.L field accesses in parse.jl with affine(term).
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01478URHkYh8YPDHBLznsAR7
Behavior-preserving cleanups (guarded by the Phase 1 regression suite) plus the value-level absorption invariant test that would have caught 1.1-1.4. - terms_extract.jl (2.2): merge the parallel extract_operators/extract_affines and sort_and_extract_* families into one accessor-parameterized implementation; expand now preserves Term.repr; drop unused xt bindings. - terms_properties.jl (2.2): use affine(term) instead of term.A.L. - precomposeNonlinear.jl (2.5): drop the redundant ArrayPartition gradient! method (ArrayPartition <: AbstractArray). - addition.jl (2.5): collapse the duplicated +/- and broadcasted +/- bodies into _addsub / _broadcasted_addsub helpers. - test/test_phase2_absorption.jl (2.1): merge_function_with_operator satisfies absorbed_f(x) ~= λ·f(A·x + d) for eye/diagonal/AAᴴ-diagonal/general-linear/nonlinear. - AGENTS.md (2.5): document the retained-but-unused DifferentiationInterface (Phase 5) and AbstractFFTs deps. (The get_structure first-operator fix, another 2.5 item, shipped in the Phase 1 commit alongside the addition_tricky_part.jl length fix.) Deferred to follow-ups: scored algorithm matching (2.3), an explicit rejecting trait ruleset (2.4), and a repo-wide Runic pass. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01478URHkYh8YPDHBLznsAR7
… diagnostics Completes the deferred Phase 2 items (behavior-preserving, guarded by the Phase 1 regression tests and the full suite: 475 pass / 1 broken / 0 fail): - 2.2: `variables(::Variable)` now returns a 1-tuple `(x,)`, consistent with `variables(::Expression)`, removing the latent `Iterators.flatten` trap in `extract_variables`. Test assertion updated to the new contract. - 2.3: extract the greedy per-assumption matching loop in `parse_problem` into a named, documented, testable `match_assumption`/`candidate_term_subsets` pair. The selection score (largest-subset-first) is stated explicitly and reproduces the historical `reverse(collect(powerset(...)))` order exactly, so `parse_problem`/`suggest_algorithm`/`print_diagnostics` behavior is unchanged. - 2.4: `print_diagnostics(terms, algorithm)` now reports *why* each un-prepared term was rejected (the DCP-style unsatisfied property, e.g. `is_convex`), so a solver/problem mismatch — such as FastForwardBackward on a nonconvex problem — fails legibly at solve time instead of silently stalling. New `test/test_phase2_matching.jl` covers deterministic matching and the rejecting ruleset (parse rejection, solve error, diagnostic naming the failed property, and that ZeroFPR still accepts nonconvex smooth f). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01478URHkYh8YPDHBLznsAR7
Run Runic.jl 1.7.0 in-place over all of `src/`, ending the tab/space mix noted in the plan and normalizing spacing, argument wrapping, and explicit `return`s. Formatting-only (Runic is semantics-preserving); full suite unchanged at 475 pass / 1 broken / 0 fail. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01478URHkYh8YPDHBLznsAR7
The docs build was red once Phase 0 turned on `doctest=true`/the docs CI job. Make it build clean and extend user-facing coverage. Build fixes: - `modules = [StructuredOptimization]` (was `[…, ProximalAlgorithms]`, which forced documenting all 57 ProximalAlgorithms internals) + `checkdocs = :exports`. - Repair two latent docstring-attachment bugs: `solve`'s docstring was orphaned onto the Phase-1 `_run_solver` helper (moved the helper above it); `normalop_ls` had a blank line between its docstring and definition, so it never attached. - Drop the `@docs ZeroFPR/PANOC/PANOCplus` block (those docstrings live in ProximalAlgorithms, outside `modules`); describe the solvers in prose instead. - Add `[sources]`/`[deps]` to docs/Project.toml so the docs resolve the same locally-dev'd dependencies as the package (mirrors test/Project.toml). Content: - New `theory/parsing.md` (pipeline, trait propagation, operator-absorption cases, separable sums, matching/diagnostics) and `faq.md` (algorithm choice, "cannot parse" walkthrough, warm-starting, Float32, the fft prox trick). - Rewrote the Solvers page: auto-selection, algorithm-selection table, PANOC stepsize caveat, and a "when parsing fails" section — moving operational knowledge out of AGENTS.md into user docs. - Docstrings for `suggest_algorithm` and `print_diagnostics`. Reference hygiene (3.6): `problems`→`problem`, `maxiter`→`maxit`, `\nabla` escape, `SqrNormL2WithNormalOp` argument order, deploydocs repo URL, and stale kul-forbes links updated to JuliaFirstOrder/hakkelt. Full suite green (501 pass / 1 broken); docs build clean. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01478URHkYh8YPDHBLznsAR7
Add `test/test_phase4_coverage.jl` — value-asserting tests (computed values and captured diagnostics, not just line execution) across the worst-covered files. Raises total coverage from the 68% baseline to 86.7% (1014/1170); `parse.jl` 33% -> 73%, `term.jl` -> 99%, `sqrNormL2WithNormalOp` -> 91%. Coverage of the never-executed diagnostics/parsing paths surfaced three real bugs, fixed here: - `print_diagnostics(::TermSet, ::SimpleTerm, …)` called `findfirst` on a `TermSet`, which has no `keys`/`pairs` — it crashed instead of reporting the offending term (and cascaded through the OperatorTerm / InfimalConvolution multi-term paths that delegate to it). Search the collected vector instead. - `print_diagnostics(::Term, ::LeastSquaresTerm, …)` did `assumption.b.first`, but `assumption.b` is a bare `Symbol` — it threw a `FieldError`. Use `assumption.b`. - `a - b` for `a::Array`, `b::AbstractExpression` computed `-(b(x) + a)` = `-b(x) - a`, flipping the sign of the added constant; `c - A*x` evaluated to `-A*x - c` instead of `-A*x + c`. Negate the operator and add `a` once: `AffineAdd(-affine(b), a)`. Tests cover: displacement/variables (utils); TermSet show, scalar-mul, iteration, `==`, repr constructor, trait predicates (term); Moreau smoothing and conj (term/ prox); normal-op traits/value/`normalop_ls`; LeastSquares √λ scaling & b-sign, SquaredL2Term eye/diagonal/reject, OperatorTerm decomposition and is_eye branch; `print_diagnostics` for every assumption family (SimpleTerm/OperatorTerm/ InfimalConvolution/LeastSquares/SquaredL2/Repeated), the not-AAc-diagonal and incompatible-terms branches, single-element-TermSet delegations, and a sweep over every advertised algorithm; multi-variable separable solve; multi-variable `Usum_op` HCAT branches and array±expression. Note: `PrecomposedSlicedSeparableSum` (one variable split across several sliced terms) remains unimplemented (`prox!` undefined) — left out of scope; the sliced test uses the working multi-variable separable path instead. Full suite green; total 86.7%. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01478URHkYh8YPDHBLznsAR7
The generic multi-var/single-var Usum_op methods in addition.jl are reached whenever a nonlinear (or scaled) wrapper around a multi-variable expression is added to/subtracted from another expression -- e.g. sin(A*x + B*y) + C*z. Such a wrapper keeps several variables but is not an HCAT, so it bypasses the HCAT-specialized methods. These paths were reachable from ordinary syntax but previously unexercised by the suite. Add assertions for both operand orderings, both +/-, and the case where the single operand's variable is already present in the multi-var operand (the in-branch), verifying the resulting operator's action against a hand-computed value. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01478URHkYh8YPDHBLznsAR7
Manifest.toml is already listed in .gitignore; remove it from the index so the environment's resolved dependency versions are no longer version-controlled. The file is left in place on disk. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01478URHkYh8YPDHBLznsAR7
prox!/gradient! on the ProximalOperators functions built into term_kwargs
were called every solver iteration with no chance to preallocate their
scratch space, so any operator needing a buffer (weighted norms, LogisticLoss,
etc.) paid an allocation on every call. _run_solver now calls preallocate on
every value in term_kwargs once, using the actual x0 the solver will run
with, before the iteration starts -- matching ProximalOperators' new
preallocate(f, x) interface. Values with nothing to preallocate come back
unchanged.
PrecomposeNonlinear already builds its own eager scratch buffers (bufC etc.)
at construction time, but stored the inner g as-is; it now preallocates g for
bufC's shape too, so nested ProximalOperators functions used through a
nonlinear precomposition get the same treatment automatically at construction
time rather than needing the caller to know about it.
Verified via manual inspection that preallocate cascades correctly through
Precompose/PrecomposeDiagonal/PrecomposeNonlinear wrappers (is_preallocated
returns true on the wrapped LogisticLoss/SqrHingeLoss after a solve is set
up), and that a real least-squares + L1 solve still converges to the same
answer. The pre-existing extract_functions/SeparableSum test failures (1
fail + 3 errors in the full suite) are unrelated: they reproduce identically
against unmodified ProximalOperators and StructuredOptimization code --
extract_functions(t::TermSet) broadcasts to a Vector{Any}, and
SeparableSum(fs::Vararg) wraps a single Vector argument in a 1-tuple instead
of splatting it, a bug independent of this change.
Depends on ProximalOperators.jl's `preallocate`/`threaded` branch (not yet
merged to master).
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XAY78b7aJmkhatBJiW3iYd
…ableSum splat bug Merges the value/gradient-consistency and adjoint-scaling fixes for SqrNormL2WithNormalOp from the unmerged branches and the vendored MriReconstructionToolbox fork, generalized to keep the array-λ (weighted) gradient support already on this branch. ls now auto-detects the normal-op optimization for a single-variable expression with a non-identity linear operator, instead of requiring the separate normalop_ls function. Multi-variable expressions keep the plain path, since a multi-variable normal-op term's operator has to stay the identity on its own joint domain and so cannot be combined with unrelated-variable terms afterwards via the generic Term-extraction machinery (expand/extract_operators). Also fixes a real, pre-existing bug: extract_functions(t::TermSet) and extract_functions_nodisp(t::TermSet) passed SeparableSum a Vector instead of splatting it, so SeparableSum wrapped the whole vector as a single 1-tuple element instead of building one function per term. This was silently broken for any TermSet extraction with 2+ terms. Fixed test/Project.toml's [sources] to match the root Project.toml so the workspace Manifest.toml resolves the same dependency paths in both environments (was pointing ProximalOperators/AbstractOperators at stale paths lacking the preallocate feature). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Fixes the remaining CPU-array assumptions found while wiring up
end-to-end GPU testing through GPUEnv:
- ls's normal-op path and convert(Expression, ::Variable) both built an
Eye operator from (domain_type, size) alone, which defaults to a CPU
array type and silently dropped the real operator/array's GPU storage.
Now built from a real domain array/operator instead
(Eye(AbstractOperators.allocate_in_domain(...)) / Eye(~x)), so the
identity operator always matches the actual array type in play.
- PrecomposeNonlinear's scratch buffers were allocated with zeros(t, s)
(always a CPU Array); now use the existing
AbstractOperators.allocate_in_domain/allocate_in_codomain helpers,
matching the pattern already used in SqrNormL2WithNormalOp.
- hingeloss/sqrhingeloss/crossentropy were pinned to b::Array{R,1}, so a
GPU-array label vector couldn't dispatch at all; loosened to
AbstractVector{R} to match the already-generic logisticloss.
- _weighted_sqnorm's explicit indexing loop over lambda/d was scalar
GPU indexing; rewritten as a broadcast + reduction.
Adds test/test_gpu.jl (wired into runtests.jl), which runs ls/norm/
hingeloss and full problem()/solve() round trips through GPUEnv on
every GPUArrays-compatible backend found on the host (JLArrays always,
plus real CUDA where available), checking against the CPU result.
Requires the AbstractOperators.latest-stable / ProximalAlgorithms.jl /
ProximalOperators.jl checkouts this repo's [sources] point at to
already carry GPU support.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
normalop_ls was folded into ls() in a928555; the @docs block in docs/src/functions.md still referenced it, breaking Documenter's strict checkdocs build (no docstring exists for it anymore). ls()'s own docstring already covers the normal-op auto-detection behavior. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HTc1cqf5hWyg9CtYEfmiSP
`ls` can only fold a linear operator into `SqrNormL2WithNormalOp` for a single-variable expression. A multi-variable normal-op term would have to carry an identity over its own joint `ArrayPartition` domain, which collapses several variables into one operator domain and breaks the one-variable-per-domain invariant that `expand` and `_sort_and_extract` rely on when they pad and permute a term's operator up to the problem's full variable tuple. Retry the rewrite in `merge_function_with_operator` instead. By then the operator has been expanded to the problem's full domain and nothing is composed with it afterwards, so the joint domain is no longer a problem. It is applied on the general-linear branch, which already documents that prox is invalid there, so only the gradient is at stake and nothing is given up by folding. Two conditions gate the rewrite, so it never makes things worse: * `Lᴴ * L` has to fuse into a single operator rather than stay a `Compose(Lᴴ, L)`, which is what tells us a cheaper normal operator exists at all. An `HCAT` does not fuse on its own, so its block Gram `[LᵢᴴLⱼ]` is assembled explicitly as a `VCAT` of `HCAT` rows over the joint domain — and only when every one of the N² block products fuses, since the block form costs N² applications against the 2N of the operator and its adjoint. * `L` has to map into a codomain at least as large as its domain. `LᴴL` acts on the domain, so on a wide `L` it is both slower and worse conditioned. A least-squares term over several variables is the usual way to end up wide, its domain being the sum of the blocks' domains. The displacement is re-attached to the normal operator as `Aᴴd`, per block row for a block Gram: `AffineAdd` compares `size(d)`, a flat length for an `ArrayPartition`, against the operator's codomain size, which for a `VCAT` is a tuple of block sizes, so wrapping the whole operator is rejected. The product built for the applicability test is handed to the constructor through a new `pureAᴴA` keyword rather than computed a second time. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HTc1cqf5hWyg9CtYEfmiSP
`ls` folded a single-variable operator into a `SqrNormL2WithNormalOp` the
moment the term was built. That hid the operator from every later decision:
- the diagonal and AAᴴ-diagonal absorptions never saw it, so `ls(fft(x) - b)`
lost the exact prox the DFT case exists to provide, and a diagonal operator
was never folded into the weight;
- the normal operator was formed even where it is slower and worse
conditioned (a wide operator), which the parse-time path already declines
via `normal_op_worthwhile`;
- the resulting term advertised `is_proximable`, inherited from convexity,
although the function implements only `gradient!`.
`ls` now returns a plain `SqrNormL2` over the expression it was given, and
`merge_function_with_operator` — the one place that knows both the expanded
operator and what the selected algorithm asks of the term — picks the
formulation. The rule the syntax layer follows from here: build `λ·f(A·x + d)`
triples, nothing else.
Two defects surfaced once the operator reached the absorption:
- the diagonal branch folded a `SqrNormL2` into the weighted form while
silently dropping the displacement. `½‖diag(a)·x - b‖²` came out as
`½‖diag(a)·x‖²`. Reachable before this commit through `norm(D*x - b, 2)^2`.
- `SqrNormL2WithNormalOp` claimed proximability by default; it now declares
`is_proximable == false`.
Documented in docs/src/theory/parsing.md, which gains the normal-operator row
of the absorption table and the reason the rows are ordered as they are.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HTc1cqf5hWyg9CtYEfmiSP
The package stores a term as `λ · f(A·x + d)`. Two extractors implemented two different conventions for taking the function back out: `extract_functions` folded the displacement into a `PrecomposeDiagonal`, `extract_functions_nodisp` did not. Pairing the folding one with `extract_affines` (which keeps the displacement in the operator) counts `d` twice — the defect fixed in Phase 1 for the `prepare` methods, but left in place in three `print_diagnostics` methods, which therefore printed a decomposition that is not the one that would be solved. Closes PLAN.md 1.1/1.2 residue and 2.1: - `print_diagnostics` for `OperatorTerm`/`TermSet`, `OperatorTermWithInfimal- Convolution`/`Term` and `.../TermSet` now use the no-displacement convention, matching their `prepare` counterparts. - `extract_functions_nodisp` is renamed `weighted_function` (it applies λ and nothing else) and `extract_functions` is deleted. PLAN.md's instruction to delete `_nodisp` instead was written before the Phase 1 fixes made it the correct convention. - The one site that genuinely wants the displacement inside the function, `PrecomposedSlicedSeparableSum` (which is handed the linear blocks and precomposes them itself), gets an explicit local `fold_displacement` helper with a comment stating why it differs. - Two inline `λ == 1 ? f : Postcompose(f, λ)` copies in `prepare`/ `print_diagnostics` for `OperatorTerm` now call `weighted_function`. Suite: 778 pass / 1 broken (the expected Aqua ambiguities check), up from 773 by the five new assertions. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HTc1cqf5hWyg9CtYEfmiSP
Padding a term with `Zeros` blocks for the variables it does not mention was implemented twice: `expand` did it at the expression level by repeated `ex += Zeros(...) * x`, and `add_missing_vars` did it at the operator level for `Usum_op`. Same rule, two places to keep in step. `expand` now calls `add_missing_vars` and wraps the widened `(variables, operator)` pair back into an `Expression`, keeping the direct-HCAT construction PLAN.md 2.2 asks for. The `Term` method is a one-liner that carries λ, `f` and `repr` across unchanged. When nothing is missing, `add_missing_vars` returns its input, so `expand` is now the identity rather than a rebuild. New tests assert the padded block is a genuine zero block (the widened operator agrees with the original for every value of the added variable), that `repr`, λ and `f` survive, and that the no-op case returns the converted expression itself. Suite: 785 pass / 1 broken. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HTc1cqf5hWyg9CtYEfmiSP
`==(ex, b)` reached into the expression, demanded a `MatrixOp`, folded it into an `IndAffine` and returned a term over `variables(ex)[1]`. Everything else errored with "Currently affine equality supported only with `MatrixOp`" — including `DiagOp(a)*x == b`, a trivial projection, and `fft(x) == b`, which is AAᴴ-diagonal — and a multi-variable equality silently lost every variable after the first. The syntax layer now builds `Term(IndPoint(b), ex)` and leaves the choice of formulation to `merge_function_with_operator`, per PLAN.md 2.6. The existing diagonal and AAᴴ-diagonal branches cover the two cases that used to error, and a new `IndPoint` + `MatrixOp` rule reproduces today's `IndAffine` exactly. The proximability gate is restated to match. `is_AAc_diagonal(affine(term))` was a proxy for "absorbing the operator keeps an exact prox"; the `IndAffine` rule is a case the proxy does not cover. `keeps_exact_prox(op, f)` states the predicate directly, mirroring the branch table of `merge_function_with_operator`, and is called from `prepare(::Term, ::SimpleTerm, _)`, `prepare(::TermSet, ::SimpleTerm, _)` and `is_proximable(::Term)`. The two diagnostics messages that named the old proxy are reworded. Tests assert projections, not values: an indicator is `Inf` at almost every point, which made the previous `IndAffine` value comparisons vacuous. The diagonal and DFT cases are checked against their hand-computed feasible points (`b ./ a` and `x0`), the `MatrixOp` path against `prox(IndAffine(A, b), ...)` from both spellings of the constraint, and a multi-variable equality against the joint operator it is supposed to build. Suite: 810 pass / 1 broken. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HTc1cqf5hWyg9CtYEfmiSP
Two fixed orderings decided how a problem was parsed: the `if`-chain in `merge_function_with_operator`, and "first subset that prepares" in `match_assumption`. Both are now scored searches, and they are scored together — a cheaper formulation is only better if the algorithm that gets selected can use it (PLAN.md 2.3 and the first Phase 5 bullet). Formulation layer. `best_formulation(op, f, disp, λ, needs)` ranks the eight formulations this package can build by `(keeps an exact prox ? 0 : 1, cost)`, with the old branch order breaking exact ties, and `merge_function_with_operator` builds only the winner. Costs are in units of one application of `op` plus one of `opᴴ`, normalised so the generic `Precompose` costs 2; the full table is in the `best_formulation` docstring. `needs = :prox` restricts the search to formulations whose `prox!` is exact, which is how a caller states what the selected algorithm will ask of the term — the same question the `keeps_exact_prox` gate asks, so gate and filter cannot disagree. Algorithm layer. `match_assumption` scores every subset that prepares by `(-size, formulation cost, powerset position)` instead of returning the first, and prunes: subsets come largest-first, so the search stops at the end of the size class that first succeeded. `parse_problem(terms)` picks the algorithm whose complete parse is cheapest rather than the first that parses. The cost constraint was the design driver, and it exposed two real problems: - `fused_normal_op` answers "does `Lᴴ L` fuse?" by *building* the product — for a `MatrixOp` that is the Gram matrix, O(n²m), more than several solver iterations. `normal_op_fuses` now answers it from the types alone, via inference on `adjoint` and `*`, and is the only one scoring may call; `fused_normal_op` is reached once, for the winner. The type-level predicate is conservative (an inference result of `Any` counts as not fusing) and `merge_function_with_operator` falls back to `Precompose` if it is ever optimistic, so a wrong prediction costs a suboptimal formulation, never a wrong one. A test pins the two against each other on seven operator shapes. - `is_AAc_diagonal(::MatrixOp)` is `isdiag(A*Aᴴ)` upstream: O(m²n) and an m×m temporary. The parser has always paid this, once per proximability check. `is_aac_diagonal` first tries to *disprove* row orthogonality on a sample of row pairs, which settles anything not genuinely AAᴴ-diagonal in O(n), and falls through to the upstream check otherwise — so the answer is identical, not an approximation. It is also asked last and only when it can still change the winner. The full suite runs in 2m42s, down from 13m10s. Scoring now allocates nothing and its cost is independent of operator size: the tests assert `@allocated best_formulation == 0` for both a 10x8 and an 800x600 operator, and that scoring the whole problem costs under a fifth of a five-iteration PANOCplus pass (measured ~20x cheaper). Suite: 850 pass / 1 broken, with the pre-existing parse results unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HTc1cqf5hWyg9CtYEfmiSP
`solve` printed a full per-term diagnostic before failing, then threw "Sorry, I cannot parse this problem for solver of type ...". The report goes to stdout, so a caller that catches the error learns nothing from it — the opposite of the rejecting ruleset PLAN.md 2.4 asks for. `parse_failure_message` builds the message from the same `unsatisfied_reasons` the report uses, naming each unparseable term by its `repr` (what the user wrote, not the desugared operator graph) and the DCP-style property that blocked it. All three `solve` error paths use it. The solver-list path diagnoses against the solvers it was given rather than the whole registry: `ZeroFPR` parses the nonconvex example the test uses, so diagnosing globally would have produced a message claiming nothing was wrong. `closest_algorithm(terms, algorithms)` carries that choice and also replaces the inline loop in `print_diagnostics(terms)`. Suite: 856 pass / 1 broken. (An earlier run of this same tree had `Aqua.test_persistent_tasks` fail; it spawns a subprocess under a timeout and was flaky under node load, and passes on re-run with no change to the tree.) Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HTc1cqf5hWyg9CtYEfmiSP
…eshold `benchmark/benchmarks.jl` exports `SUITE` in AirspeedVelocity's convention, with `benchmark/Project.toml` mirroring the package's `[sources]` so it measures the same code the tests run against. Note the singular directory: `benchmarks/` (plural) holds the documentation demo scripts and is untouched. Four groups, each guarding a claim the code or the documentation makes: `formulation/` (the normal operator against `Precompose` across tall, square, mildly wide and wide operators, plus a non-fusing operator), `block_gram/` (the assembled block Gram against the two-pass `HCAT`), `absorption/` (the diagonal and AAᴴ-diagonal prox tricks against the naive forms) and `parse/` (the scoring budget against a five-iteration solve). `.github/workflows/benchmark.yml` runs it PR-vs-base through `benchpkg`, informational only (`continue-on-error`), with the table in the job summary. It installs the *registered* AirspeedVelocity rather than the local fork at /project/c_mrrecon/AirspeedVelocity.jl: that fork is six commits ahead of upstream with table-formatting and emoji changes only, and CI must not depend on a checkout that exists on one machine. The measurements replace "set from one observed regression" in the `normal_op_worthwhile` docstring with a table. They confirm the `n <= m` threshold rather than moving it: the per-iteration saving collapses to a few percent as soon as `n > m` while the one-off Gram construction keeps growing, so break-even moves from ~44-50 iterations (tall, square) to ~400-950 (mildly wide, wide) — before counting the squared condition number, which the timings do not capture at all. AGENTS.md gains a Benchmarks section, including the singular/plural distinction. Suite: 856 pass / 1 broken (this commit changes one docstring in `src/`). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HTc1cqf5hWyg9CtYEfmiSP
`DifferentiationInterface` and `AbstractFFTs` were declared in `[deps]` and `[compat]` but referenced nowhere in `src/`. AGENTS.md justified keeping them: `DifferentiationInterface` was reserved for the Phase 5 differentiable-solvers / unrolling work, which has been dropped by decision, and `AbstractFFTs` was said to be what the FFT bindings resolve against — they in fact go through `FFTW`/`FFTWOperators`, which pull it in transitively anyway. Both are removed and AGENTS.md updated to say so. Formatting can no longer regress silently: `.github/workflows/format.yml` runs Runic in `--check --diff` mode over `src/`, `test/` and `benchmark/`. `test/` was never Runic-formatted (tabs throughout `runtests.jl`, mixed spacing elsewhere), so this commit formats it, which is most of the diff. PLAN.md asked for a `.JuliaFormatter.toml` "pinning the Runic style". There is no such thing — JuliaFormatter has no Runic style — and a config file would only point editors at a second formatter that disagrees with the one CI enforces. The check job is the part that has teeth, and AGENTS.md now records why the config file is deliberately absent. Also gitignores `.serena/`, the Serena symbol-index cache. Suite: 856 pass / 1 broken. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HTc1cqf5hWyg9CtYEfmiSP
Closes PLAN.md Phase 3 and the part of Phase 0.3 that was configured but never enabled by any content. New theory pages. `theory/problem_form.md` defines the vocabulary the parser uses — convex, strongly convex, smooth, proximable, generalized quadratic, set indicator — with three proximal operators worked out in closed form, and renders the algorithm-assumption table from `ProximalAlgorithms.get_assumptions` in an `@example` block so it cannot drift from the code. `theory/matrix_free.md` explains why no matrix is stuffed, tabulates what each expression composes into, and states the normal-operator trick together with both of its caveats: the adjoint-scaling correction that keeps the returned value consistent with the gradient, and the squared condition number. New `docs/examples/`, seven Literate.jl pages executed during the build, so an example that stops working fails CI: lasso with a warm-started regularization path, TV denoising, audio declipping, a multi-variable problem, Rosenbrock, FFT deconvolution with the formulation choice measured, and a "when parsing fails" walkthrough. Writing them found two things worth recording in the pages themselves: `norm(A*x, 1)` is refused by proximal-gradient methods but accepted by algorithms with an operator slot (the earlier claim that nothing parses it was simply wrong), and a scalar index like `x[1]` produces an expression the gradient path cannot handle, so Rosenbrock uses `x[1:1]`. Doctests now exist. Nine `jldoctest` blocks in `build_solve.jl`, `minimize.jl` and `proximalOperators_bind.jl` replace `julia` blocks that were never executed — one of which claimed `Variable(4)` prints as `Variable(Float64, (4,))`, which it has not done for some time. `DocMeta.setdocmeta!` supplies the imports and a fixed seed, and `doctestfilters` absorb float noise and type parameters. The assertions are structural (types, counts, satisfied constraints) rather than printed floats. `checkdocs = :all` makes an undocumented-and-unreferenced docstring fail the build. That required docstrings for `TermSet`, `rank`, `^` on a `NormL2` term and `PrecomposeNonlinear`, and a new `internals.md` page collecting the parser and normal-operator internals — which is a page worth having anyway. Also: `docs/Project.toml` `[sources]` now match the package's exactly (they pointed at different checkouts of AbstractOperators and ProximalOperators, so the documented API was not the tested one), and README's badges and links point at `hakkelt` rather than the stale `JuliaFirstOrder`. Docs build clean with `checkdocs = :all` and doctests executing. Suite: 856 pass / 1 broken. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HTc1cqf5hWyg9CtYEfmiSP
Closes PLAN.md Phase 4, measured against the finished code rather than the stale Jul 3 figure the plan quoted. Measured with LocalCoverage before and after. `parse.jl`, the file the plan singles out, goes from 74.78% to 93.45%; the total from 87.91% (1127/1282) to 95.19% (1208/1269), which meets the 90% target and the ~95% stretch. The new tests, all value- or message-asserting: - Sliced separable sums, both the disjoint case the parser accepts and the overlapping one it refuses, with the diagnostic checked. - `RepeatedSimpleTerm` and `RepeatedOperatorTerm`, single-term and TermSet, accepting and rejecting, with the collected functions checked against the terms they came from. - The infimal-convolution paths: the `func₁` and `func₂` branches, the multi-variable fallback that attaches a block identity, and the diagnostics. - The rejection paths of `LeastSquaresTerm` and `SquaredL2Term`, including a term already folded into the normal-operator formulation, whose target is read back out of the operator's displacement. - `PrecomposeNonlinear`'s value, and the adjoint-scaling probe's fallback branch (an operator that annihilates the constant probe vector). Two things the coverage measurement exposed, both recorded rather than papered over: - The multiple-terms-per-variable branch of `prepare_proximable_single_var_per_term` is **unreachable**: its only caller enters it exactly when every variable has one term. Deleted, with a comment saying where the sliced case is actually handled. - `PrecomposedSlicedSeparableSum` does not handle the shape this package builds for one variable with several sliced terms: its value disagrees with the sum of its own pieces, and `prox!` throws while iterating the per-variable operator. The defect is in the pinned ProximalOperators, in a path nothing exercised before this commit. Two `@test_broken` assertions pin it so a fix flips them green; the surrounding assertions check that the parser hands over the right pieces, which it does. The two per-file targets not met are `sqrNormL2WithNormalOp.jl` at 94.23% (against ≥95%) and `utils.jl` at 8/9 lines: what is left in both is single-line trait definitions that are const-folded at the call site and never counted. Also: `codecov.yml`'s target moves from the 68% baseline to the measured 95%, and the AGENTS.md coverage recipe is corrected — it said `--project=test`, but LocalCoverage runs the suite in a subprocess and must *not* be in the test environment, which is why it is absent from `test/Project.toml`. The recipe now sets up an environment of its own and shows how to read `coverage_gaps`, which is what tells you which lines to write a test for. Suite: 917 pass / 3 broken (Aqua ambiguities plus the two above) / 0 fail. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HTc1cqf5hWyg9CtYEfmiSP
normalop_ls previously crashed on multi-variable expressions
(MethodError: no method matching length(::AffineAdd{HCAT{...}})).
- expression.jl: allow ndoms(L,1) > 1 when is_eye(L) (DCAT(Eye,Eye) is
a provable no-op, needed for the multi-variable eye_op below).
- proximalOperators_bind.jl: build eye_op as
Eye(ArrayPartition((~xi for xi in ex.x)...)) instead of iterating
the (non-iterable, AffineAdd-wrapped) operator ex.L. This preserves
each variable's own domain block (a block-identity DCAT), not a
collapsed shared codomain.
- term.jl: add _scalar_codomain_type helper since codomain_type of a
multi-domain operator returns a nested Tuple, not a scalar Type;
Term's inner constructor needs a scalar for its real()/eltype checks.
Depends on AbstractOperators HCAT.has_optimized_normalop/get_normal_op
(fork PR hakkelt/AbstractOperators.jl#6) for the fast Toeplitz-NFFT
normal-operator path to actually fire on multi-variable problems;
without it the fix is still correctness-only (no crash) but the fast
path silently falls back to a plain Compose.
Verified: multi-variable normalop_ls builds, gradient matches ls
formulation, FastForwardBackward reaches same minimizer as ls; HCAT
normal-op fusion confirmed exact + ~3x faster against a real NFFTOp;
single-variable path untouched (bit-identical). Full test suite:
430/432 (2 non-passing are pre-existing unrelated flakiness/known-
broken, confirmed via isolated re-run, not caused by this change).
`ls(A*x - y)` carries the displacement `-y`, but the solvers that consume the `LeastSquaresTerm` assumption (ADMM, CG, CGNR) minimize `‖Lx - b‖²`, so the displacement has to reach them negated. Passing it through unchanged made those solvers minimize `‖Lx + y‖²`, returning `-x`: the correct magnitude with a flipped sign. The `SimpleTerm` path is unaffected, because `merge_function_with_operator` feeds the displacement to `Precompose`/`PrecomposeDiagonal`, which use the `f(x + disp)` convention. Add signed correctness assertions for both affected paths in test_usage_small.jl: CGNR against the closed-form ridge solution, and ADMM against a PANOCplus reference on an overdetermined l1 problem. The existing assertions there only checked `!isnothing(sol)`, which a sign flip cannot fail, and the one numerical assertion covers ZeroFPR — a solver on the unaffected path. ADMM needs an explicit `rho` here because its default adaptive penalty sequence stalls on this problem. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Rsct6ZoJyGhXz89g7eEQmC
`SqrNormL2WithNormalOp` builds `AᴴA = A' * A` eagerly in its constructor, but `prepare(::LeastSquaresTerm, ...)` passed the solver only `f.A`, so an algorithm that needs the normal operator built a second, identical `Compose` chain of its own — two independent chains, each with its own buffers. `LeastSquaresTerm` gains an optional `AHA` assumption symbol, so only algorithms that actually form the normal operator ask for one (ADMM does; the CG family does not and is unchanged). `prepare` fills it with `remove_displacement(f.AᴴA)` — the linear part, matching the `remove_displacement` already applied to `op`, since ADMM carries `b` separately — and only when `lambda == 1`, because a different `lambda` rescales `op` and the cached `AᴴA` would no longer match. Pairs with the ProximalAlgorithms side (`ADMMIteration`'s `AHA` field). On MRT's benchmark TV problem the saving is ~4 MiB and ~7.6 ms on a one-outer-iteration solve, and inside noise over 30 iterations: this is a setup-cost change, not a steady-state one. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…tor skip the size test `fused_normal_op` gated everything on `normal_op_worthwhile`, a dense-matrix cost estimate that requires the codomain to be at least as large as the domain. That estimate cannot model a structured operator, and it vetoes exactly the case the normal-operator form exists for: a subsampled Fourier encoding maps into a codomain *smaller* than its domain, and its normal operator is still much the cheaper of the two. When `AbstractOperators.has_optimized_normalop(L)` is true the operator has already answered the question, so `L' * L` is returned without consulting the estimate. That product need not collapse to a single operator -- `get_normal_op(::Compose)` fuses the innermost adjoint pair and keeps the outer factors, so an encoding operator `S`-then-`F` becomes `Sᴴ·(FᴴF)·S`, one transform where the naive form needs two. Without this, an MRI reconstruction silently loses the fused normal operator and pays two transforms per gradient. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…advertise one `normal_op_applicable` is the predicate `best_formulation` scores candidates with, and it did not carry the shortcut `fused_normal_op` has, so an operator advertising an optimized normal operator was never considered for the fold at all. Both of its remaining tests reject exactly that operator: `normal_op_worthwhile` because such an operator typically maps into a smaller codomain than its domain, and `normal_op_fuses` because `get_normal_op(::Compose)` fuses only the innermost adjoint pair and stays a `Compose`. For an MRI encoding operator (32x32 image, 3x undersampled, 4 coils) the parser was choosing `:precompose` at cost 2.0 over a normal-operator form costing n/m = 0.73, and the reconstruction came out measurably worse: relative error 0.33 against 0.24 for the same seed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.