Begin example usage doc

Author: SamuelBadr

docs/make.jl

@@ -4,7 +4,7 @@ using Documenter
 DocMeta.setdocmeta!(SparseIR, :DocTestSetup, :(using SparseIR); recursive=true)
 
 makedocs(;
-    modules=[SparseIR],
+    # modules=[SparseIR],
     authors="Samuel Badr <samuel.badr@gmail.com> and contributors",
     repo="https://github.com/SpM-lab/SparseIR.jl/blob/{commit}{path}#{line}",
     sitename="SparseIR.jl",

docs/src/guide.md

@@ -0,0 +1,36 @@
+\newcommand{\var}{\mathrm{Var}}
+
+# Example usage and detailed explanation
+
+We will explain the inner workings of `SparseIR.jl` by means of an example use case, adapted from the `sparse-ir` paper.
+<!-- TODO: link to paper once it's released -->
+
+## Problem statement
+
+> Let us perform self-consistent second-order perturbation theory for the single impurity Anderson model at finite temperature.
+> Its Hamiltonian is given by
+> $$
+>     H = U c^\dagger_\uparrow c^\dagger_\downarrow c_\downarrow c_\uparrow + \sum_{p\sigma} \big(V_{p\sigma}  f_{p\sigma}^\dagger c_\sigma + V_{p\sigma}^* c_\sigma^\dagger c_\sigma^\dagger\big) + \sum_{p\sigma} \epsilon_{p} f_{p\sigma}^\dagger f_{p\sigma}
+> $$
+> where $U$ is the electron interaction strength, $c_\sigma$ annihilates an electron on the impurity, $f_{p\sigma}$ annihilates an electron in the bath, $\dagger$ denotes the Hermitian conjugate, $p\in\mathbb R$ is bath momentum, and $\sigma\in\{\uparrow, \downarrow\}$ is spin.
+> The hybridization strength $V_{p\sigma}$ and bath energies $\epsilon_p$ are chosen such that the non-interacting density of states is semi-elliptic with a half-bandwidth of one, $\rho_0(\omega) = \frac2\pi\sqrt{1-\omega^2}$, $U=1.2$, $\beta=10$, and the system is assumed to be half-filled.
+
+## Treatment
+
+We first import `SparseIR` and construct an appropriate basis ($\omega_\mathrm{max} = 8$ should be more than enough for this example):
+```julia-repl
+julia> using SparseIR
+
+julia> basis = FiniteTempBasis(fermion, 10, 8)
+FiniteTempBasis{LogisticKernel, Float64}(fermion, 10.0, 8.0)
+```
+There's quite a lot happening behind the scenes in this first innocuous-looking statement, so let's break it down:
+Because we did not specify otherwise, the constructor chose the analytic continuation kernel for fermions, `LogisticKernel(Λ=80.0)`, defined by
+$$
+    K(x, y) = \frac{e^{-Λ y (x + 1) / 2}}{1 + e^{-Λ y}},
+$$
+for us.
+
+Central is the _singular value expansion_'s (SVE) computation, which is handled by the function `compute_sve`:
+
+It first constructs 

docs/src/index.md

@@ -11,4 +11,5 @@ Documentation for [SparseIR](https://github.com/SpM-lab/SparseIR.jl).
 
 ```@autodocs
 Modules = [SparseIR]
+Private = false
 ```

src/_linalg.jl

@@ -11,7 +11,7 @@ Truncated rank-revealing QR decomposition with full column pivoting.
 
 Decomposes a `(m, n)` matrix `A` into the product:
 
-    A[:,piv] == Q @ R
+    A[:,piv] == Q * R
 
 where `Q` is an `(m, k)` isometric matrix, `R` is a `(k, n)` upper
 triangular matrix, `piv` is a permutation vector, and `k` is chosen such
@@ -92,9 +92,9 @@ end
 """
 Truncated singular value decomposition.
 
-Decomposes a `(m, n)` matrix `A` into the product:
+Decomposes an `(m, n)` matrix `A` into the product:
 
-    A == U @ (s[:,None] * VT)
+    A == U * (s .* VT)
 
 where `U` is a `(m, k)` matrix with orthogonal columns, `VT` is a `(k, n)`
 matrix with orthogonal rows and `s` are the singular values, a set of `k`

src/augment.jl

@@ -7,11 +7,11 @@ they used:
 
 where P_l[x] is the $l$-th Legendre polynomial.
 
-In this class, the basis functions are defined by
+In this type, the basis functions are defined by
 
     U_l(\tau) \equiv c_l (\sqrt{2l+1}/beta) * P_l[x(\tau)],
 
-where c_l are additional l-depenent constant factors.
+where c_l are additional l-dependent constant factors.
 By default, we take c_l = 1, which reduces to the original definition.
 """
 struct LegendreBasis{T<:AbstractFloat} <: AbstractBasis

src/basis.jl

@@ -10,7 +10,7 @@ getstatistics(basis::AbstractBasis) = basis.statistics
 Intermediate representation (IR) basis in reduced variables.
 
 For a continuation kernel `K` from real frequencies, `ω ∈ [-ωmax, ωmax]`, to
-imaginary time, `τ ∈ [0, β]`, this class stores the truncated singular
+imaginary time, `τ ∈ [0, β]`, this type stores the truncated singular
 value expansion or IR basis:
 
     K(x, y) ≈ sum(u[l](x) * s[l] * v[l](y) for l in range(L))
@@ -66,7 +66,6 @@ struct DimensionlessBasis{K<:AbstractKernel,T<:AbstractFloat} <: AbstractBasis
     uhat::PiecewiseLegendreFTVector{T}
     s::Vector{T}
     v::PiecewiseLegendrePolyVector{T}
-    sampling_points_v::Vector{T}
     statistics::Statistics
 end
 
@@ -91,9 +90,7 @@ function DimensionlessBasis(
     # since it has lower relative error.
     even_odd = Dict(fermion => :odd, boson => :even)[statistics]
     û = hat.(u, even_odd; n_asymp=conv_radius(kernel))
-    rts = roots(last(v))
-    sampling_points_v = [v.xmin; (rts[begin:(end - 1)] .+ rts[(begin + 1):end]) / 2; v.xmax]
-    return DimensionlessBasis(kernel, u, û, s, v, sampling_points_v, statistics)
+    return DimensionlessBasis(kernel, u, û, s, v, statistics)
 end
 
 """
@@ -114,7 +111,7 @@ end
 Intermediate representation (IR) basis for given temperature.
 
 For a continuation kernel `K` from real frequencies, `ω ∈ [-ωmax, ωmax]`, to
-imaginary time, `τ ∈ [0, beta]`, this class stores the truncated singular
+imaginary time, `τ ∈ [0, beta]`, this type stores the truncated singular
 value expansion or IR basis:
 
     K(τ, ω) ≈ sum(u[l](τ) * s[l] * v[l](ω) for l in 1:L)

src/basis_set.jl

@@ -1,9 +1,9 @@
 """
     FiniteTempBasisSet
 
-Class for holding IR bases and sparse-sampling objects.
+Type for holding IR bases and sparse-sampling objects.
 
-An object of this class holds IR bases for fermions and bosons
+An object of this type holds IR bases for fermions and bosons
 and associated sparse-sampling objects.
 
 # Fields

src/kernel.jl

@@ -111,7 +111,7 @@ on ``[0, 1] × [0, 1]`` that is given as either:
 ```math
     K_\mathrm{red}(x, y) = K(x, y) \pm K(x, -y)
 ```
-This kernel is what this class represents.  The full singular functions can
+This kernel is what this type represents.  The full singular functions can
 be reconstructed by (anti-)symmetrically continuing them to the negative
 axis.
 """

src/poly.jl

@@ -8,7 +8,6 @@ intervals ``S[i] = [a[i], a[i+1]]``, where on each interval the function
 is expanded in scaled Legendre polynomials.
 """
 struct PiecewiseLegendrePoly{T} <: Function
-    nsegments::Int
     polyorder::Int
     xmin::T
     xmax::T
@@ -24,20 +23,21 @@ struct PiecewiseLegendrePoly{T} <: Function
     norm::Vector{T}
 
     function PiecewiseLegendrePoly(
-        nsegments, polyorder, xmin, xmax, knots, Δx, data, symm, l, xm, inv_xs, norm
+        polyorder::Integer, xmin::AbstractFloat, xmax::AbstractFloat, knots::AbstractVector,
+        Δx::AbstractVector, data::AbstractMatrix, symm::Integer, l::Integer,
+        xm::AbstractVector, inv_xs::AbstractVector, norm::AbstractVector,
     )
         !any(isnan, data) || error("data contains NaN")
-        size(knots) == (nsegments + 1,) || error("Invalid knots array")
         issorted(knots) || error("knots must be monotonically increasing")
         Δx ≈ diff(knots) || error("Δx must work with knots")
         return new{eltype(knots)}(
-            nsegments, polyorder, xmin, xmax, knots, Δx, data, symm, l, xm, inv_xs, norm
+            polyorder, xmin, xmax, knots, Δx, data, symm, l, xm, inv_xs, norm
         )
     end
 end
 
 function PiecewiseLegendrePoly(data, p::PiecewiseLegendrePoly; symm=0)
-    return PiecewiseLegendrePoly(p.nsegments, p.polyorder, p.xmin, p.xmax, p.knots, p.Δx,
+    return PiecewiseLegendrePoly(p.polyorder, p.xmin, p.xmax, p.knots, p.Δx,
         data, symm, p.l, p.xm, p.inv_xs, p.norm)
 end
 
@@ -50,7 +50,7 @@ function PiecewiseLegendrePoly(
     inv_xs = 2 ./ Δx
     norm = sqrt.(inv_xs)
 
-    return PiecewiseLegendrePoly(nsegments, polyorder, first(knots), last(knots), knots,
+    return PiecewiseLegendrePoly(polyorder, first(knots), last(knots), knots,
         Δx, data, symm, l, xm, inv_xs, norm)
 end
 
@@ -236,7 +236,7 @@ struct PowerModel{T<:AbstractFloat}
     moments::Vector{T}
 end
 
-const _DEFAULT_GRID = [
+const DEFAULT_GRID = [
     range(0; length=2^6)
     trunc.(Int, exp2.(range(6, 25; length=16 * (25 - 6) + 1)))
 ]
@@ -332,11 +332,11 @@ function hat(poly::PiecewiseLegendrePoly, freq; n_asymp=Inf)
 end
 
 """
-    findextrema(polyFT::PiecewiseLegendreFT, part=nothing, grid=_DEFAULT_GRID)
+    findextrema(polyFT::PiecewiseLegendreFT, part=nothing, grid=DEFAULT_GRID)
 
 Obtain extrema of fourier-transformed polynomial.
 """
-function findextrema(polyFT::PiecewiseLegendreFT, part=nothing, grid=_DEFAULT_GRID)
+function findextrema(polyFT::PiecewiseLegendreFT, part=nothing, grid=DEFAULT_GRID)
     f = _func_for_part(polyFT, part)
     x₀ = _discrete_extrema(f, grid)
     x₀ .= 2x₀ .+ polyFT.ζ
@@ -359,7 +359,7 @@ function _func_for_part(polyFT::PiecewiseLegendreFT, part=nothing)
 end
 
 function _discrete_extrema(f::Function, xgrid)
-    fx = f.(xgrid)
+    fx = Float64.(f.(xgrid))
     absfx = abs.(fx)
 
     # Forward differences: derivativesignchange[i] now means that the secant changes sign

src/sampling.jl

@@ -1,7 +1,7 @@
 """
     AbstractSampling
 
-Abstract class for sparse sampling.
+Abstract type for sparse sampling.
 
 Encodes the "basis transformation" of a propagator from the truncated IR
 basis coefficients `G_ir[l]` to time/frequency sampled on sparse points