Merge pull request #37 from SpM-lab/WIP_augment

Replace composite bases with augmented basis

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

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}

src/SparseIR.jl

@@ -3,15 +3,14 @@ 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
 export overlap
-export LegendreBasis, MatsubaraConstBasis
 export FiniteTempBasisSet
 export LogisticKernel, RegularizedBoseKernel
-export CompositeBasis, CompositeBasisFunction, CompositeBasisFunctionFT
+export AugmentedBasis, TauConst, TauLinear, MatsubaraConst
 export TauSampling, MatsubaraSampling, evaluate, fit, evaluate!, fit!,
        MatsubaraSampling64F, MatsubaraSampling64B, TauSampling64
 
@@ -33,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")
@@ -44,7 +40,6 @@ include("kernel.jl")
 include("sve.jl")
 include("basis.jl")
 include("augment.jl")
-include("composite.jl")
 include("sampling.jl")
 include("spr.jl")
 include("basis_set.jl")

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

src/_roots.jl

@@ -1,4 +1,4 @@
-function find_all(f, xgrid::AbstractVector)
+function find_all(f, xgrid::AbstractVector{T}) where {T}
     fx = f.(xgrid)
     hit = iszero.(fx)
     x_hit = @view xgrid[hit]
@@ -12,26 +12,33 @@ function find_all(f, xgrid::AbstractVector)
     a = @view xgrid[where_a]
     b = @view xgrid[where_b]
     fa = @view fx[where_a]
-    fb = @view fx[where_b]
 
-    x_bisect = bisect_cont.(f, a, b, fa, fb)
+    ϵ_x = if T <: AbstractFloat
+        eps(T) * maximum(abs, xgrid)
+    else
+        0
+    end
+    x_bisect = bisect.(f, a, b, fa, ϵ_x)
 
     return sort!([x_hit; x_bisect])
 end
 
-function bisect_cont(f, a, b, fa, fb)
+function bisect(f, a, b, fa, ϵ_x)
     while true
-        mid = (a + b) / 2
+        mid = midpoint(a, b)
+        closeenough(a, mid, ϵ_x) && return mid
         fmid = f(mid)
         if signbit(fa) ≠ signbit(fmid)
-            b, fb = mid, fmid
+            b = mid
         else
             a, fa = mid, fmid
         end
-        isapprox(a, b; rtol=0, atol=1e-10) && return mid
     end
 end
 
+@inline closeenough(a::T, b::T, ϵ) where {T<:AbstractFloat} = isapprox(a, b; rtol=0, atol=ϵ)
+@inline closeenough(a::T, b::T, _) where {T<:Integer} = a == b
+
 function refine_grid(grid, ::Val{α}) where {α}
     n = length(grid)
     newn = α * (n - 1) + 1
@@ -49,3 +56,59 @@ function refine_grid(grid, ::Val{α}) where {α}
     newgrid[end] = last(grid)
     return newgrid
 end
+
+function discrete_extrema(f::Function, xgrid)
+    fx::Vector{Float64} = f.(xgrid)
+    absfx               = abs.(fx)
+
+    # Forward differences: derivativesignchange[i] now means that the secant 
+    # changes sign fx[i+1]. This means that the extremum is STRICTLY between 
+    # x[i] and x[i+2].
+    gx                     = diff(fx)
+    sgx                    = signbit.(gx)
+    derivativesignchange   = @views (sgx[begin:(end - 1)] .≠ sgx[(begin + 1):end])
+    derivativesignchange_a = [derivativesignchange; false; false]
+    derivativesignchange_b = [false; false; derivativesignchange]
+
+    a      = xgrid[derivativesignchange_a]
+    b      = xgrid[derivativesignchange_b]
+    absf_a = absfx[derivativesignchange_a]
+    absf_b = absfx[derivativesignchange_b]
+    res    = bisect_discr_extremum.(f, a, b, absf_a, absf_b)
+
+    # We consider the outer point to be extremua if there is a decrease
+    # in magnitude or a sign change inwards
+    sfx = signbit.(fx)
+    if absfx[begin] > absfx[begin + 1] || sfx[begin] ≠ sfx[begin + 1]
+        pushfirst!(res, first(xgrid))
+    end
+    if absfx[end] > absfx[end - 1] || sfx[end] ≠ sfx[end - 1]
+        push!(res, last(xgrid))
+    end
+
+    return res
+end
+
+function bisect_discr_extremum(f, a, b, absf_a, absf_b)
+    d = b - a
+
+    d <= 1 && return absf_a > absf_b ? a : b
+    d == 2 && return a + 1
+
+    m      = midpoint(a, b)
+    n      = m + 1
+    absf_m = abs(f(m))
+    absf_n = abs(f(n))
+
+    if absf_m > absf_n
+        return bisect_discr_extremum(f, a, n, absf_a, absf_n)
+    else
+        return bisect_discr_extremum(f, m, b, absf_m, absf_b)
+    end
+end
+
+# This implementation of `midpoint` is performance-optimized but safe
+# only if `lo <= hi`.
+@inline midpoint(lo::T, hi::T) where {T<:Integer} = lo + ((hi - lo) >>> 0x01)
+@inline midpoint(lo::T, hi::T) where {T<:AbstractFloat} = lo + ((hi - lo) * 0.5)
+@inline midpoint(lo, hi) = midpoint(promote(lo, hi)...)

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)
 
 """
@@ -65,7 +64,7 @@ Default sampling points on the imaginary time/x axis.
 default_tau_sampling_points(::AbstractBasis) = error("unimplemented")
 
 """
-    default_matsubara_sampling_points(basis::AbstractBasis; mitigate=true)
+    default_matsubara_sampling_points(basis::AbstractBasis)
 
 Default sampling points on the imaginary frequency axis.
 """
@@ -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)
@@ -188,6 +187,8 @@ basis coefficients `G_ir[l]` to time/frequency sampled on sparse points
 """
 abstract type AbstractSampling{T,Tmat,F<:SVD} end
 
+Base.broadcastable(sampling::AbstractSampling) = Ref(sampling)
+
 function cond(sampling::AbstractSampling)
     return first(sampling.matrix_svd.S) / last(sampling.matrix_svd.S)
 end