From 9418f506e8e2932018959dcb03c7f9f4d6792076 Mon Sep 17 00:00:00 2001 From: John Lapeyre Date: Sun, 15 Sep 2024 21:48:28 -0400 Subject: [PATCH] Squash commits to prepare for v2 The main change in the API is that functions in `IsApprox` with the same name as those in Base and `LinearAlgebra` are now different methods of the same function. So you can do `using` of both `LinearAlgebra` and `IsApprox` withouth qualifying any symbols. Another API change making the signatures uniform, rather than having the same argument in different places in different functions or methods. Making this change was enabled by the previous item. * A few details * Change argument order for Base.isapprox (#13) Previously, I had `Base.isapprox(::AbstractApprox, x, y)`. Now we use `Base.isapprox(x, y, ::AbstractApprox)`. This order agrees with all other testing functions included in this package, such as, `isone`, `isreal`, `ishermitian`, etc. * All predicates either extend `Base`/`LinearAlgebra` or are new to this package Previously, some predicates were copies of functions in `Base`/`LinearAlgebra`, but did not extend these functions with new methods, but rather were new functions under `IsApprox`. Other predicates in `IsApprox` had no function in `Base`/`LinearAlgebra` of the same name. This commit changes this so that all predicate functions of the same name as those in `Base`/`LinearAlgebra` now extend those functions. All predicate functions, including those extended from `LinearAlgebra` are exported. With this commit, you can do `using IsApprox` and suffer no name collisions with `Base`/`LinearAlgebra`. `JET` and `Aqua` find no ambiguities or piracy. * This commit and the previous, precursor commit, include a breaking change to the API. Methods like `Base.isapprox(::AbstractApprox, x, y)` are now `Base.isapprox(x, y, ::AbstractApprox)` Methods like `Base.isapprox(x, y; approx::AbstractApprox)` are now `Base.isapprox(x, y, ::AbstractApprox)` In other words, the calling convention is now uniform. * Other changes * Add tests and fix small bug * Use 42.0 rather than 42 in isinteger test A method for rationalize(::Int) does not exist in older versions (at leat 1.6) of Julia. So we use a float instead so that the test runs on more versions. * Fix bug in isunitary for Equal and IsApprox Don't compare result to zero. Add one and compare to one. This is the proper scale because we expect the vectors to be normalized. * Explain new API * Fix some incorrect statements in comments `Base.isapprox` does not use matrix norms. It treats multidimensional arrays as vectors. Closes #12 * Try updating codecov action Added repository secret to GH repo config as well. * Disable test using feature not in v1.0 * Add tests, remove two nonsense methods Looks like I tried to add methods to satisfy JET or Aqua in the past. But they don't have arg combinations that could ever be called, even because of ambiguity. I removed them and all tests, including JET and Aqua still pass. --- .github/workflows/CI.yml | 4 +- README.md | 91 ++++++++-------- src/IsApprox.jl | 2 + src/base_applications.jl | 90 ++++++++-------- src/core.jl | 35 +++--- src/other_applications.jl | 98 ++++++++++------- src/precompile.jl | 2 +- test/runtests.jl | 216 ++++++++++++++++++++++++++++++++++++-- 8 files changed, 380 insertions(+), 158 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 85eaa12..5619bba 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -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 diff --git a/README.md b/README.md index 8033355..cd2877b 100644 --- a/README.md +++ b/README.md @@ -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)) , @@ -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`. @@ -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 @@ -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). diff --git a/src/IsApprox.jl b/src/IsApprox.jl index 59536e1..13c2e68 100644 --- a/src/IsApprox.jl +++ b/src/IsApprox.jl @@ -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 """ diff --git a/src/base_applications.jl b/src/base_applications.jl index 5a86ff6..9f865c1 100644 --- a/src/base_applications.jl +++ b/src/base_applications.jl @@ -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 @@ -56,66 +59,67 @@ end # 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) 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) +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 @@ -140,13 +144,13 @@ end ### 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 @@ -154,14 +158,14 @@ isinteger(x::Integer, ::AbstractApprox=Equal()) = true ## 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) @@ -172,7 +176,7 @@ isinteger(x::Rational, approx_test::AbstractApprox) = isinteger(float(x), approx _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) @@ -182,9 +186,9 @@ function istriu(A::AbstractMatrix, k::Integer = 0; approx::AbstractApprox=Equal( 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 @@ -194,13 +198,13 @@ function istril(A::AbstractMatrix, k::Integer = 0; approx::AbstractApprox=Equal( 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) diff --git a/src/core.jl b/src/core.jl index 400e795..cb1d385 100644 --- a/src/core.jl +++ b/src/core.jl @@ -37,10 +37,9 @@ EachApprox(; kws...) = EachApprox(kws) """ Approx(; kws...) <: AbstractApprox -Specifies using the legacy `isapprox` interface. For example, -for `AbstractMatrix`, matrix norms are used to test closeness. -`kw` are keyword pairs that are forwarded to `isapprox`. -For example, `Approx(atol=1e-9)`. +Specifies using the legacy `isapprox` interface. For example, for +`AbstractMatrix`, vector norms are used to test closeness. `kw` are keyword +pairs that are forwarded to `isapprox`. For example, `Approx(atol=1e-9)`. """ struct Approx{T} <: AbstractApprox kw::T @@ -52,33 +51,33 @@ end Approx(; kws...) = Approx(kws) """ - isapprox(::Equal, x, y) + isapprox(x, y, ::Equal) Return `true` if `x == y`. """ -Base.isapprox(::Equal, x, y) = (x == y) +Base.isapprox(x, y, ::Equal) = (x == y) """ - isapprox(a::Union{EachApprox, Approx}, x, y) + isapprox(x, y, a::Union{EachApprox, Approx}) Use the definition of approximate equality specified by `a` to determine if `x` is approximately `y`. """ -Base.isapprox(a::Union{EachApprox, Approx}, x, y) = isapprox(x, y; pairs(a.kw)...) +Base.isapprox(x, y, a::Union{EachApprox, Approx}) = isapprox(x, y; a.kw...) """ - isapprox(a::EachApprox, A::AbstractArray, B::AbstractArray) + isapprox(A::AbstractArray, B::AbstractArray, a::EachApprox) Compute element-wise approximate equality. """ -function Base.isapprox(a::EachApprox, A::AbstractArray, B::AbstractArray) +function Base.isapprox(A::AbstractArray, B::AbstractArray, approx::EachApprox) n1 = length(A) n2 = length(B) if n1 != n2 return false end for (x, y) in zip(A, B) - if ! isapprox(a, x, y) + if ! isapprox(x, y, approx) return false end end @@ -100,17 +99,17 @@ struct UpToPhase{T} <: AbstractApprox end UpToPhase(; kws...) = UpToPhase(kws) -function Base.isapprox(a::UpToPhase, x::Number, y::Number) +function Base.isapprox(x::Number, y::Number, a::UpToPhase) aa = Approx(;a.kw...) - if isapprox(aa, x, zero(x)) - return isapprox(aa, y, zero(y)) - elseif isapprox(aa, y, zero(y)) - return isapprox(aa, x, zero(x)) + if isapprox(x, zero(x), aa) + return isapprox(y, zero(y), aa) + elseif isapprox(y, zero(y), aa) + return isapprox(x, zero(x), aa) end return isunitary(x / y, aa) end -function Base.isapprox(_app::UpToPhase, A::AbstractArray, B::AbstractArray) +function Base.isapprox(A::AbstractArray, B::AbstractArray, _app::UpToPhase) n1 = length(A) n2 = length(B) if n1 != n2 @@ -130,7 +129,7 @@ function Base.isapprox(_app::UpToPhase, A::AbstractArray, B::AbstractArray) isunitary(z, app) || return false seen_non_zero_flag = true else - isapprox(app, a/b, z) || return false + isapprox(a/b, z, app) || return false end end end diff --git a/src/other_applications.jl b/src/other_applications.jl index 9a07255..0fc91c8 100644 --- a/src/other_applications.jl +++ b/src/other_applications.jl @@ -3,12 +3,16 @@ import LinearAlgebra ## We use ispossemidef for numbers as well as matrices. ## This follows `ishermitian`, etc. +import LinearAlgebra: isposdef + +ispossemidef(x) = ispossemidef(x, Equal()) + """ - ispossemidef(x::Number, approx_test::AbstractApprox=Equal()) + ispossemidef(x::Number, approx_test::AbstractApprox) Return `true` if `x` is (approximately) non-negative. """ -function ispossemidef(x::Real, approx_test::AbstractApprox=Equal()) +function ispossemidef(x::Real, approx_test::AbstractApprox) return x > zero(x) || iszero(x, approx_test) end @@ -17,21 +21,21 @@ function _isposdef(z::Complex, approx_test::AbstractApprox, posdeffunc) return posdeffunc(real(z), approx_test) && iszero(imag(z), approx_test) end -function ispossemidef(z::Complex, approx_test::AbstractApprox=Equal()) +function ispossemidef(z::Complex, approx_test::AbstractApprox) return _isposdef(z, approx_test, ispossemidef) end -function isposdef(z::Complex, approx_test::AbstractApprox=Equal()) +function isposdef(z::Complex, approx_test::AbstractApprox) return _isposdef(z, approx_test, isposdef) end """ - isposdef(x::Number, approx_test::AbstractApprox=Equal()) + isposdef(x::Number, approx_test::AbstractApprox) Return `true` if `x` is approximately greater than zero. """ -isposdef(x::Number) = isposdef(x, Equal()) isposdef(x::Real, ::Equal) = x > zero(x) +#isposdef(x::Number) = isposdef(x, Equal()) ## For methods other than `Equal`, x can be zero or negative isposdef(x::Number, approx_test::AbstractApprox) = ispossemidef(x, approx_test) @@ -43,16 +47,16 @@ function _isposdef(m::AbstractMatrix, approx_test::AbstractApprox, posdeffunc) end """ - ispossemidef(m::AbstractMatrix, approx_test::AbstractApprox=Equal()) + ispossemidef(m::AbstractMatrix, approx_test::AbstractApprox) Return `true` if `m` is positive semidefinite. """ -function ispossemidef(m::AbstractMatrix, approx_test::AbstractApprox=Equal()) +function ispossemidef(m::AbstractMatrix, approx_test::AbstractApprox) return _isposdef(m, approx_test, ispossemidef) end """ - isposdef(m::AbstractMatrix, approx_test::AbstractApprox=Equal()) + isposdef(m::AbstractMatrix, approx_test::AbstractApprox) Return `true` if `m` is positive definite. """ @@ -60,13 +64,13 @@ function isposdef(m::AbstractMatrix, approx_test::AbstractApprox) return _isposdef(m, approx_test, isposdef) end -isposdef(A::AbstractMatrix) = isposdef(A, Equal()) +# TODO: I think I no longer need this, after chaning arg order, +# and importing all predicate functions +#isposdef(A::AbstractMatrix) = isposdef(A, Equal()) # copied from dense.jl isposdef(A::AbstractMatrix, ::Equal) = ishermitian(A, Equal()) && LinearAlgebra.isposdef(LinearAlgebra.cholesky(LinearAlgebra.Hermitian(A); check = false)) - - ## Compared two methods: ## a) Allocate, ie m' * m. b) iterate over columns ## Found: @@ -75,35 +79,38 @@ isposdef(A::AbstractMatrix, ::Equal) = ## 3. Doing the allocation, eg m' * m can be slightly faster. Eg for 100x100 dense identity matrix. ## 4. For rand(100, 100), iterating over columns is 1000 times faster. Fails on first column. ## `approx_test` is `Equal` or `EachApprox`. -function _isunitary(m::AbstractMatrix, approx_test::AbstractApprox, dotf, transposef) +function _isunitary(m::AbstractMatrix, approx_test::AbstractApprox, dotf::F1, transposef::F2) where {F1, F2} + _one = one(eltype(m)) rowinds = axes(m)[2] for i in rowinds - isapprox(approx_test, dotf(view(m, :, i), view(transposef(m), :, i)), 1) || return false + isapprox(dotf(view(m, :, i), view(transposef(m), :, i)), _one, approx_test) || return false for j in i+1:last(rowinds) - isapprox(approx_test, dotf(view(m, :, i), view(transposef(m), :, j)), 0) || return false + isapprox(dotf(view(m, :, i), view(transposef(m), :, j)) + _one, _one, approx_test) || return false end end return true end -## This is slower even for small matrices. -## Then why am I using it ? (May 2021) +isunitary(x) = isunitary(x, Equal()) + +## Use vector norm. +## Slower, but more generally useful. function isunitary(m::AbstractMatrix, approx_test::Approx) - return isapprox(approx_test, m' * m, LinearAlgebra.I) + return isapprox(m' * m, LinearAlgebra.I, approx_test) end _identity(x) = x """ - isunitary(m::AbstractMatrix, approx_test::AbstractApprox=Equal()) + isunitary(m::AbstractMatrix, approx_test::AbstractApprox) Return `true` if `m` is unitary. If `m` is real, this tests orthogonality. """ -isunitary(m::AbstractMatrix, approx_test::AbstractApprox=Equal()) = +isunitary(m::AbstractMatrix, approx_test::AbstractApprox) = _isunitary(m, approx_test, LinearAlgebra.dot, _identity) # abs2 is much faster, but we would need to use sqrt to adjust the tolerance, thus losing any advantage. -isunitary(x::Number, approx_test::AbstractApprox=Equal()) = isone(abs(x), approx_test) -isunitary(J::LinearAlgebra.UniformScaling, approx_test::AbstractApprox=Equal()) = isunitary(J.λ, approx_test) +isunitary(x::Number, approx_test::AbstractApprox) = isone(abs(x), approx_test) +isunitary(J::LinearAlgebra.UniformScaling, approx_test::AbstractApprox) = isunitary(J.λ, approx_test) """ _dotu(x::AbstractVector, y::AbstractVector) @@ -121,51 +128,62 @@ function _dot(x::AbstractVector{<:Complex}, y::AbstractVector{<:Complex}) return sum(*(z...) for z in zip(x, y)) end +isinvolution(x) = isinvolution(x, Equal()) + """ - isinvolution(m::AbstractMatrix, approx_test::AbstractApprox=Equal()) + isinvolution(m::AbstractMatrix, approx_test::AbstractApprox) Return `true` if `m * m == I` """ -isinvolution(m::AbstractMatrix, approx_test::AbstractApprox=Equal()) = _isunitary(m, approx_test, _dot, transpose) +isinvolution(m::AbstractMatrix, approx_test::AbstractApprox) = _isunitary(m, approx_test, _dot, transpose) function isinvolution(m::AbstractMatrix, approx_test::Approx) - return isapprox(approx_test, m * m, LinearAlgebra.I) + return isapprox(m * m, LinearAlgebra.I, approx_test) end -function isidempotent(m::AbstractMatrix, approx_test::AbstractApprox=Equal()) - return isapprox(approx_test, m * m, m) +isidempotent(x) = isidempotent(x, Equal()) + +function isidempotent(m::AbstractMatrix, approx_test::AbstractApprox) + return isapprox(m * m, m, approx_test) end -function isnormal(m::AbstractMatrix, approx_test::AbstractApprox=Equal()) - return isapprox(approx_test, m * m', m' * m) +isnormal(x) = isnormal(x, Equal()) + +function isnormal(m::AbstractMatrix, approx_test::AbstractApprox) + return isapprox(m * m', m' * m, approx_test) end """ - commutes(X, Y, approx_test::AbstractApprox=Equal()) + commutes(X, Y, approx_test::AbstractApprox) Return `true` if `X` and `Y` commute. """ -commutes(X, Y, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, X * Y, Y * X) +commutes(X, Y, approx_test::AbstractApprox) = isapprox(X * Y, Y * X, approx_test) + +commutes(X, Y) = commutes(X, Y, Equal()) """ - anticommutes(X, Y, approx_test::AbstractApprox=Equal()) + anticommutes(X, Y, approx_test::AbstractApprox) Return `true` if `X` and `Y` anticommute. """ -anticommutes(X, Y, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, X * Y, -(Y * X)) +anticommutes(X, Y, approx_test::AbstractApprox) = isapprox(X * Y, -(Y * X), approx_test) + +anticommutes(X, Y) = anticommutes(X, Y, Equal()) """ - isnormalized(itr, approx_test::AbstractApprox=Equal()) + isnormalized(itr, approx_test::AbstractApprox) Return `true` if `itr` is normalized, that is, if the items sum to one. If `itr` is a dictionary, the items are the values. """ -isnormalized(itr, approx_test::AbstractApprox=Equal()) = isnormalized(Base.IteratorEltype(itr), itr, approx_test) -isnormalized(::Base.EltypeUnknown, itr, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, sum(itr), 1) -isnormalized(::Base.HasEltype, itr, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, sum(itr), one(eltype(itr))) -isnormalized(d::_AbstractDict{<:Any, V}, approx_test::AbstractApprox=Equal()) where V = isapprox(approx_test, sum(values(d)), one(V)) -isnormalized(x::Base.HasEltype, approx_test::AbstractApprox) = throw(MethodError(isnormalized, (x, approx_test))) -isnormalized(x::Base.EltypeUnknown, approx_test::AbstractApprox) = throw(MethodError(isnormalized, (x, approx_test))) +isnormalized(itr, approx_test::AbstractApprox) = isnormalized(Base.IteratorEltype(itr), itr, approx_test) +isnormalized(::Base.EltypeUnknown, itr, approx_test::AbstractApprox) = isapprox(sum(itr), 1, approx_test) +isnormalized(::Base.HasEltype, itr, approx_test::AbstractApprox) = isapprox(sum(itr), one(eltype(itr)), approx_test) +isnormalized(d::_AbstractDict{<:Any, V}, approx_test::AbstractApprox) where V = isapprox(sum(values(d)), one(V), approx_test) +# isnormalized(x::Base.HasEltype, approx_test::AbstractApprox) = throw(MethodError(isnormalized, (x, approx_test))) +# isnormalized(x::Base.EltypeUnknown, approx_test::AbstractApprox) = throw(MethodError(isnormalized, (x, approx_test))) +isnormalized(x) = isnormalized(x, Equal()) # A Vector has no algebraic interpretation wrt isposdef, ispossemidef. So we don't define # a method. To help isprobdist, we do the following. diff --git a/src/precompile.jl b/src/precompile.jl index 1466146..dc7f163 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -13,7 +13,7 @@ IsApprox.iszero(m, approx) IsApprox.isone(m, approx) IsApprox.isreal(m, approx) - IsApprox.isdiag(m, approx=approx) + IsApprox.isdiag(m, approx) IsApprox.isunitary(m, approx) IsApprox.isidempotent(m, approx) IsApprox.isnormalized(m, approx) diff --git a/test/runtests.jl b/test/runtests.jl index feda169..dd874c3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using IsApprox: IsApprox, isone, iszero, isreal, isinteger, ispossemidef, isposd using Dictionaries: Dictionary using Test import LinearAlgebra +using LinearAlgebra: Hermitian, Symmetric if VERSION >= v"1.7" && VERSION <= v"1.11" @testset "JET" begin @@ -12,6 +13,18 @@ end include("aqua_test.jl") +@testset "isapprox" begin + m = rand(2,2) + for atest = (EachApprox(), UpToPhase()) + wrong_shape = rand(2, 3) + @test ! isapprox(m, wrong_shape, atest) + @test isapprox(m, m, atest) + end + @test isapprox(0.0, 0.0, UpToPhase()) + @test !isapprox(0.0, 1.0, UpToPhase()) + @test !isapprox(1.0, 0.0, UpToPhase()) +end + @testset "isnormalized, isprobdist" begin v1 = collect(1:10) ./ sum(1:10) g1 = (x for x in v1) @@ -37,6 +50,8 @@ include("aqua_test.jl") @test ! isprobdist(v4) @test isprobdist(v4, Approx()) @test ! isprobdist(v4, Approx(atol=1e-16)) + + @test_throws MethodError isnormalized(["dog"], Approx()) end @@ -51,6 +66,22 @@ end @test isone(1.0 + 1e-7, Approx(atol=1e-6)) @test ! isone(1.0 + 1e-7, Approx()) + @test ! isone(rand(600, 600), Approx()) + + if VERSION >= v"1.6" + bigi = collect(float(LinearAlgebra.I(600))) + @test isone(bigi, Approx()) + bigi_noise = bigi + 1e-14 * randn(600,600) + @test isone(bigi_noise, Approx()) + @test isone(bigi_noise, EachApprox(atol=1e-12)) + @test !isone(bigi_noise, Approx(atol=1e-15)) + @test !isone(bigi_noise, EachApprox(atol=1e-15)) + end + + m = [1. 0.; 1e-14 1.] + @test isone(m, Approx()) + @test !isone(m, Approx(atol=1e-16)) + # maximum element is 9.9e-11 # norm(m) == 1.8892148837710255e-10 m = [9.914830625828892e-11 4.473996666408231e-11 3.066184129948095e-11; @@ -66,12 +97,28 @@ end end @testset "isreal" begin + isreal = IsApprox.isreal + @test isreal(1) @test isreal(1.0) @test isreal(1.0, Approx()) @test ! isreal(1.0 + 1e-10im) @test ! isreal(1.0 + 1e-7im, Approx()) @test isreal(1.0 + 1e-10im, Approx(atol=1e-9)) + + m = real.(rand(ComplexF64, 3, 3)) + @test isreal(m) + m_noisy = m + 1e-14 * randn(ComplexF64, 3, 3) + @test !isreal(m_noisy) + + for approx in (Approx, EachApprox) + @test isreal(m_noisy, approx()) + @test isreal(m_noisy, approx(atol=1e-8)) + @test !isreal(m_noisy, approx(atol=1e-16)) + @test isreal(m_noisy, approx()) + @test isreal(m_noisy, approx(atol=1e-8)) + @test !isreal(m_noisy, approx(atol=1e-16)) + end end @testset "isinteger" begin @@ -84,6 +131,13 @@ end @test ! isinteger(1 + im * 1e-5, Approx()) @test isinteger(exp(im*2), UpToPhase()) @test ! isinteger(1.1 * exp(im*2), UpToPhase()) + @test isinteger(BigFloat(10), Equal()) + @test !isinteger(BigFloat(10//12), Equal()) + @test isinteger(10//1, Equal()) + @test !isinteger(10//12, Equal()) + @test !isinteger(10//12, Approx()) + @test isinteger(1000000001//1000000000, Approx()) + @test !isinteger(1000000001//1000000000, Equal()) end @testset "isdiag" begin @@ -91,7 +145,8 @@ end m = [1.0 6.108385298833888e-20; 7.691926633708195e-20 1.0] @test isdiag(m0) @test ! isdiag(m) - @test isdiag(m; approx=Approx(atol=1e-10)) + @test isdiag(m, Approx(atol=1e-10)) + @test isdiag(3, Approx()) end @testset "isposdef" begin @@ -106,6 +161,9 @@ end @test ispossemidef(1) @test ! ispossemidef(-1) @test isposdef(1 + 0im) + @test ispossemidef(1 + 0im) + @test ispossemidef(1. + 0.0 * im) + @test ispossemidef(1. + im*1e-14, Approx(atol=1e-10)) @test ispossemidef(-1e-5, Approx(atol=1e-3)) @test isposdef(1) @@ -150,7 +208,8 @@ end @test ! isunitary(merr) @test ! isunitary(merr, Equal()) @test isunitary(merr, Approx()) - @test ! isunitary(merr, EachApprox()) + @test isunitary(merr, EachApprox()) + @test ! isunitary(merr, EachApprox(rtol=1e-20)) @test ! isunitary(merr, EachApprox(atol=1e-10)) @test isunitary(merr, EachApprox(atol=5e-10)) @test ! isunitary(merr, Approx(atol=1e-10)) @@ -175,14 +234,20 @@ end end @testset "UpToPhase" begin - @test isapprox(UpToPhase(), 2, 2) - @test isapprox(UpToPhase(), 2.1, 2.1 + 1e-10) - @test ! isapprox(UpToPhase(atol=1e-12), 2.1, 2.1 + 1e-10) - @test isapprox(UpToPhase(), 2.1, 2.1 * cis(2*pi*3.1)) + @test isapprox(2, 2, UpToPhase()) + @test isapprox(2.1, 2.1 + 1e-10, UpToPhase()) + @test ! isapprox(2.1, 2.1 + 1e-10, UpToPhase(atol=1e-12)) + @test isapprox(2.1, 2.1 * cis(2*pi*3.1), UpToPhase()) m = rand(5, 5) - @test isapprox(UpToPhase(), m, m .* cis(2*pi*1.3)) - @test isapprox(UpToPhase(), m, m .* (1+1e-10)*cis(2*pi*1.3)) - @test ! isapprox(UpToPhase(), m, m .* (1+1e-7)*cis(2*pi*1.3)) + @test isapprox(m, m .* cis(2*pi*1.3), UpToPhase()) + @test isapprox(m, m .* (1+1e-10)*cis(2*pi*1.3), UpToPhase()) + @test ! isapprox(m, m .* (1+1e-7)*cis(2*pi*1.3), UpToPhase()) + + m = rand(2,2) + mz = copy(m) + mz[2, 1] = 0 + @test !isapprox(m, mz, UpToPhase()) + @test !isapprox(mz, m, UpToPhase()) end @testset "Dictionary" begin @@ -203,3 +268,136 @@ end @test_throws MethodError EachApprox("dog") @test_throws MethodError UpToPhase(1e-10) end + +# These methods only exist to avoid method ambiguity +@testset "BigInt BigFloat Rational" begin + isone = IsApprox.isone + iszero = IsApprox.iszero + isinteger = IsApprox.isinteger + @test isone(big(1)) + @test !isone(big(0)) + @test !iszero(big(1)) + @test iszero(big(0)) + @test isone(big(1), Equal()) + @test !isone(big(0), Equal()) + @test !iszero(big(1), Equal()) + @test iszero(big(0), Equal()) + @test isinteger(rationalize(42.0)) + @test ! isinteger(1//2) +end + +@testset "Hermitian Symmetric" begin + ishermitian = IsApprox.ishermitian + isreal = IsApprox.isreal + issymmetric = IsApprox.issymmetric + isdiag = IsApprox.isdiag + + m = rand(ComplexF64, 4, 4); + m_herm = m + m' + @test !ishermitian(m) + @test ishermitian(m_herm) + noise = 1e-14 * randn(ComplexF64, 4, 4); + m_herm_noisy = m_herm + noise + @test ! ishermitian(m_herm_noisy) + @test ! ishermitian(m_herm_noisy, Equal()) + @test ishermitian(m_herm_noisy, Approx()) + @test ishermitian(m_herm_noisy, Approx(atol=1e-8)) + @test ! ishermitian(m_herm_noisy, Approx(atol=1e-16)) + + @test ! ishermitian(rand(2, 3), EachApprox()) + @test_throws DimensionMismatch ishermitian(rand(2, 3), Approx()) + @test ! issymmetric(rand(2, 3), EachApprox()) + + for uplo in (:U, :L) + m_Herm = Hermitian(m, uplo) + @test ishermitian(m_Herm) + @test ishermitian(m_Herm, Approx()) + @test !isreal(m_herm) + @test !isreal(m_herm, Approx()) + @test !isreal(m_Herm) + @test !isreal(m_Herm, Approx()) + @test !isdiag(m_herm) + @test !isdiag(m_Herm) + @test !isdiag(m_Herm, Approx()) + end + + m = rand(Float64, 4, 4); + m_symm = m + m' + @test !ishermitian(m) + @test ishermitian(m_symm) + noise = 1e-14 * randn(Float64, 4, 4); + m_symm_noisy = m_symm + noise + @test ! ishermitian(m_symm_noisy) + @test ! ishermitian(m_symm_noisy, Equal()) + @test ishermitian(m_symm_noisy, Approx()) + @test ishermitian(m_symm_noisy, Approx(atol=1e-8)) + @test ! ishermitian(m_symm_noisy, Approx(atol=1e-16)) + + mc = rand(ComplexF64, 4,4) + @test ishermitian(Hermitian(mc), Approx()) + @test ishermitian(Hermitian(mc), EachApprox()) + @test issymmetric(Symmetric(mc), Approx()) + @test issymmetric(Symmetric(mc), EachApprox()) + mr = rand(4,4) + @test ishermitian(Hermitian(mr), Approx()) + @test ishermitian(Hermitian(mr), EachApprox()) + @test issymmetric(Symmetric(mr), Approx()) + @test issymmetric(Symmetric(mr), EachApprox()) + m_Symm = Symmetric(m) + @test ishermitian(m_Symm) + @test ishermitian(m_Symm, Approx()) + @test ishermitian(m_Symm, EachApprox()) + + @test issymmetric(1.0, Approx()) + + @test ishermitian(1. + im*1e-15, Approx()) + @test !ishermitian(1. + im*1e-7, Approx()) + + # Wrong shape returns `false` for complex, but throws for real. + # This should be fixed somehow. + wrong_shape = rand(ComplexF64, 2, 3) + @test !issymmetric(wrong_shape, Approx()) + @test !issymmetric(wrong_shape, EachApprox()) + + @test isreal(m_symm) + @test isreal(m_symm, Approx()) + @test isreal(m_Symm) + @test isreal(m_Symm, Approx()) + @test !isdiag(m_symm) + @test !isdiag(m_Symm) + @test !isdiag(m_Symm, Approx()) +end + +@testset "istriu, istril" begin + istriu = IsApprox.istriu + istril = IsApprox.istril + for istri in (istril, istriu) + @test istri(1) + @test istri(1, Equal()) + @test istri(1, Approx()) + @test istri(1, EachApprox()) + end +end + +@testset "isnormal" begin + m = rand(ComplexF64, 3, 3); + @test !isnormal(m) + @test isnormal(m * m') + @test isnormal(m * m', Approx()) + @test isidempotent([1 0; 0 0]) + @test !isidempotent(m) +end + +@testset "commutes" begin + @test commutes(1., 2., Approx()) + @test commutes(1., 2., EachApprox()) + @test !anticommutes(1., 2., Approx()) + @test !anticommutes(1., 2., EachApprox()) + @test !anticommutes(1., 2., Equal()) + @test !anticommutes(1., 2.) + X = [0. 1.; 1. 0.]; + Xn = X + 1e-13*rand(2,2) + @test commutes(X, Xn, Approx()) + @test !commutes(X, Xn, Equal()) + @test !commutes(X, Xn) +end