Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Squash commits to prepare for v2 #14

Merged
merged 1 commit into from
Oct 4, 2024
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
4 changes: 3 additions & 1 deletion .github/workflows/CI.yml
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,8 @@ jobs:
- uses: julia-actions/julia-buildpkg@v1
- uses: julia-actions/julia-runtest@v1
- uses: julia-actions/julia-processcoverage@v1
- uses: codecov/codecov-action@v3
- uses: codecov/codecov-action@v4
with:
files: lcov.info
token: ${{ secrets.CODECOV_TOKEN }}
fail_ci_if_error: false
91 changes: 45 additions & 46 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -10,35 +10,50 @@
`IsApprox` implements an interface for applying different definitions of "approximate" in tests for approximate (or exact) equality.
It is also fun and hip.

### Note: API change for v2

Methods in this package extend predicate functions defined in `Base` and `LinearAlgebra`.
In addition, a few predicate functions that do not have an analogue in `Base` or `LinearAlgebra`
are exported. All predicate methods defined here take an argument of type `AbstractApprox` in the final
position.

Previously, `Base.isapprox` was extended by taking an `AbstractApprox` in the first position.
All other predicate functions were duplicated or new in this package and could conflict with
names in `Base` and `LinearAlgebra`. With the new API, `using IsApprox` causes no conflicts
or piracy and does not interfere with any exsisting code using predicates, `isapprox`,
`issymmetric`, `ishermitian`, etc.

Design requirements of `IsApprox` are:

* It should provide a drop-in replacement for (and extend) `isapprox` as well as
several application functions, such as `isone` and `issymmetric`. In
particular, many functions that currently check for a property exactly (to
machine precision) will instead use `IsApprox` to implement both exact and
approximate comparison.
* It extends `Base.isapprox` as well as several application functions, such as `isone` and
`issymmetric`. In particular, many functions that currently check for a property exactly
(to machine precision) will instead use `IsApprox` to implement both exact and approximate
comparison.

* Replacements of existing methods (eg. `isone(::Float64)`) must incur no run-time penalty.
In practice, this means specifying the notion of "approximate" via types, eg `Equal` and `Approx`
so that the compiler inlines the comparison code.
* It should incur no run-time performance penalty for existing calls of `isapprox` and other
predicate functions, such as `isone`.

## Examples


The Jupyter notebook is out of date.

