rtol is way too large, changing default to eps

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)))
 

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)

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))

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

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

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))