Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions src/ScaleInvariantAnalysis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
30 changes: 30 additions & 0 deletions src/penalties.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

"""
Expand Down
74 changes: 65 additions & 9 deletions src/support.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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))"))
Expand Down