See this [Jupyter notebook](https://github.com/jlapeyre/IsApprox.jl/blob/master/Notebooks/IsApprox.ipynb)
for examples. See also the [test suite](https://github.com/jlapeyre/IsApprox.jl/blob/master/test/runtests.jl).
for examples.

See also the [test suite](https://github.com/jlapeyre/IsApprox.jl/blob/master/test/runtests.jl).

## Motivation

For some applications, `LinearAlgebra` wants to know if a matrix is exactly
Hermitian. Quantum information packages, on the other hand, might want to know
if a matrix is approximately (or exactly) Hermitian. Furthermore, many functions that check
whether a property (approximately) holds are interdependent. For example
`isdiag` calls functions that eventually call `iszero`. And `isposdef` calls
`ishermitian`. Furthermore again, one might want to check approximate equality
in norm; or elementwise. One might want to specify a tolerance and have it
propagate. In practice, packages
tend to reimplement tests in ways that do not satisfy all these criteria,
and fail to be composable.
Hermitian. Quantum information packages, on the other hand, might want to know if a matrix
is approximately (or exactly) Hermitian. Furthermore, many functions that check whether a
property (approximately) holds are interdependent. For example `isdiag` calls functions
that eventually call `iszero`. And `isposdef` calls `ishermitian`. Furthermore again, one
might want to check approximate equality in norm; or elementwise. One might want to
specify a tolerance and have it propagate. In practice, packages tend to reimplement tests
in ways that do not satisfy all these criteria, and fail to be composable.

Such packages include [`QuantumInformation`](https://github.com/iitis/QuantumInformation.jl)(
[code example](https://github.com/iitis/QuantumInformation.jl/blob/b47400ebb09d10cc1eba5f7bf06badeb6cfe5429/src/utils.jl#L93-L113))
,
Expand All @@ -52,12 +67,12 @@ equality is needed.

## Description

`IsApprox` allows users to specify different definitions of
closeness, via a zero-cost abstraction. That is, specifying the definition of closeness
need not incur a run-time cost. The code that implements tests for properties such as
symmetry or positivity may then be somewhat decoupled from the specification of
closeness. Furthermore, a simple, small, collection of closeness measures should be
adequate for the vast majority of use cases.
`IsApprox` allows users to specify different definitions of closeness, via a zero-cost
abstraction. That is, specifying the definition of closeness need not incur a run-time
cost. The code that implements tests for properties such as symmetry or positivity may
then be somewhat decoupled from the specification of closeness. Furthermore, a simple,
small, collection of closeness measures should be adequate for the vast majority of use
cases.

Four subtypes of `AbstractApprox` are included, `Equal`, `Approx`, `EachApprox`, and `UpToPhase`.

Expand All @@ -69,21 +84,17 @@ Four subtypes of `AbstractApprox` are included, `Equal`, `Approx`, `EachApprox`,
Consider `ishermitian`.

* `ishermitian(A)` or equivalently `ishermitian(A, Equal())` demands exact equality.
This implementation and the function of the same name in `LinearAlgebra` lower to the same code.
That is, the `IsApprox` interface adds no performance penalty.
`ishermitian(A, Equal())` calls `LinearAlgebra.ishermitian(A)`.

* `ishermitian(A, Approx(kws...))` has the same semantics as `Base.isapprox`. In this
case, we test that `A` is close to Hermitian in some norm. In this case, a separate code
path is required, namely
* `ishermitian(A, Approx(kws...))` has the same semantics as `Base.isapprox`. Namely,

```julia
ishermitian(A::AbstractMatrix, approx::Approx) = isapprox(approx, A, adjoint(A))
ishermitian(A::AbstractMatrix, approx::Approx) = isapprox(A, adjoint(A), approx)
```

* `ishermitian(A, EachApprox(kws...))`. `EachApprox` specifies element-wise closeness.
If `A` is not close to Hermitian, this test is much faster than `Approx` because
only order `1` elements must be tested. This implementation shares a code path
with that for `Equal`.
only order `1` elements must be tested.

## API

Expand All @@ -93,22 +104,10 @@ with that for `Equal`.

### `isapprox`

This extends `Base.isapprox` with methods that take an initial argument of type `AbstractApprox`.
The application functions below take an optional argument of type `AbstractApprox` in the final
This extends `Base.isapprox` with methods that take a final positional argument of type `AbstractApprox`.
The application functions below take an optional argument of type `AbstractApprox` also in the final
position and (may) forward this argument to `isapprox`.

### `isone`, `iszero`, `ishermitian`, etc.

These are not exported, and do not extend the `Base` and `LinearAlgebra` functions of the same names.
They take an optional final argument of type `AbstractApprox`. They are not exported because they
would overwrite existing definitions. However, the `AbstractApprox` interface could be moved into
`Base`.

There are also functions, which *are* exported, that are in neither `Base` nor the standard library, such
as `IsApprox.isunitary`. These follow the parameter ordering and calling conventions
as `IsApprox.isone`, etc.


## Style

This package will probably try to follow the [Blue Style Guide](https://github.com/invenia/BlueStyle).
Expand Down
2 changes: 2 additions & 0 deletions src/IsApprox.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,8 @@ using Dictionaries: Dictionaries
export AbstractApprox, Equal, EachApprox, Approx, UpToPhase
export isposdef, ispossemidef, isunitary, isinvolution, isidempotent, isnormal, commutes, anticommutes
export isnormalized, isprobdist
# From LinearAlgebra
export ishermitian, issymmetric, istriu, istril, isbanded, isdiag

# This is from DictTools.jl which is not yet registered
"""
Expand Down
90 changes: 47 additions & 43 deletions src/base_applications.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,22 +4,25 @@

using LinearAlgebra: Hermitian, Symmetric, HermOrSym

import Base: isone, iszero, isreal, isinteger
import LinearAlgebra: ishermitian, issymmetric, istriu, istril, isbanded, isdiag

### isone, iszero

isone(x, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, x, one(x))
iszero(x, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, x, zero(x))
Base.isone(x, approx_test::AbstractApprox) = isapprox(x, one(x), approx_test)
Base.iszero(x, approx_test::AbstractApprox) = isapprox(x, zero(x), approx_test)
iszero(approx_test::AbstractApprox) = x -> iszero(x, approx_test)
iszero(x::AbstractArray, approx_test::AbstractApprox=Equal()) = all(iszero(approx_test), x)
iszero(x::AbstractArray, approx_test::AbstractApprox) = all(iszero(approx_test), x)
isone(x::BigInt, ::Equal) = Base.isone(x)
iszero(x::BigInt, ::Equal) = Base.iszero(x)

isone(A::StridedMatrix, approx_test::Approx) = isapprox(approx_test, A, one(A))
iszero(A::StridedMatrix, approx_test::Approx) = isapprox(approx_test, A, zero(A))
isone(A::StridedMatrix, approx_test::Approx) = isapprox(A, one(A), approx_test)
iszero(A::StridedMatrix, approx_test::Approx) = isapprox(A, zero(A), approx_test)

# dense.jl
const ISONE_CUTOFF = 2^21 # 2M

function isone(A::StridedMatrix, approx_test::AbstractApprox=Equal())
function isone(A::StridedMatrix, approx_test::AbstractApprox)
m, n = size(A)
m != n && return false # only square matrices can satisfy x == one(x)
if sizeof(A) < ISONE_CUTOFF
Expand Down Expand Up @@ -56,66 +59,67 @@

# The call `ishermitian(A::AbstractMatrix, B::AbstractMatrix) lowers
# to exactly the same code as that in LinearAlgebra
function ishermitian(A::AbstractMatrix, approx_test::AbstractApprox=Equal())
function ishermitian(A::AbstractMatrix, approx_test::AbstractApprox)
indsm, indsn = axes(A)
if indsm != indsn
return false
end
for i = indsn, j = i:last(indsn)
if ! isapprox(approx_test, A[i,j], adjoint(A[j,i]))
if ! isapprox(A[i,j], adjoint(A[j,i]), approx_test)
return false
end
end
return true
end

# This uses the isapprox interface, which compares using a norm
ishermitian(A::AbstractMatrix, approx_test::Approx) = isapprox(approx_test, A, adjoint(A))
ishermitian(x::Number, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, x, conj(x))
ishermitian(A::Hermitian, ::AbstractApprox=Equal()) = true
# This method uses the isapprox interface, which compares using the Frobenius norm.
ishermitian(A::AbstractMatrix, approx_test::Approx) = isapprox(A, adjoint(A), approx_test)

ishermitian(x::Number, approx_test::AbstractApprox) = isapprox(x, conj(x), approx_test)
ishermitian(A::Hermitian, ::AbstractApprox) = true
ishermitian(A::Hermitian, ::Approx) = true
ishermitian(A::Symmetric{<:Real}, ::AbstractApprox=Equal()) = true
ishermitian(A::Symmetric{<:Real}, ::AbstractApprox) = true
ishermitian(A::Symmetric{<:Real}, ::Approx) = true
ishermitian(A::Symmetric{<:Complex}, approx_test::AbstractApprox=Equal()) = isreal(A, approx_test)
ishermitian(A::Symmetric{<:Complex}, approx_test::AbstractApprox) = isreal(A, approx_test)

Check warning on line 83 in src/base_applications.jl

View check run for this annotation

Codecov / codecov/patch

src/base_applications.jl#L83

Added line #L83 was not covered by tests
ishermitian(A::Symmetric{<:Complex}, approx_test::Approx) = isreal(A, approx_test)
issymmetric(A::Hermitian{<:Real}, ::AbstractApprox=Equal()) = true
issymmetric(A::Hermitian{<:Complex}, approx_test::AbstractApprox=Equal()) = isreal(A, approx_test)
issymmetric(A::Symmetric, ::AbstractApprox=Equal()) = true
issymmetric(A::Hermitian{<:Real}, ::AbstractApprox) = true
issymmetric(A::Hermitian{<:Complex}, approx_test::AbstractApprox) = isreal(A, approx_test)

Check warning on line 86 in src/base_applications.jl

View check run for this annotation

Codecov / codecov/patch

src/base_applications.jl#L85-L86

Added lines #L85 - L86 were not covered by tests
issymmetric(A::Symmetric, ::AbstractApprox) = true

issymmetric(A::AbstractMatrix{<:Real}, approx_test::AbstractApprox=Equal()) =
issymmetric(A::AbstractMatrix{<:Real}, approx_test::AbstractApprox) =
ishermitian(A, approx_test)

# Copied from LinearAlgebra. Why does the iteration over i differ slightly from ishermitian ?
function issymmetric(A::AbstractMatrix, approx_test::AbstractApprox=Equal())
function issymmetric(A::AbstractMatrix, approx_test::AbstractApprox)
indsm, indsn = axes(A)
if indsm != indsn
return false
end
for i = first(indsn):last(indsn), j = (i):last(indsn)
if ! isapprox(approx_test, A[i,j], transpose(A[j,i]))
if ! isapprox(A[i,j], transpose(A[j,i]), approx_test)
return false
end
end
return true
end

issymmetric(x::Number, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, x, x)
issymmetric(x::Number, approx_test::AbstractApprox) = isapprox(x, x, approx_test)

### isreal

# complex.jl
isreal(x::Real, approx_test::AbstractApprox=Equal()) = true
isreal(z::Complex, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, real(z), z)
isreal(x::Real, approx_test::AbstractApprox) = true
isreal(z::Complex, approx_test::AbstractApprox) = isapprox(real(z), z, approx_test)
# Old way
#isreal(z::Complex, approx_test::AbstractApprox=Equal()) = iszero(imag(z), approx_test)
#isreal(z::Complex, approx_test::AbstractApprox) = iszero(imag(z), approx_test)

isreal(x::AbstractArray{<:Real}, approx_test::AbstractApprox=Equal()) = true
isreal(x::AbstractArray{<:Real}, approx_test::AbstractApprox) = true

isreal(approx_test::AbstractApprox) = x -> isreal(x, approx_test)
isreal(x::AbstractArray, approx_test::AbstractApprox=Equal()) = all(isreal(approx_test),x)
isreal(x::AbstractArray, approx_test::AbstractApprox) = all(isreal(approx_test),x)

isreal(A::HermOrSym{<:Real}, approx_test::AbstractApprox=Equal()) = true
function isreal(A::HermOrSym, approx_test::AbstractApprox=Equal())
isreal(A::HermOrSym{<:Real}, approx_test::AbstractApprox) = true
function isreal(A::HermOrSym, approx_test::AbstractApprox)
n = size(A, 1)
@inbounds if A.uplo == 'U'
for j in 1:n
Expand All @@ -140,28 +144,28 @@
### isinteger

# This must be changed if this is integrated into Base
isinteger(x) = isinteger(x, Equal())
# isinteger(x) = isinteger(x, Equal())
# We use union or explicit types to avoid method ambiguity. Is there a way around this ?
isinteger(x::BigFloat, ::Equal) = Base.isinteger(x)
isinteger(x::Rational, ::Equal) = Base.isinteger(x)

# number.jl
isinteger(x::Integer, ::AbstractApprox=Equal()) = true
isinteger(x::Integer, ::AbstractApprox) = true

# floatfuncs.jl
## The original is x - trunc(x) == 0. So this implementation might differ
## from the stock version of isinteger for some user's subtype of AbstractFloat.
## We choose to this implementation because the default relative tolerance is
## reasonable. That is, `isapprox(approx_test, x - trunc(x), 0)` requires
## specifying `atol`.
isinteger(x::AbstractFloat, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, x, trunc(x))
isinteger(x::AbstractFloat, approx_test::AbstractApprox) = isapprox(x, trunc(x), approx_test)

# complex.jl
isinteger(z::Complex, approx_test::AbstractApprox=Equal()) =
isinteger(z::Complex, approx_test::AbstractApprox) =
isreal(z, approx_test) && isinteger(real(z), approx_test)

# TODO: Need to think about difference between real(z) and abs(z) regarding tolerance in all
# methods for complex numbers (not just isingeger)
# methods for complex numbers (not just isinteger)
isinteger(z::Complex, approx_test::UpToPhase) = isinteger(abs(z), approx_test)

isinteger(x::Rational, approx_test::AbstractApprox) = isinteger(float(x), approx_test)
Expand All @@ -172,7 +176,7 @@
_require_one_based_indexing(A...) = !Base.has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))

# TODO: Why is approx a kw arg here and below?
function istriu(A::AbstractMatrix, k::Integer = 0; approx::AbstractApprox=Equal())
function istriu(A::AbstractMatrix, k::Integer, approx::AbstractApprox)
_require_one_based_indexing(A)
m, n = size(A)
for j in 1:min(n, m + k - 1)
Expand All @@ -182,9 +186,9 @@
end
return true
end
istriu(x::Number, ::AbstractApprox) = true
istriu(::Number, ::AbstractApprox) = true

function istril(A::AbstractMatrix, k::Integer = 0; approx::AbstractApprox=Equal())
function istril(A::AbstractMatrix, k::Integer, approx::AbstractApprox)
_require_one_based_indexing(A)
m, n = size(A)
for j in max(1, k + 2):n
Expand All @@ -194,13 +198,13 @@
end
return true
end
istril(x::Number, ::AbstractApprox) = true
istril(::Number, ::AbstractApprox) = true

isbanded(A::AbstractMatrix, kl::Integer, ku::Integer; approx::AbstractApprox=Equal()) =
istriu(A, kl; approx=approx) && istril(A, ku; approx=approx)
isbanded(A::AbstractMatrix, kl::Integer, ku::Integer, approx::AbstractApprox) =
istriu(A, kl, approx) && istril(A, ku, approx)

isdiag(A::AbstractMatrix; approx::AbstractApprox=Equal()) = isbanded(A, 0, 0; approx=approx)
isdiag(x::Number; approx::AbstractApprox=Equal()) = true
isdiag(A::HermOrSym; approx::AbstractApprox=Equal()) =
isdiag(A::AbstractMatrix, approx::AbstractApprox) = isbanded(A, 0, 0, approx)
isdiag(x::Number, approx::AbstractApprox) = true
isdiag(A::HermOrSym, approx::AbstractApprox) =
isdiag(A.uplo == 'U' ? LinearAlgebra.UpperTriangular(A.data) :
LinearAlgebra.LowerTriangular(A.data); approx=approx)
LinearAlgebra.LowerTriangular(A.data), approx)
Loading
Loading