From 1abfb48f1d6719b39b87db3e50275215b2143ba3 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 17 Jul 2026 08:30:06 -0500 Subject: [PATCH] Add iscover, the cover predicate a[i]*b[j] >= abs(A[i,j]) is the inequality this package is organized around, but there was no function to ask it. The solvers satisfy it by construction, so the gap went unnoticed; soft covers do not, which makes "does this actually cover A?" a real question a caller cannot currently put to the package. The test suite had already answered it by carrying its own predicate, which is the usual sign one belongs in the package. Coverage is a walk over the support: an entry that is zero constrains nothing, so only nonnegative scales make it sound to skip it. That precondition is enforced rather than documented -- a negative scale raises rather than silently reporting a cover that isn't one -- and it lets the check run in time proportional to the support, through the same foreach_support hook the solvers use. rtol and atol supply the slack a producing algorithm warrants, and default to none, so the bare call is the exact inequality. atol is subtracted only when nonzero, since abs(A[i,j]) - atol is undefined when the two carry different units. The tests now exercise the shipped function instead of a private copy. Assisted-by: Claude Opus 4.8 --- src/ScaleInvariantAnalysis.jl | 3 +- src/iscover.jl | 85 ++++++++++++++++++++++++++++++++ test/helpers.jl | 24 --------- test/iscover.jl | 92 +++++++++++++++++++++++++++++++++++ test/runtests.jl | 3 +- 5 files changed, 181 insertions(+), 26 deletions(-) create mode 100644 src/iscover.jl create mode 100644 test/iscover.jl diff --git a/src/ScaleInvariantAnalysis.jl b/src/ScaleInvariantAnalysis.jl index 947b0a8..c7396c7 100644 --- a/src/ScaleInvariantAnalysis.jl +++ b/src/ScaleInvariantAnalysis.jl @@ -6,7 +6,7 @@ using PrecompileTools: PrecompileTools, @compile_workload using Random: Random, AbstractRNG, MersenneTwister export AbsLog, AbsLinear -export cover_objective +export cover_objective, iscover export cover, cover!, symcover, symcover! export soft_symcover, soft_symcover!, soft_cover, soft_cover! export initialize_cover, initialize_cover!, initialize_symcover, initialize_symcover! @@ -20,6 +20,7 @@ end include("penalties.jl") include("support.jl") +include("iscover.jl") include("heuristic_covers.jl") include("initializers.jl") # the start menu; consumed by both solver families below include("soft_covers.jl") diff --git a/src/iscover.jl b/src/iscover.jl new file mode 100644 index 0000000..2a68554 --- /dev/null +++ b/src/iscover.jl @@ -0,0 +1,85 @@ +# The cover predicate. Kept beside the traversal it is built on: checking coverage is +# exactly a walk over the support, since an entry that is zero constrains nothing. + +""" + iscover(a, b, A; rtol=0, atol=0) + iscover(a, A; rtol=0, atol=0) + +Test whether `a` and `b` cover `A`, that is, whether `a[i]*b[j] >= abs(A[i,j])` for every +entry. The two-argument form tests the symmetric cover `a*a'`, and requires `A` to be +square. + +This is the inequality the whole package is organized around. [`cover`](@ref), +[`symcover`](@ref), and the `*_min` solvers all satisfy it by construction — up to the +roundoff noted below — so the predicate earns its keep mainly on covers that carry no such +guarantee: those from [`soft_cover`](@ref) and [`soft_symcover`](@ref), which penalize +under-coverage rather than forbidding it, and covers a caller has adjusted by hand. + +`rtol` and `atol` supply the slack the producing algorithm warrants, testing +`a[i]*b[j] >= abs(A[i,j])*(1 - rtol) - atol`. Neither defaults to any slack at all, so the +bare call is the exact inequality. Use `rtol` for roundoff that scales with entry magnitude +(the log-domain arithmetic behind the heuristics warrants a few multiples of `eps`), and +`atol` for the convergence tolerance of an iterative solver. `atol` is subtracted only when +it is nonzero: `abs(A[i,j]) - atol` is undefined when the two carry different units, so an +`rtol`-only check never forms it. That makes `rtol` the only one of the two meaningful for a +dimensional `A`, whose entries need not share units — no single scalar `atol` is +commensurate with all of them. + +`a` and `b` must be nonnegative; a negative scale raises an `ArgumentError`. Zero is +allowed, and is what every solver here returns for a row or column with no support. The +requirement is not cosmetic: only nonnegativity makes it sound to skip the zero entries of +`A`, which is what lets this run in time proportional to the support rather than to +`length(A)`. + +`eachindex(a)` must match `axes(A, 1)` and `eachindex(b)` must match `axes(A, 2)`. + +See also: [`cover_objective`](@ref), [`cover`](@ref), [`symcover`](@ref). + +# Examples + +```jldoctest +julia> A = [1.0 2.0; 3.0 4.0]; + +julia> a, b = cover(A); + +julia> iscover(a, b, A; rtol=8eps()) +true + +julia> iscover([1.0, 1.0], [1.0, 1.0], A) # a*b' = ones, which does not reach A[2,2] +false +``` +""" +function iscover(a::AbstractVector, b::AbstractVector, A::AbstractMatrix; rtol=0, atol=0) + eachindex(a) == axes(A, 1) || + throw(DimensionMismatch("indices of `a` must match row-indexing of `A`, got eachindex(a)=$(eachindex(a)), axes(A, 1)=$(axes(A, 1))")) + eachindex(b) == axes(A, 2) || + throw(DimensionMismatch("indices of `b` must match column-indexing of `A`, got eachindex(b)=$(eachindex(b)), axes(A, 2)=$(axes(A, 2))")) + _require_nonneg(a, "a") + _require_nonneg(b, "b") + # Zero entries of `A` are skipped by `foreach_support`, and need no check: they demand + # `a[i]*b[j] >= 0`, which nonnegative scales satisfy outright. + covered = Ref(true) + foreach_support(A) do i, j, v + covered[] &= _iscovered(a[i] * b[j], v, rtol, atol) + end + return covered[] +end + +function iscover(a::AbstractVector, A::AbstractMatrix; kwargs...) + axes(A, 1) == axes(A, 2) || + throw(DimensionMismatch("iscover(a, A) requires a square matrix, got axes $(axes(A))")) + return iscover(a, a, A; kwargs...) +end + +_iscovered(p, v, rtol, atol) = iszero(atol) ? p >= v * (1 - rtol) : p >= v * (1 - rtol) - atol + +function _require_nonneg(x::AbstractVector, name::String) + for i in eachindex(x) + # `zero(x[i])`, not `zero(eltype(x))`: a dimensional scale carries its units in the + # value, and `zero(Quantity{Float64})` is undefined. Also rejects NaN, which fails + # every comparison. + x[i] >= zero(x[i]) || + throw(ArgumentError("iscover requires nonnegative scales, got $name[$i] = $(x[i])")) + end + return nothing +end diff --git a/test/helpers.jl b/test/helpers.jl index 27dd8a5..436e01b 100644 --- a/test/helpers.jl +++ b/test/helpers.jl @@ -3,30 +3,6 @@ # The four penalty instances, for loops over properties that must hold for every ϕ. const PENALTIES = (AbsLog{1}(), AbsLog{2}(), AbsLinear{1}(), AbsLinear{2}()) -""" - iscover(a, b, A; rtol=0, atol=0) - iscover(a, A; rtol=0, atol=0) - -Cover feasibility: `a[i]*b[j] >= abs(A[i,j])` for every entry, with slack for -the accuracy the producing algorithm warrants. `rtol` absorbs roundoff that -scales with entry magnitude (the log-domain arithmetic used by the heuristics); -`atol` absorbs the convergence tolerance of iterative solvers. Heuristic covers -are feasible by construction and should pass with `rtol` a few multiples of -`eps` and `atol=0`; the `*_min` solvers need the `atol` they converge to. - -Only `rtol` applies to a dimensional `A`, whose entries need not share units: -no one scalar `atol` is commensurate with every entry. -""" -function iscover(a, b, A; rtol=0, atol=0) - return all(_iscovered(a[i] * b[j], abs(A[i, j]), rtol, atol) - for i in axes(A, 1), j in axes(A, 2)) -end - -# `atol` is subtracted only when it is nonzero, so that an `rtol`-only check never -# forms `abs(A[i,j]) - atol` -- undefined when the two carry different units. -_iscovered(p, v, rtol, atol) = iszero(atol) ? p >= v * (1 - rtol) : p >= v * (1 - rtol) - atol -iscover(a, A; kwargs...) = iscover(a, a, A; kwargs...) - """ isbalanced(a, b, A; atol=1e-8) diff --git a/test/iscover.jl b/test/iscover.jl new file mode 100644 index 0000000..a35e447 --- /dev/null +++ b/test/iscover.jl @@ -0,0 +1,92 @@ +@testset "iscover" begin + A = [1.0 2.0; 3.0 4.0] + + @testset "the exact inequality, with no default slack" begin + @test iscover([1.0, 3.0], [1.0, 2.0], A) # a*b' == A exactly + @test !iscover([1.0, 1.0], [1.0, 1.0], A) # ones cannot reach A[2,2] + # Exactly tight is covered; a hair under is not. + @test iscover([1.0], [1.0], fill(1.0, 1, 1)) + @test !iscover([1.0], [prevfloat(1.0)], fill(1.0, 1, 1)) + @test iscover([1.0], [prevfloat(1.0)], fill(1.0, 1, 1); rtol=1e-15) + end + + @testset "slack" begin + Aone = fill(1.0, 1, 1) + @test !iscover([0.9], [1.0], Aone) + @test iscover([0.9], [1.0], Aone; rtol=0.2) + @test iscover([0.9], [1.0], Aone; atol=0.2) + @test !iscover([0.9], [1.0], Aone; rtol=0.05) + @test !iscover([0.9], [1.0], Aone; atol=0.05) + end + + @testset "symmetric form" begin + Asym = [4.0 1.0; 1.0 3.0] + @test iscover([2.0, 2.0], Asym) + @test !iscover([1.0, 1.0], Asym) + @test iscover(symcover(Asym), Asym; rtol=8eps()) + # The two-argument form is the symmetric cover a*a'. + a = [2.0, 2.0] + @test iscover(a, Asym) == iscover(a, a, Asym) + @test_throws DimensionMismatch iscover([1.0, 1.0], [1.0 2.0 3.0; 4.0 5.0 6.0]) + end + + @testset "zero entries constrain nothing" begin + # A zero entry is satisfied by any nonnegative scales, including zero ones. + Az = [1.0 0.0; 0.0 0.0] + @test iscover([1.0, 0.0], [1.0, 0.0], Az) + @test iscover([1.0, 0.0], Az) + # ... but a nonzero entry backed by a zero scale is not covered. + @test !iscover([1.0, 0.0], [1.0, 0.0], [1.0 0.0; 0.0 1.0]) + end + + @testset "negative scales are rejected" begin + @test_throws "iscover requires nonnegative scales" iscover([-1.0, 3.0], [1.0, 2.0], A) + @test_throws "iscover requires nonnegative scales" iscover([1.0, 3.0], [-1.0, 2.0], A) + @test_throws "iscover requires nonnegative scales" iscover([-1.0, 1.0], [4.0 1.0; 1.0 3.0]) + # NaN fails every comparison, so it is rejected too rather than silently passing. + @test_throws "iscover requires nonnegative scales" iscover([NaN, 3.0], [1.0, 2.0], A) + # The check runs before the traversal, so it fires even when no entry is violated. + @test_throws "iscover requires nonnegative scales" iscover([-1.0, -1.0], [-1.0, -1.0], zeros(2, 2)) + end + + @testset "axes must match" begin + @test_throws DimensionMismatch iscover([1.0], [1.0, 2.0], A) + @test_throws DimensionMismatch iscover([1.0, 3.0], [1.0], A) + Ao = OffsetArray(A, -1:0, -1:0) + @test_throws DimensionMismatch iscover([1.0, 3.0], [1.0, 2.0], Ao) + @test iscover(OffsetArray([1.0, 3.0], -1:0), OffsetArray([1.0, 2.0], -1:0), Ao) + end + + @testset "storage types agree with the dense reference: $(typeof(S))" for S in ( + sparse([1.0 2.0; 3.0 4.0]), + Diagonal([2.0, 3.0]), + SymTridiagonal([2.0, 3.0], [1.0]), + Tridiagonal([1.0], [2.0, 3.0], [0.5]), + Bidiagonal([2.0, 3.0], [1.0], :U), + Symmetric(sparse([4.0 1.0; 1.0 3.0]))) + D = Matrix(S) + for a in ([1.0, 1.0], [2.0, 2.0], [3.0, 3.0]) + @test iscover(a, a, S) == iscover(a, a, D) + end + end + + @testset "dimensional scales" begin + # A scale carries its units in the value, not the type: `zero(Quantity{Float64})` + # is undefined, so the nonnegativity check must ask the element, not the eltype. + L = u"m" + Au = [1.0/L^2 2.0/L^2; 2.0/L^2 4.0/L^2] + au = [2.0/L, 2.0/L] + @test iscover(au, Au) + @test !iscover([1.0/L, 1.0/L], Au) + @test iscover(symcover(Au), Au; rtol=8eps()) + @test_throws "iscover requires nonnegative scales" iscover([-1.0/L, 2.0/L], Au) + end + + @testset "a soft cover is not guaranteed to cover" begin + # The predicate's reason for existing: soft covers penalize under-coverage + # rather than forbidding it, so whether one covers is a real question. + Asoft = [1.0 4.0; 4.0 1.0] + @test iscover(symcover(Asoft), Asoft; rtol=8eps()) + @test !iscover(soft_symcover(AbsLinear{2}(), Asoft), Asoft) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index adce275..4390217 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,12 +12,13 @@ using Aqua using ExplicitImports using Test -include("helpers.jl") # iscover, covaries, PENALTIES +include("helpers.jl") # isbalanced, covaries, PENALTIES @testset "ScaleInvariantAnalysis.jl" begin include("penalties.jl") # cover_objective include("support.jl") # foreach_support(_sym) traversal + include("iscover.jl") # the cover predicate include("heuristic_covers.jl") # symcover/cover and their internals include("soft_covers.jl") # soft_symcover/soft_cover multistart descent include("initializers.jl") # initialize_symcover/initialize_cover strategies