From d86e9942eba9ce6fc8d21e0c9b88f027d3388e44 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 17 Jul 2026 08:16:09 -0500 Subject: [PATCH] Declare the extension hooks public and document their contracts foreach_support and foreach_support_sym are how a cover algorithm reads a matrix: specializing them is what lets a storage type be covered in time proportional to its support, as the bundled SparseArrays extension does. AbstractCoverPenalty is the supertype of the exported penalties and the first-argument slot of most of the API. All three were reachable and relied upon, but none declared its status, so users had no way to tell a supported extension point from an internal that could move. Each now states what an implementer must guarantee. Two points that the code enforced only by convention are written down: - foreach_support_sym requires A to be symmetric in value, not merely square. It cannot check this cheaply and does not try; the precondition is what makes reporting one member of each index pair sufficient. - Subtyping AbstractCoverPenalty buys cover_objective and nothing more. The solvers dispatch on the concrete built-in penalties, so a custom penalty can score a cover but cannot be minimized. Saying so beats letting a MethodError be the documentation. `public` is parsed as a keyword only from Julia 1.11, so the declaration goes through a VERSION gate rather than a new Compat dependency. Assisted-by: Claude Opus 4.8 --- src/ScaleInvariantAnalysis.jl | 5 +++ src/penalties.jl | 30 ++++++++++++++ src/support.jl | 74 ++++++++++++++++++++++++++++++----- 3 files changed, 100 insertions(+), 9 deletions(-) diff --git a/src/ScaleInvariantAnalysis.jl b/src/ScaleInvariantAnalysis.jl index b1db7bc..0ab7c16 100644 --- a/src/ScaleInvariantAnalysis.jl +++ b/src/ScaleInvariantAnalysis.jl @@ -12,6 +12,11 @@ export initialize_cover, initialize_cover!, initialize_symcover, initialize_symc export symcover_min, symcover_min!, cover_min, cover_min! export soft_symcover_min, soft_symcover_min!, soft_cover_min, soft_cover_min! +# `public` is parsed as a keyword only from Julia 1.11; this package supports 1.10. +@static if VERSION >= v"1.11" + eval(Meta.parse("public AbstractCoverPenalty, foreach_support, foreach_support_sym")) +end + include("penalties.jl") include("support.jl") include("heuristic_covers.jl") diff --git a/src/penalties.jl b/src/penalties.jl index 785103d..8286f3a 100644 --- a/src/penalties.jl +++ b/src/penalties.jl @@ -2,6 +2,36 @@ # φ types # ============================================================ +""" + AbstractCoverPenalty <: Function + +Supertype of the penalty functions `ϕ` that score a cover, and the type of the +first argument of most of this package's API. The built-in subtypes are +[`AbsLog`](@ref) and [`AbsLinear`](@ref). + +A penalty is a function of the single ratio `r = |A[i,j]| / (a[i]*b[j])`, and +[`cover_objective`](@ref) sums it over the entries of `A`. Because `ϕ` sees only +that ratio, and every diagonal rescaling of `A` leaves it fixed, any objective +built from a penalty is automatically scale-invariant. + +# Extending + +A subtype must be callable on a nonnegative real: + + (::MyPenalty)(r::Real) + +`r` ranges over `[0, Inf]`. Both endpoints occur and neither may error: `r = 0` +whenever `A[i,j]` is zero, and `cover_objective` passes `typemax` for an entry +left uncovered by a zero scale. Penalties are conventionally singleton structs. + +That call is the whole contract, and it buys exactly one thing: +[`cover_objective`](@ref) works for any subtype. **The solvers do not.** Every +solver in this package dispatches on a concrete built-in penalty — `AbsLog{2}` +is solved natively, the `AbsLinear` penalties through JuMP — so a custom subtype +passed to [`symcover_min`](@ref), [`soft_symcover`](@ref), or any other solver +raises a `MethodError`. Scoring covers with your own penalty is supported; +minimizing it is not. +""" abstract type AbstractCoverPenalty<:Function end """ diff --git a/src/support.jl b/src/support.jl index ef8c86b..b5e9cdf 100644 --- a/src/support.jl +++ b/src/support.jl @@ -2,19 +2,38 @@ # (geometric-mean init, feasibility boost, tightening) so each is written once # instead of once per storage type. # -# `foreach_support(f, A)` calls `f(i, j, v)` for every stored entry with -# `v = abs(A[i, j]) != 0`, in a storage-friendly order. -# -# `foreach_support_sym(f, A)` calls `f(i, j, v)` once per unordered index pair -# in a canonical triangle (including the diagonal), for use when `A` is known -# to be symmetric-valued. For `Tridiagonal`/`Bidiagonal`, whose two off-diagonal -# bands need not agree, the pair's value is `max(|A[i,i+1]|, |A[i+1,i]|)` -# (matching what a full-grid tighten already enforces on both entries). -# # Both are higher-order functions rather than iterators so that `f` is # specialized and inlined into a tight loop at each call site; every index # used is the matrix's own (`axes`, `eachindex`), so offset axes are honored. +""" + foreach_support(f, A) + +Call `f(i, j, v)` once for every entry of `A` whose magnitude +`v = abs(A[i, j])` is nonzero, and return `nothing`. Entries that are zero are +skipped, so `f` never sees `v == 0`. The order is whatever suits `A`'s storage +and is not part of the contract; `i` and `j` are `A`'s own indices, so offset +axes are honored. + +This is the hook through which cover algorithms read a matrix. Specializing it +is what lets a storage type be covered in time proportional to its support +rather than to `length(A)` — the package's own `SparseMatrixCSC` methods, which +walk `nzrange` instead of the full grid, are the model. + +# Extending + +To support a new matrix type, define + + ScaleInvariantAnalysis.foreach_support(f, A::MyMatrix) + +which must call `f(i, j, abs(A[i, j]))` exactly once for each `(i, j)` with +`abs(A[i, j]) != 0`, must not call `f` for any other entry (a stored zero is +still a zero), and must return `nothing`. Emitting an entry twice double-counts +it in the objective; omitting one silently drops a constraint, yielding a +"cover" that does not cover. + +See also: [`foreach_support_sym`](@ref). +""" function foreach_support(f, A::AbstractMatrix) for j in axes(A, 2) for i in axes(A, 1) @@ -25,6 +44,43 @@ function foreach_support(f, A::AbstractMatrix) return nothing end +""" + foreach_support_sym(f, A) + +Symmetric counterpart of [`foreach_support`](@ref): call `f(i, j, v)` once per +unordered index pair rather than once per entry, and return `nothing`. Pairs are +reported in the canonical orientation `i <= j`, the diagonal included, with +`v = abs(A[i, j])`; zero pairs are skipped. `A` must be square, or a +`DimensionMismatch` is thrown. + +`A` must also be **symmetric in value**, not merely square — this is a +precondition the function cannot check cheaply and does not try to. It is what +makes reporting one member of each pair sufficient: a symmetric cover +constrains `a[i]*a[j]` by a single magnitude, so visiting `(j, i)` as well would +only duplicate it. Handing this an asymmetric matrix does not error; it silently +covers a symmetrization of it, and which one is not specified. + +(The `Bidiagonal`/`Tridiagonal` methods, whose two off-diagonal bands are stored +separately and need not agree, report `max(|A[i,j]|, |A[j,i]|)` — the value a +full-grid tighten would enforce on both entries. Under the precondition the two +bands agree and this is just `abs(A[i,j])`; outside it, the choice is +robustness, not a promise.) + +# Extending + +To support a new matrix type, define + + ScaleInvariantAnalysis.foreach_support_sym(f, A::MyMatrix) + +which must call `f(i, j, v)` exactly once for each pair `i <= j` with +`v = abs(A[i, j]) != 0`, must not call `f` for zero pairs, and must return +`nothing`. Reporting the same pair in both orientations double-counts it. A +type whose storage is triangular (`Symmetric{<:Any,<:SparseMatrixCSC}` in this +package's own extension) must map stored `(i, j)` with `i > j` back to `(j, i)` +rather than emit it as found. + +See also: [`foreach_support`](@ref). +""" function foreach_support_sym(f, A::AbstractMatrix) ax = axes(A, 1) axes(A, 2) == ax || throw(DimensionMismatch("foreach_support_sym requires a square matrix, got axes $(axes(A))"))