From 184a9407041f243dc0d69f6979305b4f77927a7d Mon Sep 17 00:00:00 2001 From: araujoms Date: Fri, 7 Jun 2024 10:36:22 +0200 Subject: [PATCH] rtol is way too large, changing default to eps --- ext/KetExact/KetExact.jl | 4 +++- src/basic.jl | 10 +++++++--- src/entropy.jl | 4 ++-- src/measurements.jl | 14 +++++++------- src/nonlocal.jl | 2 +- test/basic.jl | 6 +++--- 6 files changed, 23 insertions(+), 17 deletions(-) diff --git a/ext/KetExact/KetExact.jl b/ext/KetExact/KetExact.jl index 0cd13d3..add9cb7 100644 --- a/ext/KetExact/KetExact.jl +++ b/ext/KetExact/KetExact.jl @@ -6,7 +6,9 @@ import LinearAlgebra as LA Ket._root_unity(::Type{CN.Cyc{R}}, n::Integer) where {R<:Real} = CN.E(Int64(n)) Ket._sqrt(::Type{CN.Cyc{R}}, n::Integer) where {R<:Real} = CN.root(Int64(n)) -Ket._tol(::Type{CN.Cyc{R}}) where {R<:Real} = Base.rtoldefault(R) +Ket._rtol(::Type{CN.Cyc{R}}) where {R<:Real} = Base.rtoldefault(R) +Ket._eps(::Type{CN.Cyc{R}}) where {R<:Real} = 0 +Ket._eps(::Type{CN.Cyc{R}}) where {R<:AbstractFloat} = eps(R) LA.norm(v::AbstractVector{CN.Cyc{R}}) where {R<:Real} = Ket._sqrt(CN.Cyc{R}, Int64(sum(abs2, v))) diff --git a/src/basic.jl b/src/basic.jl index 690758a..8d91c6b 100644 --- a/src/basic.jl +++ b/src/basic.jl @@ -174,20 +174,24 @@ function gell_mann!(res::AbstractMatrix{T}, i::Integer, j::Integer, d::Integer = end export gell_mann! -_tol(::Type{T}) where {T<:Number} = Base.rtoldefault(real(T)) +_rtol(::Type{T}) where {T<:Number} = Base.rtoldefault(real(T)) + +_eps(::Type{T}) where {T<:Number} = _realeps(real(T)) +_realeps(::Type{T}) where {T<:AbstractFloat} = eps(T) +_realeps(::Type{<:Real}) = 0 """ cleanup!(M::AbstractArray{T}; tol = Base.rtoldefault(real(T))) Zeroes out real or imaginary parts of `M` that are smaller than `tol`. """ -function cleanup!(M::AbstractArray{T}; tol = _tol(T)) where {T<:Number} # SD: is it type stable? +function cleanup!(M::AbstractArray{T}; tol = _eps(T)) where {T<:Number} # SD: is it type stable? wrapper = Base.typename(typeof(M)).wrapper cleanup!(parent(M); tol) return wrapper(M) end -function cleanup!(M::Array{T}; tol = _tol(T)) where {T<:Number} +function cleanup!(M::Array{T}; tol = _eps(T)) where {T<:Number} if isbitstype(T) M2 = reinterpret(real(T), M) #this is a no-op when T<:Real _cleanup!(M2; tol) diff --git a/src/entropy.jl b/src/entropy.jl index c7fa7c2..44fc303 100755 --- a/src/entropy.jl +++ b/src/entropy.jl @@ -94,8 +94,8 @@ Computes the Shannon entropy -Σᵢpᵢlog(pᵢ) of a non-negative vector `p` us Reference: [Entropy (information theory)](https://en.wikipedia.org/wiki/Entropy_(information_theory)). """ -function entropy(b::Real, p::AbstractVector) - if any(p .< -Base.rtoldefault(eltype(p))) +function entropy(b::Real, p::AbstractVector{T}) where {T<:Real} + if any(p .< -Base.rtoldefault(T)) throw(ArgumentError("p must be non-negative.")) end h = -sum(p[i] * _log(b, p[i]) for i = 1:length(p)) diff --git a/src/measurements.jl b/src/measurements.jl index 9d48275..6087aff 100644 --- a/src/measurements.jl +++ b/src/measurements.jl @@ -122,7 +122,7 @@ function test_mub(B::Vector{Matrix{T}}) where {T<:Number} sc2_exp = inv_d end sc2 = abs2(LA.dot(B[x][:, a], B[y][:, b])) - if abs2(sc2 - sc2_exp) > _tol(T) + if abs2(sc2 - sc2_exp) > _eps(T) return false end end @@ -166,16 +166,16 @@ Tests whether `vecs` is a vector of `d²` vectors |vᵢ⟩ such that |vᵢ⟩⟨ function test_sic(vecs::Vector{Vector{T}}) where {T<:Number} d = length(vecs[1]) length(vecs) == d^2 || throw(ArgumentError("Number of vectors must be d² = $(d^2), got $(length(vecs)).")) - m = zeros(T, d^2, d^2) + normalization = inv(T(d^2)) + symmetry = inv(T(d^2 * (d + 1))) for j in 1:d^2, i in 1:j - # expected scalar product + inner_product = abs2(LA.dot(vecs[i], vecs[j])) if i == j - sc2_exp = inv(T(d^2)) + deviation = abs2(inner_product - normalization) else - sc2_exp = inv(T(d^2 * (d + 1))) + deviation = abs2(inner_product - symmetry) end - sc2 = abs2(LA.dot(vecs[i], vecs[j])) - if abs2(sc2 - sc2_exp) > _tol(T) + if deviation > _eps(T) return false end end diff --git a/src/nonlocal.jl b/src/nonlocal.jl index 96e70a0..9eddc20 100755 --- a/src/nonlocal.jl +++ b/src/nonlocal.jl @@ -191,7 +191,7 @@ function correlation_tensor(p::AbstractArray{T,N2}; marg::Bool = true) where {T< res[x] = sum((-1)^sum(a[n] for n in 1:N if x[n] ≤ m[n]; init = 0) * sum(p[a, x_colon...]) for a in cia) / prod(m[n] for n in 1:N if x[n] > m[n]; init = 1) - if abs2(res[x]) < _tol(T) + if abs2(res[x]) < _eps(T) res[x] = 0 end end diff --git a/test/basic.jl b/test/basic.jl index 1215b3d..2944ecb 100644 --- a/test/basic.jl +++ b/test/basic.jl @@ -33,7 +33,7 @@ @testset "Cleanup" begin for R in [Float64, Double64, Float128, BigFloat] a = zeros(R, 2, 2) - a[1] = 0.5 * Ket._tol(R) + a[1] = 0.5 * Ket._eps(R) a[4] = 1 b = Hermitian(copy(a)) c = UpperTriangular(copy(a)) @@ -48,8 +48,8 @@ @test d == [0 0; 0 1] T = Complex{R} a = zeros(T, 2, 2) - a[1] = 0.5 * Ket._tol(T) + im - a[3] = 1 + 0.5 * Ket._tol(T) * im + a[1] = 0.5 * Ket._eps(T) + im + a[3] = 1 + 0.5 * Ket._eps(T) * im a[4] = 1 b = Hermitian(copy(a)) c = UpperTriangular(copy(a))