diff --git a/README.md b/README.md index ca295c1..9dd660f 100644 --- a/README.md +++ b/README.md @@ -17,11 +17,11 @@ pkg> add SparseIR ```julia using SparseIR -beta = 10.0 +β = 10.0 ωmax = 1.0 -eps = 1e-7 -basis_f = FiniteTempBasis(fermion, beta, ωmax, eps) -basis_b = FiniteTempBasis(boson, beta, ωmax, eps) +ε = 1e-7 +basis_f = FiniteTempBasis(Fermionic(), β, ωmax, ε) +basis_b = FiniteTempBasis(Bosonic(), β, ωmax, ε) ``` ## Tutorial and sample codes diff --git a/docs/src/guide.md b/docs/src/guide.md index 5467e66..3731d04 100644 --- a/docs/src/guide.md +++ b/docs/src/guide.md @@ -35,7 +35,7 @@ for us, where 80.0 is the value of the scale parameter ``\Lambda = \beta\omega_\ ### SVE -Central is the _singular value expansion_'s (SVE) computation, which is handled by the function `compute_sve`: +Central is the _singular value expansion_'s (SVE) computation, which is handled by the function `SVEResult`: Its purpose is constructing the decomposition ```math \begin{equation}\label{SVE} diff --git a/src/SparseIR.jl b/src/SparseIR.jl index 8984914..dbb9cda 100644 --- a/src/SparseIR.jl +++ b/src/SparseIR.jl @@ -3,7 +3,7 @@ Intermediate representation (IR) for many-body propagators. """ module SparseIR -export fermion, boson +export Fermionic, Bosonic export MatsubaraFreq, BosonicFreq, FermionicFreq, pioverbeta export FiniteTempBasis export SparsePoleRepresentation, to_IR, from_IR @@ -32,9 +32,6 @@ include("_specfuncs.jl") using ._LinAlg: tsvd include("freq.jl") -const boson = Bosonic() -const fermion = Fermionic() - include("abstract.jl") include("svd.jl") include("gauss.jl") diff --git a/src/_linalg.jl b/src/_linalg.jl index c471984..3a6de2b 100644 --- a/src/_linalg.jl +++ b/src/_linalg.jl @@ -27,7 +27,7 @@ function rrqr!(A::AbstractMatrix{T}; rtol=eps(T)) where {T<:AbstractFloat} taus = Vector{T}(undef, k) swapcol = Vector{T}(undef, m) - xnorms = map(norm, eachcol(A)) + xnorms = norm.(eachcol(A)) pnorms = copy(xnorms) sqrteps = sqrt(eps(T)) @@ -296,7 +296,7 @@ function svd_jacobi!(U::AbstractMatrix{T}; rtol=eps(T), maxiter=20) where {T} offd < rtol * Unorm && break end - s = map(norm, eachcol(U)) + s = norm.(eachcol(U)) U ./= transpose(s) return SVD(U, s, VT) end diff --git a/src/abstract.jl b/src/abstract.jl index 4401da1..b0768b2 100644 --- a/src/abstract.jl +++ b/src/abstract.jl @@ -18,9 +18,10 @@ as follows: where `basis.uhat[l]` is now the Fourier transform of the basis function. """ -abstract type AbstractBasis end +abstract type AbstractBasis{S <: Statistics} end Base.size(::AbstractBasis) = error("unimplemented") +Base.broadcastable(b::AbstractBasis) = Ref(b) """ Base.getindex(basis::AbstractBasis, I) @@ -30,10 +31,8 @@ Return basis functions/singular values for given index/indices. This can be used to truncate the basis to the `n` most significant singular values: `basis[1:3]`. """ -Base.getindex(::AbstractBasis, I) = error("unimplemented") - +Base.getindex(::AbstractBasis, _) = error("unimplemented") Base.firstindex(::AbstractBasis) = 1 - Base.length(basis::AbstractBasis) = length(basis.s) """ @@ -76,7 +75,7 @@ default_matsubara_sampling_points(::AbstractBasis) = error("unimplemented") Quantum statistic (Statistics instance, Fermionic() or Bosonic()). """ -statistics(basis::AbstractBasis) = basis.statistics +statistics(::AbstractBasis{S}) where {S<:Statistics} = S() """ Λ(basis::AbstractBasis) diff --git a/src/augment.jl b/src/augment.jl index df2a614..e142ed4 100644 --- a/src/augment.jl +++ b/src/augment.jl @@ -1,3 +1,19 @@ +""" + AbstractAugmentation + +Scalar function in imaginary time/frequency. + +This represents a single function in imaginary time and frequency, together +with some auxiliary methods that make it suitable for augmenting a basis. + +See also: [`AugmentedBasis`](@ref) +""" +abstract type AbstractAugmentation <: Function end + +const AugmentationTuple = Tuple{Vararg{<:AbstractAugmentation}} + +create(aug::AbstractAugmentation, ::AbstractBasis) = aug + """ AugmentedBasis <: AbstractBasis @@ -21,41 +37,42 @@ that serves as a base for multi-point functions [^shinaoka2018]. Bases augmented with `TauConst` and `TauLinear` tend to be poorly conditioned. Care must be taken while fitting and compactness should be enforced if possible to regularize the problem. - + While vertex bases, i.e. bases augmented with `MatsubaraConst`, stay reasonably well-conditioned, it is still good practice to treat the Hartree--Fock term separately rather than including it in the basis, if possible. -See also: [`MatsubaraConst`](@ref) for vertex basis [^wallerberger2021], [`TauConst`](@ref), [`TauLinear`](@ref) for multi-point [^shinaoka2018] +See also: [`MatsubaraConst`](@ref) for vertex basis [^wallerberger2021], +[`TauConst`](@ref), +[`TauLinear`](@ref) for multi-point [^shinaoka2018] [^wallerberger2021]: https://doi.org/10.1103/PhysRevResearch.3.033168 [^shinaoka2018]: https://doi.org/10.1103/PhysRevB.97.205111 """ -struct AugmentedBasis{B<:AbstractBasis} <: AbstractBasis - basis::B - augmentations::Any - u::Any - uhat::Any +struct AugmentedBasis{S<:Statistics,B<:FiniteTempBasis{S},A<:AugmentationTuple,F,FHAT} <: AbstractBasis{S} + basis :: B + augmentations :: A + u :: F + uhat :: FHAT end function AugmentedBasis(basis::AbstractBasis, augmentations...) - augmentations = Tuple(augmentation_factory(basis, augmentations...)) + augmentations = Tuple(create(aug, basis) for aug in augmentations) u = AugmentedTauFunction(basis.u, augmentations) - û = AugmentedMatsubaraFunction(basis.uhat, [n -> hat(aug, n) for aug in augmentations]) + û = AugmentedMatsubaraFunction(basis.uhat, augmentations) return AugmentedBasis(basis, augmentations, u, û) end -statistics(basis::AugmentedBasis) = statistics(basis.basis) naug(basis::AugmentedBasis) = length(basis.augmentations) -function getindex(basis::AugmentedBasis, index) - stop = range_to_size(index) +function getindex(basis::AugmentedBasis, index::AbstractRange) + stop = range_to_length(index) stop > naug(basis) || error("Cannot truncate to only augmentation.") return AugmentedBasis(basis.basis[begin:(stop - naug(basis))], basis.augmentations) end -Base.size(basis::AugmentedBasis) = (length(basis), ) +Base.size(basis::AugmentedBasis) = (length(basis),) Base.length(basis::AugmentedBasis) = naug(basis) + length(basis.basis) significance(basis::AugmentedBasis) = significance(basis.basis) accuracy(basis::AugmentedBasis) = accuracy(basis.basis) @@ -77,51 +94,50 @@ function iswellconditioned(basis::AugmentedBasis) return wbasis && waug end - +############################################################################################ +# Augmented Functions # +############################################################################################ abstract type AbstractAugmentedFunction <: Function end -struct AugmentedFunction <: AbstractAugmentedFunction - fbasis::Any - faug::Any +struct AugmentedFunction{FB,FA} <: AbstractAugmentedFunction + fbasis :: FB + faug :: FA end -augmentedfunction(a::AugmentedFunction) = a +augmentedfunction(a::AugmentedFunction) = a fbasis(a::AbstractAugmentedFunction) = augmentedfunction(a).fbasis faug(a::AbstractAugmentedFunction) = augmentedfunction(a).faug naug(a::AbstractAugmentedFunction) = length(faug(a)) Base.length(a::AbstractAugmentedFunction) = naug(a) + length(fbasis(a)) -Base.size(a::AbstractAugmentedFunction) = (length(a), ) +Base.size(a::AbstractAugmentedFunction) = (length(a),) function (a::AbstractAugmentedFunction)(x) fbasis_x = fbasis(a)(x) faug_x = [faug_l(x) for faug_l in faug(a)] - return fbasis_x .+ faug_x + return vcat(faug_x, fbasis_x) end function (a::AbstractAugmentedFunction)(x::AbstractArray) fbasis_x = fbasis(a)(x) - faug_x = (faug_l.(x) for faug_l in faug(a)) - return sum(fbasis_x .+ transpose(faug_xi) for faug_xi in faug_x) + faug_x = [faug_l.(transpose(x)) for faug_l in faug(a)] + return vcat(faug_x..., fbasis_x) end function Base.getindex(a::AbstractAugmentedFunction, r::AbstractRange) - stop = range_to_size(r) + stop = range_to_length(r) stop > naug(a) || error("Don't truncate to only augmentation") - return AugmentedFunction(fbasis(a)[begin:(stop-naug(a))], faug(a)) + return AugmentedFunction(fbasis(a)[begin:(stop - naug(a))], faug(a)) end function Base.getindex(a::AbstractAugmentedFunction, l::Integer) - if l < naug(a) - return faug(a)[l] - else - return fbasis(a)[l-naug(a)] - end + return l ≤ naug(a) ? faug(a)[l] : fbasis(a)[l - naug(a)] end - -struct AugmentedTauFunction <: AbstractAugmentedFunction - a::AugmentedFunction +### AugmentedTauFunction + +struct AugmentedTauFunction{FB,FA} <: AbstractAugmentedFunction + a::AugmentedFunction{FB,FA} end augmentedfunction(aτ::AugmentedTauFunction) = aτ.a @@ -137,19 +153,25 @@ function deriv(aτ::AugmentedTauFunction, n=1) return AugmentedTauFunction(dbasis, daug) end +### AugmentedMatsubaraFunction -struct AugmentedMatsubaraFunction <: AbstractAugmentedFunction - a::AugmentedFunction +struct AugmentedMatsubaraFunction{FB,FA} <: AbstractAugmentedFunction + a::AugmentedFunction{FB,FA} end augmentedfunction(amat::AugmentedMatsubaraFunction) = amat.a -AugmentedMatsubaraFunction(fbasis, faug) = AugmentedMatsubaraFunction(AugmentedFunction(fbasis, faug)) +function AugmentedMatsubaraFunction(fbasis, faug) + AugmentedMatsubaraFunction(AugmentedFunction(fbasis, faug)) +end zeta(amat::AugmentedMatsubaraFunction) = zeta(fbasis(amat)) +############################################################################################ +# Augmentations # +############################################################################################ -abstract type AbstractAugmentation end +### TauConst struct TauConst <: AbstractAugmentation β::Float64 @@ -159,30 +181,28 @@ struct TauConst <: AbstractAugmentation end end -function create(::Type{TauConst}, basis::AbstractBasis) - statistics(basis)::Bosonic - return TauConst(β(basis)) -end +create(::Type{TauConst}, basis::AbstractBasis{Bosonic}) = TauConst(β(basis)) function (aug::TauConst)(τ) 0 ≤ τ ≤ aug.β || throw(DomainError(τ, "τ must be in [0, β].")) return 1 / √(aug.β) end +function (aug::TauConst)(n::BosonicFreq) + iszero(n) || return 0.0 + return √(aug.β) +end +(::TauConst)(::FermionicFreq) = error("TauConst is not a Fermionic basis.") function deriv(aug::TauConst, n=1) - iszero(n) && return aug + !iszero(n) || return aug return τ -> 0.0 end -function hat(aug::TauConst, n::BosonicFreq) - iszero(n) || return 0.0 - return √(aug.β) -end - +### TauLinear struct TauLinear <: AbstractAugmentation - β::Float64 - norm::Float64 + β :: Float64 + norm :: Float64 function TauLinear(β) β > 0 || throw(DomainError(β, "Temperature must be positive.")) norm = √(3 / β) @@ -190,16 +210,19 @@ struct TauLinear <: AbstractAugmentation end end -function create(::Type{TauLinear}, basis::AbstractBasis) - statistics(basis)::Bosonic - return TauLinear(β(basis)) -end +create(::Type{TauLinear}, basis::AbstractBasis{Bosonic}) = TauLinear(β(basis)) function (aug::TauLinear)(τ) 0 ≤ τ ≤ aug.β || throw(DomainError(τ, "τ must be in [0, β].")) x = 2 / aug.β * τ - 1 return aug.norm * x end +function (aug::TauLinear)(n::BosonicFreq) + inv_w = value(n, aug.β) + inv_w = iszero(n) ? inv_w : 1 / inv_w + return aug.norm * 2 / im * inv_w +end +(::TauLinear)(::FermionicFreq) = error("TauLinear is not a Fermionic basis.") function deriv(aug::TauLinear, n=1) iszero(n) && return aug @@ -207,12 +230,7 @@ function deriv(aug::TauLinear, n=1) return τ -> 0.0 end -function hat(aug::TauLinear, n::BosonicFreq) - inv_w = value(n, aug.β) - inv_w = iszero(n) ? inv_w : 1 / inv_w - return aug.norm * 2/im * inv_w -end - +### MatsubaraConst struct MatsubaraConst <: AbstractAugmentation β::Float64 @@ -228,21 +246,8 @@ function (aug::MatsubaraConst)(τ) 0 ≤ τ ≤ aug.β || throw(DomainError(τ, "τ must be in [0, β].")) return NaN end - -deriv(aug::MatsubaraConst, n=1) = aug - -function hat(::MatsubaraConst, ::MatsubaraFreq) +function (::MatsubaraConst)(::MatsubaraFreq) return 1.0 end - -augmentation_factory(basis::AbstractBasis, augs...) = - Iterators.map(augs) do aug - if aug isa AbstractAugmentation - return aug - else - return create(aug, basis) - end - end - -create(aug, ) \ No newline at end of file +deriv(aug::MatsubaraConst, _=1) = aug diff --git a/src/basis.jl b/src/basis.jl index f980e93..69ad08e 100644 --- a/src/basis.jl +++ b/src/basis.jl @@ -4,7 +4,7 @@ Intermediate representation (IR) basis for given temperature. For a continuation kernel `K` from real frequencies, `ω ∈ [-ωmax, ωmax]`, to -imaginary time, `τ ∈ [0, beta]`, this type stores the truncated singular +imaginary time, `τ ∈ [0, β]`, this type stores the truncated singular value expansion or IR basis: K(τ, ω) ≈ sum(u[l](τ) * s[l] * v[l](ω) for l in 1:L) @@ -40,10 +40,9 @@ the variables. points `w`, you can call the function `v(w)`. To obtain a single basis function, a slice or a subset `l`, you can use `v[l]`. """ -struct FiniteTempBasis{K,T,S,TP} <: AbstractBasis +struct FiniteTempBasis{S,K,T,TP} <: AbstractBasis{S} kernel :: K - sve_result :: SVEResult{T, K, TP} - statistics :: S + sve_result :: SVEResult{T,K,TP} accuracy :: T β :: T u :: PiecewiseLegendrePolyVector{T} @@ -55,14 +54,14 @@ end """ FiniteTempBasis(statistics, β, ωmax, ε=nothing; - kernel=LogisticKernel(β * ωmax), sve_result=compute_sve(kernel; ε)) + kernel=LogisticKernel(β * ωmax), sve_result=SVEResult(kernel; ε)) Construct a finite temperature basis suitable for the given `statistics` and cutoffs `β` and `ωmax`. """ function FiniteTempBasis(statistics::Statistics, β::Number, ωmax::Number, ε=nothing; - max_size = typemax(Int), kernel=LogisticKernel(β * ωmax), - sve_result=compute_sve(kernel; ε)) + max_size=typemax(Int), kernel=LogisticKernel(β * ωmax), + sve_result=SVEResult(kernel; ε)) β > 0 || throw(DomainError(β, "Inverse temperature β must be positive")) ωmax ≥ 0 || throw(DomainError(ωmax, "Frequency cutoff ωmax must be non-negative")) @@ -75,9 +74,9 @@ function FiniteTempBasis(statistics::Statistics, β::Number, ωmax::Number, ε=n end # The polynomials are scaled to the new variables by transforming the - # knots according to: tau = beta/2 * (x + 1), w = ωmax * y. Scaling + # knots according to: tau = β/2 * (x + 1), w = ωmax * y. Scaling # the data is not necessary as the normalization is inferred. - ωmax = kernel.Λ / β + ωmax = Λ(kernel) / β u_knots = β / 2 * (u.knots .+ 1) v_knots = ωmax * v.knots u_ = PiecewiseLegendrePolyVector(u, u_knots; Δx=β / 2 * u.Δx, symm=u.symm) @@ -94,29 +93,31 @@ function FiniteTempBasis(statistics::Statistics, β::Number, ωmax::Number, ε=n û_full = PiecewiseLegendreFTVector(û_base_full, statistics; n_asymp=conv_radius(kernel)) û = û_full[1:length(s)] - return FiniteTempBasis(kernel, sve_result, statistics, accuracy, float(β), - u_, v_, s_, û, û_full) + return FiniteTempBasis(kernel, sve_result, accuracy, float(β), u_, v_, s_, û, û_full) end -const _DEFAULT_FINITE_TEMP_BASIS = FiniteTempBasis{LogisticKernel,Float64} +const DEFAULT_FINITE_TEMP_BASIS{S} = FiniteTempBasis{S,LogisticKernel,Float64} -function Base.show(io::IO, a::FiniteTempBasis{K,T}) where {K,T} - print(io, "FiniteTempBasis{$K, $T}($(statistics(a)), $(beta(a)), $(ωmax(a)))") +# TODO +function Base.show(io::IO, a::FiniteTempBasis{S,K,T}) where {S,K,T} + print(io, "FiniteTempBasis{$K, $T}($(statistics(a)), $(β(a)), $(ωmax(a)))") end -Base.getindex(basis::FiniteTempBasis, i) = - FiniteTempBasis(statistics(basis), beta(basis), ωmax(basis), nothing; - max_size=range_to_size(i), kernel=basis.kernel, sve_result=basis.sve_result) +function Base.getindex(basis::FiniteTempBasis, i) + FiniteTempBasis(statistics(basis), β(basis), ωmax(basis), nothing; + max_size=range_to_length(i), kernel=basis.kernel, + sve_result=basis.sve_result) +end significance(basis::FiniteTempBasis) = basis.s ./ first(basis.s) accuracy(basis::FiniteTempBasis) = basis.accuracy -ωmax(basis::FiniteTempBasis) = basis.kernel.Λ / beta(basis) +ωmax(basis::FiniteTempBasis) = basis.kernel.Λ / β(basis) sve_result(basis::FiniteTempBasis) = basis.sve_result kernel(basis::FiniteTempBasis) = basis.kernel -function default_tau_sampling_points(basis::FiniteTempBasis) +function default_tau_sampling_points(basis::FiniteTempBasis) x = default_sampling_points(basis.sve_result.u, length(basis)) - return β(basis)/2 * (x .+ 1) + return β(basis) / 2 * (x .+ 1) end function default_matsubara_sampling_points(basis::FiniteTempBasis) return default_matsubara_sampling_points(basis.uhat_full, length(basis)) @@ -138,50 +139,58 @@ since ``Λ == β * ωmax`` stays constant. function rescale(basis::FiniteTempBasis, new_β) new_ωmax = Λ(kernel(basis)) / new_β return FiniteTempBasis(statistics(basis), new_β, new_ωmax, nothing; - max_size=length(basis), kernel=kernel(basis), - sve_result=sve_result(basis)) + max_size=length(basis), kernel=kernel(basis), + sve_result=sve_result(basis)) end """ - finite_temp_bases(β, ωmax, ε, sve_result=compute_sve(LogisticKernel(β * ωmax); ε)) + finite_temp_bases(β, ωmax, ε, sve_result=SVEResult(LogisticKernel(β * ωmax); ε)) Construct `FiniteTempBasis` objects for fermion and bosons using the same `LogisticKernel` instance. """ function finite_temp_bases(β::AbstractFloat, ωmax::AbstractFloat, ε, - sve_result=compute_sve(LogisticKernel(β * ωmax); ε)) - basis_f = FiniteTempBasis(fermion, β, ωmax, ε; sve_result) - basis_b = FiniteTempBasis(boson, β, ωmax, ε; sve_result) + sve_result=SVEResult(LogisticKernel(β * ωmax); ε)) + basis_f = FiniteTempBasis(Fermionic(), β, ωmax, ε; sve_result) + basis_b = FiniteTempBasis(Bosonic(), β, ωmax, ε; sve_result) return basis_f, basis_b end function default_sampling_points(u::PiecewiseLegendrePolyVector, L::Integer) - (u.xmin, u.xmax) == (-1, 1) || error("expecting unscaled functions here") + (u.xmin, u.xmax) == (-1, 1) || error("Expecting unscaled functions here.") # For orthogonal polynomials (the high-T limit of IR), we know that the # ideal sampling points for a basis of size L are the roots of the L-th # polynomial. We empirically find that these stay good sampling points # for our kernels (probably because the kernels are totally positive). - L < length(u) && return roots(u[L+1]) - L < length(u) && @warn """Requesting $L sampling points but we only have - $(length(u)) basis functions in SVE.""" - - # If we do not have enough polynomials in the basis, we approximate the - # roots of the L'th polynomial by the extrema of the (L-1)'st basis - # function, which is sensible due to the strong interleaving property - # of these functions' roots. - maxima = roots(deriv(last(u))) - - # Putting the sampling points right at [0, beta], which would be the - # local extrema, is slightly worse conditioned than putting it in the - # middel. This can be understood by the fact that the roots never - # occur right at the border. - left = (first(maxima) + poly.xmin) / 2 - right = (last(maxima) + poly.xmax) / 2 - return [left; maxima; right] + if L < length(u) + x₀ = roots(u[L + 1]) + else + # If we do not have enough polynomials in the basis, we approximate the + # roots of the L'th polynomial by the extrema of the last basis + # function, which is sensible due to the strong interleaving property + # of these functions' roots. + poly = last(u) + maxima = roots(deriv(poly)) + + # Putting the sampling points right at [0, β], which would be the + # local extrema, is slightly worse conditioned than putting it in the + # middel. This can be understood by the fact that the roots never + # occur right at the border. + left = (first(maxima) + poly.xmin) / 2 + right = (last(maxima) + poly.xmax) / 2 + x₀ = [left; maxima; right] + end + + length(x₀) == L || @warn """ + Expecting to get $L sampling points for corresponding basis function, + instead got $(length(x₀)). This may happen if not enough precision is + left in the polynomial. + """ + return x₀ end -function default_matsubara_sampling_points(û::PiecewiseLegendreFTVector, L::Integer; +function default_matsubara_sampling_points(û::PiecewiseLegendreFTVector, L::Integer; fence=false) l_requested = L @@ -192,20 +201,27 @@ function default_matsubara_sampling_points(û::PiecewiseLegendreFTVector, L::In # As with the zeros, the sign changes provide excellent sampling points if l_requested < length(û) - ωn = sign_changes(û[l_requested+1]) + ωn = sign_changes(û[l_requested + 1]) else # As a fallback, use the (discrete) extrema of the corresponding # highest-order basis function in Matsubara. This turns out to be okay. - ωn = findextrema(û[L]) + ωn = find_extrema(û[L]) # For bosonic bases, we must explicitly include the zero frequency, # otherwise the condition number blows up. if statistics(û) isa Bosonic - unique!(pushfirst!(ωn, 0)) - # sort!(ωn) + pushfirst!(ωn, 0) + sort!(ωn) + unique!(ωn) end end - + + length(ωn) == l_requested || @warn """ + Expecting to get $l_requested sampling points for corresponding + $(statistics(uhat)) basis function $L, instead got $(length(ωn)). + This may happen if not enough precision is left in the polynomial. + """ + fence && fence_matsubara_sampling!(ωn) return ωn end @@ -226,7 +242,7 @@ function fence_matsubara_sampling!(ωn::Vector{<:MatsubaraFreq}) return unique!(ωn) end -function range_to_size(range::UnitRange) +function range_to_length(range::UnitRange) isone(first(range)) || error("Range must start at 1.") return last(range) end diff --git a/src/basis_set.jl b/src/basis_set.jl index 0098665..6ed466b 100644 --- a/src/basis_set.jl +++ b/src/basis_set.jl @@ -23,25 +23,25 @@ and associated sparse-sampling objects. """ struct FiniteTempBasisSet{TSF<:TauSampling,MSF<:MatsubaraSampling,TSB<:TauSampling, MSB<:MatsubaraSampling} - basis_f :: _DEFAULT_FINITE_TEMP_BASIS - basis_b :: _DEFAULT_FINITE_TEMP_BASIS + basis_f :: DEFAULT_FINITE_TEMP_BASIS{Fermionic} + basis_b :: DEFAULT_FINITE_TEMP_BASIS{Bosonic} smpl_tau_f :: TSF smpl_tau_b :: TSB smpl_wn_f :: MSF smpl_wn_b :: MSB """ - FiniteTempBasisSet(β, ωmax, ε; sve_result=compute_sve(LogisticKernel(β * ωmax); ε)) + FiniteTempBasisSet(β, ωmax, ε; sve_result=SVEResult(LogisticKernel(β * ωmax); ε)) Create basis sets for fermion and boson and associated sampling objects. Fermion and bosonic bases are constructed by SVE of the logistic kernel. """ function FiniteTempBasisSet(β::AbstractFloat, ωmax::AbstractFloat, ε; - sve_result=compute_sve(LogisticKernel(β * ωmax); ε)) + sve_result=SVEResult(LogisticKernel(β * ωmax); ε)) # Create bases using the given sve results - basis_f = FiniteTempBasis(fermion, β, ωmax, ε; sve_result) - basis_b = FiniteTempBasis(boson, β, ωmax, ε; sve_result) + basis_f = FiniteTempBasis(Fermionic(), β, ωmax, ε; sve_result) + basis_b = FiniteTempBasis(Bosonic(), β, ωmax, ε; sve_result) tau_sampling_f = TauSampling(basis_f) tau_sampling_b = TauSampling(basis_b) diff --git a/src/freq.jl b/src/freq.jl index bc45edc..5d4a262 100644 --- a/src/freq.jl +++ b/src/freq.jl @@ -51,30 +51,29 @@ imaginary-frequency axis: Ĝ(iω) = ∫ dτ exp(iωτ) G(τ) with ω = n π/β, 0 -where β is inverse temperature and by convention we include the imaginary unit -in the frequency argument, i.e, Ĝ(iω). The frequencies depend on the +where β is inverse temperature and by convention we include the imaginary unit +in the frequency argument, i.e, Ĝ(iω). The frequencies depend on the statistics of the propagator, i.e., we have that: G(τ - β) = ± G(τ) -where + is for bosons and - is for fermions. The frequencies are restricted +where + is for bosons and - is for fermions. The frequencies are restricted accordingly. - Bosonic frequency (`S == Fermionic`): `n` even (periodic in β) - Fermionic frequency (`S == Bosonic`): `n` odd (anti-periodic in β) """ struct MatsubaraFreq{S<:Statistics} <: Number - stat :: S - n :: Int + n :: Int - MatsubaraFreq(stat::Statistics, n::Integer) = new{typeof(stat)}(stat, n) + MatsubaraFreq(stat::Statistics, n::Integer) = new{typeof(stat)}(n) function MatsubaraFreq{S}(n::Integer) where {S<:Statistics} stat = S() if !allowed(stat, n) throw(ArgumentError("Frequency $(n)π/β is not $stat")) end - return new{S}(stat, n) + return new{S}(n) end end @@ -83,6 +82,8 @@ const FermionicFreq = MatsubaraFreq{Fermionic} MatsubaraFreq(n::Integer) = MatsubaraFreq(Statistics(mod(n, 2)), n) +statistics(::MatsubaraFreq{S}) where {S} = S() + """ Get prefactor `n` for the Matsubara frequency `ω = n*π/β` """ @@ -96,23 +97,23 @@ Int(a::MatsubaraFreq) = a.n """ Get value of the Matsubara frequency `ω = n*π/β` """ -value(a::MatsubaraFreq, beta::Real) = Int(a) * (π / beta) +value(a::MatsubaraFreq, β::Real) = Int(a) * (π / β) """ Get complex value of the Matsubara frequency `iω = iπ/β * n` """ -valueim(a::MatsubaraFreq, beta::Real) = 1im * value(a, beta) +valueim(a::MatsubaraFreq, β::Real) = 1im * value(a, β) """ Get statistics `ζ` for Matsubara frequency `ω = (2*m+ζ)*π/β` """ -zeta(a::MatsubaraFreq) = zeta(a.stat) +zeta(a::MatsubaraFreq) = zeta(statistics(a)) -Base.:+(a::MatsubaraFreq, b::MatsubaraFreq) = MatsubaraFreq(a.stat + b.stat, a.n + b.n) -Base.:-(a::MatsubaraFreq, b::MatsubaraFreq) = MatsubaraFreq(a.stat + b.stat, a.n - b.n) +Base.:+(a::MatsubaraFreq, b::MatsubaraFreq) = MatsubaraFreq(statistics(a) + statistics(b), a.n + b.n) +Base.:-(a::MatsubaraFreq, b::MatsubaraFreq) = MatsubaraFreq(statistics(a) + statistics(b), a.n - b.n) Base.:+(a::MatsubaraFreq) = a -Base.:-(a::MatsubaraFreq) = MatsubaraFreq(a.stat, -a.n) -Base.:*(a::BosonicFreq, c::Integer) = MatsubaraFreq(a.stat, a.n * c) +Base.:-(a::MatsubaraFreq) = MatsubaraFreq(statistics(a), -a.n) +Base.:*(a::BosonicFreq, c::Integer) = MatsubaraFreq(statistics(a), a.n * c) Base.:*(a::FermionicFreq, c::Integer) = MatsubaraFreq(a.n * c) Base.:*(c::Integer, a::MatsubaraFreq) = a * c @@ -126,14 +127,15 @@ Base.isless(a::MatsubaraFreq, b::MatsubaraFreq) = isless(a.n, b.n) # `promote_rule(<:Number, <:Number) = Number` default, together with the fact # that `@(x::Number, y::Number) = @(promote(x,y)...)` for most operations. # Let's make this error more explicit instead. -Base.promote_rule(::Type{<: MatsubaraFreq}, ::Type{<:MatsubaraFreq}) = MatsubaraFreq -Base.promote_rule(::Type{T1}, ::Type{T2}) where {T1<:MatsubaraFreq, T2<:Number} = +Base.promote_rule(::Type{<:MatsubaraFreq}, ::Type{<:MatsubaraFreq}) = MatsubaraFreq +function Base.promote_rule(::Type{T1}, ::Type{T2}) where {T1<:MatsubaraFreq,T2<:Number} throw(ArgumentError(""" Will not promote (automatically convert) $T2 and $T1. You were probably mixing a number ($T2) and a Matsubara frequency ($T1) in an additive or comparative expression, e.g., `MatsubaraFreq(0) + 1`. We disallow this. Please use `$T1(x)` explicitly.""")) +end function Base.show(io::IO, self::MatsubaraFreq) if self.n == 0 @@ -157,7 +159,8 @@ struct FreqRange{A<:Statistics} <: OrdinalRange{MatsubaraFreq{A},BosonicFreq} stop :: MatsubaraFreq{A} step :: BosonicFreq - function FreqRange(start::MatsubaraFreq{A}, stop::MatsubaraFreq{A}, step_::BosonicFreq) where {A} + function FreqRange(start::MatsubaraFreq{A}, stop::MatsubaraFreq{A}, + step_::BosonicFreq) where {A} range = Int(start):Int(step_):Int(stop) start = MatsubaraFreq{A}(first(range)) step_ = BosonicFreq(step(range)) @@ -166,9 +169,11 @@ struct FreqRange{A<:Statistics} <: OrdinalRange{MatsubaraFreq{A},BosonicFreq} end end -Base.first(self::FreqRange) = self.start -Base.last(self::FreqRange) = self.stop -Base.step(self::FreqRange) = self.step -Base.length(self::FreqRange) = Int(last(self) - first(self)) ÷ Int(step(self)) + 1 +Base.first(self::FreqRange) = self.start +Base.last(self::FreqRange) = self.stop +Base.step(self::FreqRange) = self.step +Base.length(self::FreqRange) = Int(last(self) - first(self)) ÷ Int(step(self)) + 1 Base.:(:)(start::MatsubaraFreq, stop::MatsubaraFreq) = start:BosonicFreq(2):stop -Base.:(:)(start::MatsubaraFreq, step::BosonicFreq, stop::MatsubaraFreq) = FreqRange(start, stop, step) +function Base.:(:)(start::MatsubaraFreq, step::BosonicFreq, stop::MatsubaraFreq) + FreqRange(start, stop, step) +end diff --git a/src/poly.jl b/src/poly.jl index 5341dcd..1499ca9 100644 --- a/src/poly.jl +++ b/src/poly.jl @@ -57,16 +57,11 @@ function Base.show(io::IO, p::PiecewiseLegendrePoly) return print(io, "$(typeof(p)): xmin=$(p.xmin), xmax=$(p.xmax)") end -(poly::PiecewiseLegendrePoly)(x) = _evaluate(poly, x) - -@inline function _evaluate(poly::PiecewiseLegendrePoly, x::Number) +@inline function (poly::PiecewiseLegendrePoly)(x::Number) i, x̃ = split(poly, x) return legval(x̃, view(poly.data, :, i)) * poly.norm[i] end - -@inline function _evaluate(poly::PiecewiseLegendrePoly, xs::AbstractVector) - return poly.(xs) -end +@inline (poly::PiecewiseLegendrePoly)(xs::AbstractVector) = poly.(xs) """ overlap(poly::PiecewiseLegendrePoly, f; @@ -86,7 +81,7 @@ difficulties of the integrand may occur (e.g. singularities, discontinuities). function overlap(poly::PiecewiseLegendrePoly{T}, f; rtol=eps(T), return_error=false, maxevals=10^4, points=T[]) where {T} int_result, int_error = quadgk(x -> poly(x) * f(x), - sort!(unique!([poly.knots; points]))...; + unique!(sort!([poly.knots; points]))...; rtol, order=10, maxevals) if return_error return int_result, int_error @@ -145,9 +140,9 @@ function scale(poly::PiecewiseLegendrePoly, factor) Δx=poly.Δx, symm=poly.symm) end -########################### +################################# ## PiecewiseLegendrePolyVector ## -########################### +################################# """ PiecewiseLegendrePolyVector{T} @@ -197,7 +192,7 @@ function Base.getproperty(polys::PiecewiseLegendrePolyVector, sym::Symbol) elseif sym == :data return mapreduce(poly -> poly.data, (x...) -> cat(x...; dims=3), polys) else - error("Unknown property $sym") + return getfield(polys, sym) end end @@ -258,7 +253,7 @@ end const PiecewiseLegendreFTVector{T,S} = Vector{PiecewiseLegendreFT{T,S}} -function PiecewiseLegendreFTVector(polys::PiecewiseLegendrePolyVector, +function PiecewiseLegendreFTVector(polys::PiecewiseLegendrePolyVector, stat::Statistics; n_asymp=Inf) return [PiecewiseLegendreFT(poly, stat; n_asymp) for poly in polys] end @@ -269,21 +264,21 @@ function Base.getproperty(polyFTs::PiecewiseLegendreFTVector, sym::Symbol) elseif sym == :poly return map(poly -> poly.poly, polyFTs) else - error("Unknown property $sym") + return getfield(polyFTs, sym) end end statistics(poly::PiecewiseLegendreFT) = poly.statistics statistics(polyFTs::PiecewiseLegendreFTVector) = statistics(first(polyFTs)) -# """ -# (polyFT::PiecewiseLegendreFT)(ω) +""" + (polyFT::PiecewiseLegendreFT)(ω) -# Obtain Fourier transform of polynomial for given `MatsubaraFreq` `ω`. -# """ +Obtain Fourier transform of polynomial for given `MatsubaraFreq` `ω`. +""" function (polyFT::Union{PiecewiseLegendreFT{T,S}, PiecewiseLegendreFTVector{T,S}})(ω::MatsubaraFreq{S}) where {T,S} - n = Integer(ω) + n = Int(ω) if abs(n) < polyFT.n_asymp return compute_unl_inner(polyFT.poly, n) else @@ -319,11 +314,11 @@ end Base.firstindex(::PiecewiseLegendreFT) = 1 """ - findextrema(polyFT::PiecewiseLegendreFT, part=nothing, grid=DEFAULT_GRID) + find_extrema(polyFT::PiecewiseLegendreFT; part=nothing, grid=DEFAULT_GRID) Obtain extrema of Fourier-transformed polynomial. """ -function findextrema(û::PiecewiseLegendreFT, part=nothing, grid=DEFAULT_GRID) +function find_extrema(û::PiecewiseLegendreFT; part=nothing, grid=DEFAULT_GRID) f = func_for_part(û, part) x₀ = discrete_extrema(f, grid) x₀ .= 2x₀ .+ zeta(statistics(û)) @@ -366,7 +361,7 @@ end function derivs(ppoly, x) res = [ppoly(x)] - for _ in 2:ppoly.polyorder + for _ in 2:(ppoly.polyorder) ppoly = deriv(ppoly) push!(res, ppoly(x)) end diff --git a/src/spr.jl b/src/spr.jl index 6b2106f..a99d977 100644 --- a/src/spr.jl +++ b/src/spr.jl @@ -1,26 +1,23 @@ -struct MatsubaraPoleBasis{S<:Statistics} <: AbstractBasis +struct MatsubaraPoleBasis{S<:Statistics} <: AbstractBasis{S} β :: Float64 statistics :: S poles :: Vector{Float64} end -β(basis::MatsubaraPoleBasis) = basis.β -statistics(basis::MatsubaraPoleBasis) = basis.statistics - function (basis::MatsubaraPoleBasis{S})(n::MatsubaraFreq{S}) where {S} - beta = β(basis) - iν = valueim(n, beta) + iν = valueim(n, β(basis)) if S === Fermionic return @. 1 / (iν - basis.poles) else - return @. tanh(0.5beta * basis.poles) / (iν - basis.poles) + return @. tanh(β(basis) / 2 * basis.poles) / (iν - basis.poles) end end -(basis::MatsubaraPoleBasis{S})(n::AbstractVector{MatsubaraFreq{S}}) where {S} = +function (basis::MatsubaraPoleBasis{S})(n::AbstractVector{MatsubaraFreq{S}}) where {S} mapreduce(basis, hcat, n) +end (basis::MatsubaraPoleBasis)(n::AbstractVector{<:Integer}) = basis(MatsubaraFreq.(n)) -struct TauPoleBasis{S<:Statistics} <: AbstractBasis +struct TauPoleBasis{S<:Statistics} <: AbstractBasis{S} β :: Float64 poles :: Vector{Float64} statistics :: S @@ -29,16 +26,16 @@ end ωmax(basis::TauPoleBasis) = basis.ωmax -function TauPoleBasis(beta::Real, statistics::Statistics, poles::Vector{<:AbstractFloat}) - return TauPoleBasis(beta, poles, statistics, maximum(abs, poles)) +function TauPoleBasis(β::Real, statistics::Statistics, poles::Vector{<:AbstractFloat}) + return TauPoleBasis(β, poles, statistics, maximum(abs, poles)) end function (basis::TauPoleBasis)(τ::Vector{<:AbstractFloat}) - all(τ -> 0 ≤ τ ≤ beta(basis), τ) || throw(DomainError(τ, "τ must be in [0, beta]!")) + all(τ -> 0 ≤ τ ≤ β(basis), τ) || throw(DomainError(τ, "τ must be in [0, β]!")) - x = (2 / beta(basis)) .* τ .- 1 + x = (2 / β(basis)) .* τ .- 1 y = basis.poles ./ ωmax(basis) - Λ = beta(basis) * ωmax(basis) + Λ = β(basis) * ωmax(basis) return -transpose(LogisticKernel(Λ).(x, transpose(y))) end @@ -47,27 +44,26 @@ end Sparse pole representation. """ -struct SparsePoleRepresentation{B<:AbstractBasis,T<:AbstractFloat,S<:Statistics,FMAT<:SVD} <: AbstractBasis - basis :: B - poles :: Vector{T} - u :: TauPoleBasis{S} - uhat :: MatsubaraPoleBasis{S} - statistics :: S - fitmat :: Matrix{Float64} - matrix :: FMAT +struct SparsePoleRepresentation{S<:Statistics,B<:AbstractBasis{S},T<:AbstractFloat,FMAT<:SVD + } <: AbstractBasis{S} + basis :: B + poles :: Vector{T} + u :: TauPoleBasis{S} + uhat :: MatsubaraPoleBasis{S} + fitmat :: Matrix{Float64} + matrix :: FMAT end +# TODO function Base.show(io::IO, obj::SparsePoleRepresentation) return print(io, "SparsePoleRepresentation for $(obj.basis) with poles at $(obj.poles)") end -function SparsePoleRepresentation(basis::AbstractBasis, - poles=default_omega_sampling_points(basis)) - u = TauPoleBasis(β(basis), statistics(basis), poles) - uhat = MatsubaraPoleBasis(β(basis), statistics(basis), poles) - fitmat = -basis.s .* basis.v(poles) - return SparsePoleRepresentation(basis, poles, u, uhat, - statistics(basis), fitmat, svd(fitmat)) +function SparsePoleRepresentation(b::AbstractBasis, poles=default_omega_sampling_points(b)) + u = TauPoleBasis(β(b), statistics(b), poles) + uhat = MatsubaraPoleBasis(β(b), statistics(b), poles) + fitmat = -b.s .* b.v(poles) + return SparsePoleRepresentation(b, poles, u, uhat, fitmat, svd(fitmat)) end β(obj::SparsePoleRepresentation) = β(obj.basis) @@ -80,12 +76,12 @@ function Base.getproperty(obj::SparsePoleRepresentation, d::Symbol) end end -default_tau_sampling_points(obj::SparsePoleRepresentation) = +function default_tau_sampling_points(obj::SparsePoleRepresentation) default_tau_sampling_points(obj.basis) - -default_matsubara_sampling_points(obj::SparsePoleRepresentation) = +end +function default_matsubara_sampling_points(obj::SparsePoleRepresentation) default_matsubara_sampling_points(obj.basis) - +end iswellconditioned(::SparsePoleRepresentation) = false """ diff --git a/src/sve.jl b/src/sve.jl index 80532fe..46fdbc0 100644 --- a/src/sve.jl +++ b/src/sve.jl @@ -84,8 +84,17 @@ function CentrosymmSVE(kernel, ε; InnerSVE=SamplingSVE, n_gauss, T) return CentrosymmSVE(kernel, ε, even, odd, max(even.nsvals_hint, odd.nsvals_hint)) end +struct SVEResult{T,K,TP} + u::PiecewiseLegendrePolyVector{T} + s::Vector{T} + v::PiecewiseLegendrePolyVector{T} + + kernel :: K + ε :: TP +end + """ - compute_sve(kernel::AbstractKernel; + SVEResult(kernel::AbstractKernel; Twork=nothing, ε=nothing, n_sv=typemax(Int), n_gauss=-1, T=Float64, svd_strat=:auto, sve_strat=iscentrosymmetric(kernel) ? CentrosymmSVE : SamplingSVE @@ -109,41 +118,41 @@ using a collocation). # Arguments -- `K::AbstractKernel`: Integral kernel to take SVE from. -- `ε::Real`: Accuracy target for the basis: attempt to have singular values down + - `K::AbstractKernel`: Integral kernel to take SVE from. + + - `ε::Real`: Accuracy target for the basis: attempt to have singular values down to a relative magnitude of `ε`, and have each singular value and singular vector be accurate to `ε`. A `Twork` with a machine epsilon of `ε^2` or lower is required to satisfy this. Defaults to `2.2e-16` if xprec is available, and `1.5e-8` otherwise. -- `cutoff::Real`: Relative cutoff for the singular values. A `Twork` with + - `cutoff::Real`: Relative cutoff for the singular values. A `Twork` with machine epsilon of `cutoff` is required to satisfy this. Defaults to a small multiple of the machine epsilon. - + Note that `cutoff` and `ε` serve distinct purposes. `cutoff` reprsents the accuracy to which the kernel is reproduced, whereas `ε` is the accuracy to which the singular values and vectors are guaranteed. -- `n_sv::Integer`: Maximum basis size. If given, only at most the `n_sv` most + - `n_sv::Integer`: Maximum basis size. If given, only at most the `n_sv` most significant singular values and associated singular functions are returned. -- `n_gauss (int): Order of Legendre polynomials. Defaults to kernel hinted value. -- `T`: Data type of the result. -- `Twork``: Working data type. Defaults to a data type with machine epsilon of - at most `ε^2` and at most `cutoff`, or otherwise most + - `n_gauss (int): Order of Legendre polynomials. Defaults to kernel hinted value. + - `T`: Data type of the result. + - `Twork``: Working data type. Defaults to a data type with machine epsilon of at most `ε^2`and at most`cutoff`, or otherwise most accurate data type available. -- `sve_strat::AbstractSVE`: SVE to SVD translation strategy. Defaults to `SamplingSVE`, + - `sve_strat::AbstractSVE`: SVE to SVD translation strategy. Defaults to `SamplingSVE`, optionally wrapped inside of a `CentrosymmSVE` if the kernel is centrosymmetric. -- `svd_strat` ('fast' or 'default' or 'accurate'): SVD solver. Defaults to fast - (ID/RRQR) based solution when accuracy goals are moderate, and more accurate + - `svd_strat` ('fast' or 'default' or 'accurate'): SVD solver. Defaults to fast + (ID/RRQR) based solution when accuracy goals are moderate, and more accurate Jacobi-based algorithm otherwise. Returns: An `SVEResult` containing the truncated singular value expansion. """ -function compute_sve(kernel::AbstractKernel; - Twork=nothing, cutoff=nothing, ε=nothing, n_sv=typemax(Int), - n_gauss=-1, T=Float64, svd_strat=:auto, - sve_strat=iscentrosymmetric(kernel) ? CentrosymmSVE : SamplingSVE) +function SVEResult(kernel::AbstractKernel; + Twork=nothing, cutoff=nothing, ε=nothing, n_sv=typemax(Int), + n_gauss=-1, T=Float64, svd_strat=:auto, + sve_strat=iscentrosymmetric(kernel) ? CentrosymmSVE : SamplingSVE) safe_ε, Twork, svd_strat = choose_accuracy(ε, Twork, svd_strat) sve = sve_strat(kernel, Twork(safe_ε); n_gauss, T=Twork) @@ -155,22 +164,13 @@ function compute_sve(kernel::AbstractKernel; return postprocess(sve, u, s, v, T) end -struct SVEResult{T, K, TP} - u :: PiecewiseLegendrePolyVector{T} - s :: Vector{T} - v :: PiecewiseLegendrePolyVector{T} - - kernel :: K - ε :: TP -end - function part(sve::SVEResult; ε=sve.ε, max_size=typemax(Int)) cut = count(≥(ε * first(sve.s)), sve.s) cut = min(cut, max_size) return sve.u[1:cut], sve.s[1:cut], sve.v[1:cut] end -Base.iterate(sve::SVEResult, n=1) = n ≤ 3 ? ((sve.u, sve.s, sve.v)[n], n+1) : nothing +Base.iterate(sve::SVEResult, n=1) = n ≤ 3 ? ((sve.u, sve.s, sve.v)[n], n + 1) : nothing """ matrices(sve::AbstractSVE) diff --git a/test/Project.toml b/test/Project.toml index 08af534..fc329f2 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,7 +1,6 @@ name = "SparseIR Tests" [deps] -AssociatedLegendrePolynomials = "2119f1ac-fb78-50f5-8cc0-dda848ebdb19" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5" diff --git a/test/_conftest.jl b/test/_conftest.jl index d5e0660..ad84dbd 100644 --- a/test/_conftest.jl +++ b/test/_conftest.jl @@ -1,7 +1,7 @@ if !@isdefined sve_logistic - const sve_logistic = Dict(10 => SparseIR.compute_sve(LogisticKernel(10.0)), - 42 => SparseIR.compute_sve(LogisticKernel(42.0)), - 10_000 => SparseIR.compute_sve(LogisticKernel(10_000.0)), - (10_000, 1e-12) => SparseIR.compute_sve(LogisticKernel(10_000.0); + const sve_logistic = Dict(10 => SparseIR.SVEResult(LogisticKernel(10.0)), + 42 => SparseIR.SVEResult(LogisticKernel(42.0)), + 10_000 => SparseIR.SVEResult(LogisticKernel(10_000.0)), + (10_000, 1e-12) => SparseIR.SVEResult(LogisticKernel(10_000.0); ε=1e-12)) end diff --git a/test/_util.jl b/test/_util.jl deleted file mode 100644 index aa1d99f..0000000 --- a/test/_util.jl +++ /dev/null @@ -1,15 +0,0 @@ -function typestable(@nospecialize(f), @nospecialize(t); checkonlyany=false) - v = code_typed(f, t) - stable = true - for vi in v - for (name, ty) in zip(vi[1].slotnames, vi[1].slottypes) - !(ty isa Type) && continue - if (checkonlyany && ty === Any) || - (!checkonlyany && (!Base.isdispatchelem(ty) || ty == Core.Box)) - stable = false - println("Type instability is detected! the variable is $(name) ::$ty") - end - end - end - return stable -end diff --git a/test/augment.jl b/test/augment.jl index 99fa2be..4f16cf1 100644 --- a/test/augment.jl +++ b/test/augment.jl @@ -1,13 +1,12 @@ using Test using SparseIR -using AssociatedLegendrePolynomials: Plm using LinearAlgebra @testset "augment.jl" begin @testset "Augmented bosonic basis" begin ωmax = 2 β = 1000 - basis = FiniteTempBasis(boson, β, ωmax, 1e-6) + basis = FiniteTempBasis(Bosonic(), β, ωmax, 1e-6) basis_comp = AugmentedBasis(basis, TauConst, TauLinear) # G(τ) = c - e^{-τ * pole} / (1 - e^{-β * pole}) @@ -17,23 +16,24 @@ using LinearAlgebra @test length(τ_smpl.τ) == length(basis_comp) gτ = c .+ transpose(basis.u(τ_smpl.τ)) * (-basis.s .* basis.v(pole)) magn = maximum(abs, gτ) - # @show magn # This illustrates that "naive" fitting is a problem if the fitting matrix # is not well-conditioned. gl_fit_bad = pinv(τ_smpl.matrix) * gτ gτ_reconst_bad = evaluate(τ_smpl, gl_fit_bad) @test !isapprox(gτ_reconst_bad, gτ, atol=1e-13 * magn) - # @show norm(gτ_reconst_bad - gτ) norm(gτ) 1e-13 * magn 5e-16 * cond(τ_smpl) * magn cond(τ_smpl) @test isapprox(gτ_reconst_bad, gτ, atol=5e-16 * cond(τ_smpl) * magn) + @test cond(τ_smpl) > 1e7 + @test size(τ_smpl.matrix) == (length(basis_comp), length(τ_smpl.τ)) # Now do the fit properly gl_fit = fit(τ_smpl, gτ) gτ_reconst = evaluate(τ_smpl, gl_fit) + @test isapprox(gτ_reconst, gτ, atol=1e-14 * magn) end - @testset "Vertex basis with stat = $stat" for stat in (fermion, boson) + @testset "Vertex basis with stat = $stat" for stat in (Fermionic(), Bosonic()) ωmax = 2 β = 1000 basis = FiniteTempBasis(stat, β, ωmax, 1e-6) @@ -51,4 +51,17 @@ using LinearAlgebra @test isapprox(giν_reconst, giν, atol=maximum(abs, giν) * 1e-7) end + + @testset "getindex" begin + ωmax = 2 + β = 1000 + basis = FiniteTempBasis(Bosonic(), β, ωmax, 1e-6) + basis_comp = AugmentedBasis(basis, TauConst, TauLinear) + + @test length(basis_comp.u[1:5]) == 5 + @test_throws ErrorException basis_comp.u[1:2] + @test_throws ErrorException basis_comp.u[3:7] + @test basis_comp.u[1] isa TauConst + @test basis_comp.u[2] isa TauLinear + end end diff --git a/test/composite_timings.jl b/test/composite_timings.jl index 75e3dca..07d45bd 100644 --- a/test/composite_timings.jl +++ b/test/composite_timings.jl @@ -10,11 +10,11 @@ using SparseIR #startt = time_ns() walltimes = Tuple{String,Int}[] #push!(walltimes, ("A", time_ns())) - basis = FiniteTempBasis(boson, beta, ωmax, 1e-6) + basis = FiniteTempBasis(Bosonic(), beta, ωmax, 1e-6) #push!(walltimes, ("A.1", time_ns())) - basis = FiniteTempBasis(boson, beta, ωmax, 1e-6) + basis = FiniteTempBasis(Bosonic(), beta, ωmax, 1e-6) #push!(walltimes, ("B", time_ns())) - basis_legg = LegendreBasis(boson, beta, 2) + basis_legg = LegendreBasis(Bosonic(), beta, 2) #push!(walltimes, ("C", time_ns())) basis_comp = CompositeBasis([basis_legg, basis]) #push!(walltimes, ("D", time_ns())) diff --git a/test/freq.jl b/test/freq.jl index 86ca1e8..8a56d77 100644 --- a/test/freq.jl +++ b/test/freq.jl @@ -6,10 +6,10 @@ using SparseIR @test SparseIR.zeta(MatsubaraFreq(2)) == 0 @test SparseIR.zeta(MatsubaraFreq(-5)) == 1 - @test Integer(FermionicFreq(3)) == 3 - @test Integer(BosonicFreq(-2)) == -2 + @test Int(FermionicFreq(3)) == 3 + @test Int(BosonicFreq(-2)) == -2 - @test Integer(MatsubaraFreq(Int32(4))) == 4 + @test Int(MatsubaraFreq(Int32(4))) == 4 @test_throws ArgumentError FermionicFreq(4) @test_throws ArgumentError BosonicFreq(-7) @@ -26,8 +26,8 @@ using SparseIR @test iszero(pioverbeta - pioverbeta) @test pioverbeta + oneunit(pioverbeta) == 2 * pioverbeta - @test Integer(4 * pioverbeta) == 4 - @test Integer(pioverbeta - 2 * pioverbeta) == -1 + @test Int(4 * pioverbeta) == 4 + @test Int(pioverbeta - 2 * pioverbeta) == -1 @test iszero(zero(2 * pioverbeta)) end diff --git a/test/poly.jl b/test/poly.jl index f6c1d59..f902818 100644 --- a/test/poly.jl +++ b/test/poly.jl @@ -2,15 +2,6 @@ using Test using SparseIR @testset "poly.jl" begin - @testset "typestability" begin - @test typestable(SparseIR._evaluate, - [SparseIR.PiecewiseLegendrePoly{Float64}, Float64]; - checkonlyany=true) - @test typestable(SparseIR._evaluate, - [SparseIR.PiecewiseLegendrePoly{Float64}, Vector{Float64}]; - checkonlyany=true) - end - @testset "shape" begin u, s, v = SparseIR.part(sve_logistic[42]) l = length(s) @@ -29,7 +20,7 @@ using SparseIR @testset "slice" begin sve_result = sve_logistic[42] - basis = FiniteTempBasis(fermion, 4.2, 10.0; sve_result) + basis = FiniteTempBasis(Fermionic(), 4.2, 10.0; sve_result) @test length(basis[1:4]) == 4 end @@ -44,7 +35,7 @@ using SparseIR @testset "matrix_hat" begin u, s, v = SparseIR.part(sve_logistic[42]) - uhat = SparseIR.PiecewiseLegendreFTVector(u, fermion) + uhat = SparseIR.PiecewiseLegendreFTVector(u, Fermionic()) n = MatsubaraFreq.([1, 3, 5, -1, -3, 5]) result1 = uhat[1](n) @@ -72,7 +63,7 @@ using SparseIR @testset "eval unique" begin u, s, v = SparseIR.part(sve_logistic[42]) - û = SparseIR.PiecewiseLegendreFTVector(u, fermion) + û = SparseIR.PiecewiseLegendreFTVector(u, Fermionic()) # evaluate res1 = û([1, 3, 3, 1]) diff --git a/test/runtests.jl b/test/runtests.jl index a16768c..f4383d5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,7 +4,6 @@ using Random using MultiFloats include("_conftest.jl") -include("_util.jl") @testset verbose=true "SparseIR.jl" begin include("freq.jl") @@ -20,4 +19,4 @@ include("_util.jl") include("scipost_sample_code.jl") end -nothing +nothing # otherwise we get messy output from the testset printed in the REPL diff --git a/test/sampling.jl b/test/sampling.jl index a0502cc..3d875ac 100644 --- a/test/sampling.jl +++ b/test/sampling.jl @@ -9,7 +9,7 @@ include("_conftest.jl") @testset "alias" begin β = 1 ωmax = 10 - basis = FiniteTempBasis(fermion, β, ωmax; sve_result=sve_logistic[β * ωmax]) + basis = FiniteTempBasis(Fermionic(), β, ωmax; sve_result=sve_logistic[β * ωmax]) @test TauSampling(basis) isa TauSampling64 end @@ -34,7 +34,7 @@ include("_conftest.jl") @test A \ y≈Ad \ y atol=1e-14 * norm_A rtol=0 end - @testset "fit from tau with stat = $stat, Λ = $Λ" for stat in (boson, fermion), + @testset "fit from tau with stat = $stat, Λ = $Λ" for stat in (Bosonic(), Fermionic()), Λ in (10, 42) basis = FiniteTempBasis(stat, 1, Λ; sve_result=sve_logistic[Λ]) @@ -61,7 +61,7 @@ include("_conftest.jl") end end - @testset "τ noise with stat = $stat, Λ = $Λ" for stat in (boson, fermion), Λ in (10, 42) + @testset "τ noise with stat = $stat, Λ = $Λ" for stat in (Bosonic(), Fermionic()), Λ in (10, 42) basis = FiniteTempBasis(stat, 1, Λ; sve_result=sve_logistic[Λ]) smpl = TauSampling(basis) @test issorted(smpl.sampling_points) @@ -91,12 +91,12 @@ include("_conftest.jl") @test isapprox(Gℓ, Gℓ_n, atol=12 * noise * Gℓ_magn, rtol=0) end - @testset "iω noise with stat = $stat, Λ = $Λ" for stat in (boson, fermion), + @testset "iω noise with stat = $stat, Λ = $Λ" for stat in (Bosonic(), Fermionic()), Λ in (10, 42) basis = FiniteTempBasis(stat, 1, Λ; sve_result=sve_logistic[Λ]) smpl = MatsubaraSampling(basis) - @test smpl isa (stat == fermion ? MatsubaraSampling64F : MatsubaraSampling64B) + @test smpl isa (stat == Fermionic() ? MatsubaraSampling64F : MatsubaraSampling64B) @test issorted(smpl.sampling_points) Random.seed!(1312 + 161) diff --git a/test/scipost_sample_code.jl b/test/scipost_sample_code.jl index 472b8f1..faab515 100644 --- a/test/scipost_sample_code.jl +++ b/test/scipost_sample_code.jl @@ -8,30 +8,30 @@ using FFTW β = 100 ωmax = Λ / β ϵ = 1e-8 - b = FiniteTempBasis(fermion, β, ωmax, ϵ) + b = FiniteTempBasis(Fermionic(), β, ωmax, ϵ) x = y = 0.1 - τ = 0.5β * (x + 1) + τ = β / 2 * (x + 1) ω = ωmax * y # All singular values - @show b.s - @show b.u[1](τ) - @show b.v[1](ω) + # @show b.s + # @show b.u[1](τ) + # @show b.v[1](ω) # n-th derivative of U_l(τ) and V_l(ω) for n in 1:2 u_n = SparseIR.deriv.(b.u, n) v_n = SparseIR.deriv.(b.v, n) - @show n, u_n[1](τ) - @show n, v_n[1](ω) + # @show n, u_n[1](τ) + # @show n, v_n[1](ω) end # Compute u_{ln} as a matrix for the first # 10 non-nagative fermionic Matsubara frequencies # Fermionic/bosonic frequencies are denoted by odd/even integers. hatF_t = b.uhat(1:2:19) - @show size(hatF_t) + # @show size(hatF_t) end @testset "sample 3" begin @@ -40,10 +40,10 @@ using FFTW ωmax = Λ / β ϵ = 1e-15 - @show ωmax + # @show ωmax - b = FiniteTempBasis(fermion, β, ωmax, ϵ) - @show length(b) + b = FiniteTempBasis(Fermionic(), β, ωmax, ϵ) + # @show length(b) # Sparse sampling in τ smpl_τ = TauSampling(b) diff --git a/test/spr.jl b/test/spr.jl index 6b9ba5e..b22c207 100644 --- a/test/spr.jl +++ b/test/spr.jl @@ -5,7 +5,7 @@ using Random include("_conftest.jl") @testset "spr.jl" begin - @testset "Compression with stat = $stat" for stat in (fermion, boson) + @testset "Compression with stat = $stat" for stat in (Fermionic(), Bosonic()) β = 10_000 ωmax = 1 ε = 1e-12 @@ -45,7 +45,7 @@ include("_conftest.jl") β = 2 ωmax = 21 ε = 1e-7 - basis_b = FiniteTempBasis(boson, β, ωmax, ε; sve_result=sve_logistic[β * ωmax]) + basis_b = FiniteTempBasis(Bosonic(), β, ωmax, ε; sve_result=sve_logistic[β * ωmax]) # G(iw) = sum_p coeff_p U^{SPR}(iw, omega_p) coeff = [1.1, 2.0] diff --git a/test/sve.jl b/test/sve.jl index 7042d75..99eb19d 100644 --- a/test/sve.jl +++ b/test/sve.jl @@ -23,30 +23,30 @@ a ⪅ b = leaq(a, b) @testset "sve.jl" begin @testset "smooth with Λ = $Λ" for Λ in (10, 42, 10_000) - basis = FiniteTempBasis(fermion, 1, Λ; sve_result=sve_logistic[Λ]) + basis = FiniteTempBasis(Fermionic(), 1, Λ; sve_result=sve_logistic[Λ]) check_smooth(basis.u, basis.s, 2 * maximum(basis.u(1)), 24) check_smooth(basis.v, basis.s, 50, 20) end @testset "num roots u with Λ = $Λ" for Λ in (10, 42, 10_000) - basis = FiniteTempBasis(fermion, 1, Λ; sve_result=sve_logistic[Λ]) + basis = FiniteTempBasis(Fermionic(), 1, Λ; sve_result=sve_logistic[Λ]) for ui in basis.u ui_roots = SparseIR.roots(ui) @test length(ui_roots) == ui.l end end - @testset "num roots û with stat = $stat, Λ = $Λ" for stat in (fermion, boson), + @testset "num roots û with stat = $stat, Λ = $Λ" for stat in (Fermionic(), Bosonic()), Λ in (10, 42, 10_000) basis = FiniteTempBasis(stat, 1, Λ; sve_result=sve_logistic[Λ]) for i in [1, 2, 8, 11] - x₀ = SparseIR.findextrema(basis.uhat[i]) + x₀ = SparseIR.find_extrema(basis.uhat[i]) @test i ≤ length(x₀) ≤ i + 1 end end - @testset "accuracy with stat = $stat, Λ = $Λ" for stat in (fermion, boson), + @testset "accuracy with stat = $stat, Λ = $Λ" for stat in (Fermionic(), Bosonic()), Λ in (10, 42, 10_000) basis = FiniteTempBasis(stat, 4, Λ; sve_result=sve_logistic[Λ])