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