Add iscover, the cover predicate

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")

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

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)
 

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

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