improve API, fix a bug

Author: LittleBug

README.md

@@ -33,6 +33,7 @@ This package has been registered. So, simply type `import Pkg; Pkg.add("Lehmann"
 In the following [demo](example/demo.jl), we will show how to compress a Green's function of ~10000 data point into ~20 DLR coefficients, and perform fast interpolation and fourier transform up to the accuracy ~1e-10.
 
 ```julia
+using Lehmann
 β = 100.0 # inverse temperature
 Euv = 1.0 # ultraviolt energy cutoff of the Green's function
 rtol = 1e-8 # accuracy of the representation
@@ -41,9 +42,6 @@ symmetry = :none # :ph if particle-hole symmetric, :pha is antisymmetric, :none
 
 diff(a, b) = maximum(abs.(a - b)) # return the maximum deviation between a and b
 
-# Use semicircle spectral density to generate the sample Green's function
-sample(grid, type) = Sample.SemiCircle(Euv, β, isFermi, grid, type, symmetry)
-
 dlr = DLRGrid(Euv, β, rtol, isFermi, symmetry) #initialize the DLR parameters and basis
 # A set of most representative grid points are generated:
 # dlr.ω gives the real-frequency grids
@@ -53,9 +51,9 @@ dlr = DLRGrid(Euv, β, rtol, isFermi, symmetry) #initialize the DLR parameters a
 println("Prepare the Green's function sample ...")
 Nτ, Nωn = 10000, 10000 # many τ and n points are needed because Gτ is quite singular near the boundary
 τgrid = collect(LinRange(0.0, β, Nτ))  # create a τ grid
-Gτ = sample(τgrid, :τ)
+Gτ = Sample.SemiCircle(dlr, :τ, τgrid) # Use semicircle spectral density to generate the sample Green's function in τ
 ngrid = collect(-Nωn:Nωn)  # create a set of Matsubara-frequency points
-Gn = sample(ngrid, :ωn)
+Gn = Sample.SemiCircle(dlr, :n, ngrid) # Use semicircle spectral density to generate the sample Green's function in ωn
 
 println("Compress Green's function into ~20 coefficients ...")
 spectral_from_Gτ = tau2dlr(dlr, Gτ, τgrid)
@@ -64,9 +62,9 @@ spectral_from_Gω = matfreq2dlr(dlr, Gn, ngrid)
 
 println("Prepare the target Green's functions to benchmark with ...")
 τ = collect(LinRange(0.0, β, Nτ * 2))  # create a dense τ grid to interpolate
-Gτ_target = sample(τ, :τ)
+Gτ_target = Sample.SemiCircle(dlr, :τ, τ)
 n = collect(-2Nωn:2Nωn)  # create a set of Matsubara-frequency points
-Gn_target = sample(n, :ωn)
+Gn_target = Sample.SemiCircle(dlr, :n, n)
 
 println("Interpolation benchmark ...")
 Gτ_interp = dlr2tau(dlr, spectral_from_Gτ, τ)

example/demo.jl

@@ -7,9 +7,6 @@ symmetry = :none # :ph if particle-hole symmetric, :pha is antisymmetric, :none
 
 diff(a, b) = maximum(abs.(a - b)) # return the maximum deviation between a and b
 
-# Use semicircle spectral density to generate the sample Green's function
-sample(grid, type) = Sample.SemiCircle(Euv, β, isFermi, grid, type, symmetry)
-
 dlr = DLRGrid(Euv, β, rtol, isFermi, symmetry) #initialize the DLR parameters and basis
 # A set of most representative grid points are generated:
 # dlr.ω gives the real-frequency grids
@@ -19,9 +16,9 @@ dlr = DLRGrid(Euv, β, rtol, isFermi, symmetry) #initialize the DLR parameters a
 println("Prepare the Green's function sample ...")
 Nτ, Nωn = 10000, 10000 # many τ and n points are needed because Gτ is quite singular near the boundary
 τgrid = collect(LinRange(0.0, β, Nτ))  # create a τ grid
-Gτ = sample(τgrid, :τ)
+Gτ = Sample.SemiCircle(dlr, :τ, τgrid) # Use semicircle spectral density to generate the sample Green's function in τ
 ngrid = collect(-Nωn:Nωn)  # create a set of Matsubara-frequency points
-Gn = sample(ngrid, :ωn)
+Gn = Sample.SemiCircle(dlr, :n, ngrid) # Use semicircle spectral density to generate the sample Green's function in ωn
 
 println("Compress Green's function into ~20 coefficients ...")
 spectral_from_Gτ = tau2dlr(dlr, Gτ, τgrid)
@@ -30,9 +27,9 @@ spectral_from_Gω = matfreq2dlr(dlr, Gn, ngrid)
 
 println("Prepare the target Green's functions to benchmark with ...")
 τ = collect(LinRange(0.0, β, Nτ * 2))  # create a dense τ grid to interpolate
-Gτ_target = sample(τ, :τ)
+Gτ_target = Sample.SemiCircle(dlr, :τ, τ)
 n = collect(-2Nωn:2Nωn)  # create a set of Matsubara-frequency points
-Gn_target = sample(n, :ωn)
+Gn_target = Sample.SemiCircle(dlr, :n, n)
 
 println("Interpolation benchmark ...")
 Gτ_interp = dlr2tau(dlr, spectral_from_Gτ, τ)

example/demo1.jl

@@ -6,14 +6,14 @@ d = DLRGrid(Euv = 1.0, β = 1000.0, rtol = 1e-14, isFermi = false) # Initialize
 tau_k = d.τ  # DLR imaginary time points
 
 println("Prepare the Green's function sample ...")
-G_k = Sample.SemiCircle(d.Euv, d.β, d.isFermi, tau_k, :τ)  # Evaluate known G at tau_k
+G_k = Sample.SemiCircle(d, :τ, tau_k)  # Evaluate known G at tau_k
 G_x = tau2dlr(d, G_k) # DLR coeffs from G_k
 
 println("Interpolate imaginary-time Green's function ...")
 tau_i = collect(LinRange(0, d.β, 40)) # Equidistant tau grid
 G_i = dlr2tau(d, G_x, tau_i) # Evaluate DLR at tau_i
 
-G_s = Sample.SemiCircle(d.Euv, d.β, d.isFermi, tau_i, :τ)  # Evaluate known G at tau_k
+G_s = Sample.SemiCircle(d, :τ, tau_i)  # Evaluate known G at tau_k
 
 diff(a, b) = maximum(abs.(a - b)) # return the maximum deviation between a and b
 println("Maximum difference: ", diff(G_i, G_s))

example/demo2.jl

@@ -10,7 +10,7 @@ d = DLRGrid(Euv = 1.0, β = 1000.0, isFermi = false) # Initialize DLR object
 tau_i = collect(LinRange(0, d.β, 100))
 
 println("Prepare the Green's function sample ...")
-G_i = Sample.SemiCircle(d.Euv, d.β, d.isFermi, tau_i, :τ)  # Evaluate known G at tau_k
+G_i = Sample.SemiCircle(d, :τ, tau_i)  # Evaluate known G at tau_k
 # G_i = true_G_tau(tau_i, beta)
 
 G_i_noisy = G_i .+ eta .* (rand(length(G_i)) .- 0.5)

src/operation.jl

@@ -52,14 +52,16 @@ function _weightedLeastSqureFit(Gτ, error, kernel)
         C = Gτ
     else
         @assert size(error) == size(Gτ)
-        for i = 1:size(error)[1]
-            error[i, :] /= sum(error[i, :]) / length(error[i, :])
+        w = 1.0 ./ (error .+ 1e-16)
+
+        for i = 1:size(w)[1]
+            w[i, :] /= sum(w[i, :]) / length(w[i, :])
         end
 
         # W = Diagonal(weight)
         # B = transpose(kernel) * W * kernel
         # C = transpose(kernel) * W * Gτ
-        w = 1.0 ./ error
+        # w = 1.0 ./ (error .+ 1e-16)
         B = w .* kernel
         # B = Diagonal(w) * kernel
         C = w .* Gτ

src/sample/sample.jl

@@ -3,16 +3,18 @@ using FastGaussQuadrature
 using ..Spectral
 """
     SemiCircle(Euv, β, isFermi::Bool, Grid, type::Symbol, symmetry::Symbol = :none; rtol = nothing, degree = 24, regularized::Bool = true)
+    SemiCircle(dlr, type::Symbol, Grid = dlrGrid(dlr, type); degree = 24, regularized::Bool = true)
 
 Generate Green's function from a semicircle spectral density. 
 Return the function on Grid and the systematic error.
 
 #Arguments
+- `dlr`: dlrGrid struct
 - `Euv` : ultraviolet energy cutoff
 - `β` : inverse temperature
 - `isFermi`: is fermionic or bosonic
 - `Grid`: grid to evalute on
-- `type`: imaginary-time with :τ, or Matsubara-frequency with :ωn
+- `type`: imaginary-time with :τ, or Matsubara-frequency with :n
 - `symmetry`: particle-hole symmetric :ph, particle-hole antisymmetric :pha, or :none
 - `rtol`: accuracy to achieve
 - `degree`: polynomial degree for integral
@@ -23,10 +25,10 @@ function SemiCircle(Euv, β, isFermi::Bool, Grid, type::Symbol, symmetry::Symbol
     # S(ω) = sqrt(1 - (ω / Euv)^2) / Euv # semicircle -1<ω<1
     if type == :τ
         IsMatFreq = false
-    elseif type == :ωn
+    elseif type == :n
         IsMatFreq = true
     else
-        error("not implemented!")
+        error("$type is not implemented!")
     end
 
     ##### Panels endpoints for composite quadrature rule ###
@@ -48,6 +50,20 @@ function SemiCircle(Euv, β, isFermi::Bool, Grid, type::Symbol, symmetry::Symbol
     end
     return g1
 end
+function SemiCircle(dlr, type::Symbol, Grid = dlrGrid(dlr, type); degree = 24, regularized::Bool = true)
+    return SemiCircle(dlr.Euv, dlr.β, dlr.isFermi, Grid, type, dlr.symmetry; rtol = dlr.rtol, degree = degree, regularized = regularized)
+end
+
+function dlrGrid(dlr, type::Symbol)
+    if type == :τ
+        return dlr.τ
+    elseif type == :n
+        return dlr.n
+    else
+        error("$type not implemented!")
+    end
+end
+
 
 # function getG(::Val{true}, Grid)
 #     return zeros(ComplexF64, length(Grid))
@@ -100,16 +116,18 @@ end
 
 """
     MultiPole(β, isFermi::Bool, symmetry::Symbol, Grid, type::Symbol, poles, regularized::Bool = true)
+    MultiPole(dlr, type::Symbol, poles, Grid = dlrGrid(dlr, type); regularized::Bool = true)
 
 Generate Green's function from a spectral density with delta peaks specified by the argument ``poles``. 
 Return the function on Grid and the systematic error.
 
 #Arguments
+- `dlr`: dlrGrid struct
 - `β` : inverse temperature
 - `isFermi`: is fermionic or bosonic
 - `symmetry`: particle-hole symmetric :ph, particle-hole antisymmetric :pha, or :none
 - `Grid`: grid to evalute on
-- `type`: imaginary-time with :τ, or Matsubara-frequency with :ωn
+- `type`: imaginary-time with :τ, or Matsubara-frequency with :n
 - `poles`: a list of frequencies for the delta functions
 - `regularized`: use regularized bosonic kernel if symmetry = :none
 """
@@ -119,10 +137,10 @@ function MultiPole(β, isFermi::Bool, Grid, type::Symbol, poles, symmetry::Symbo
     # poles = [0.0]
     if type == :τ
         IsMatFreq = false
-    elseif type == :ωn
+    elseif type == :n
         IsMatFreq = true
     else
-        error("not implemented!")
+        error("$type is not implemented!")
     end
 
     g = IsMatFreq ? zeros(ComplexF64, length(Grid)) : zeros(Float64, length(Grid))
@@ -143,5 +161,8 @@ function MultiPole(β, isFermi::Bool, Grid, type::Symbol, poles, symmetry::Symbo
     end
     return g
 end
+function MultiPole(dlr, type::Symbol, poles, Grid = dlrGrid(dlr, type); regularized::Bool = true)
+    return MultiPole(dlr.β, dlr.isFermi, Grid, type, poles, dlr.symmetry; regularized = regularized)
+end
 
 end

test/dlr.jl

@@ -2,12 +2,14 @@ using FastGaussQuadrature, Printf
 
 rtol(x, y) = maximum(abs.(x - y)) / maximum(abs.(x))
 
-SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr.Euv, dlr.β, dlr.isFermi, grid, type, dlr.symmetry, rtol = dlr.rtol, degree = 24, regularized = true)
+# SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr.Euv, dlr.β, dlr.isFermi, grid, type, dlr.symmetry, rtol = dlr.rtol, degree = 24, regularized = true)
+SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr, type, grid, degree = 24, regularized = true)
 
 function MultiPole(dlr, grid, type)
     Euv = dlr.Euv
     poles = [-Euv, -0.2 * Euv, 0.0, 0.8 * Euv, Euv]
-    return Sample.MultiPole(dlr.β, dlr.isFermi, grid, type, poles, dlr.symmetry; regularized = true)
+    # return Sample.MultiPole(dlr.β, dlr.isFermi, grid, type, poles, dlr.symmetry; regularized = true)
+    return Sample.MultiPole(dlr, type, poles, grid; regularized = true)
 end
 
 function compare(case, a, b, eps, requiredratio, para = "")
@@ -62,9 +64,9 @@ end
         #                            Matsubara-frequency Test                                     #
         #=========================================================================================#
         # #get Matsubara-frequency Green's function
-        Gndlr = case(dlr, dlr.n, :ωn)
+        Gndlr = case(dlr, dlr.n, :n)
         nSample = dlr10.n
-        Gnsample = case(dlr, nSample, :ωn)
+        Gnsample = case(dlr, nSample, :n)
 
         # #Matsubara frequency to dlr
         coeffn = matfreq2dlr(dlr, Gndlr)
@@ -158,7 +160,7 @@ end
     dlr = DLRGrid(Euv, β, eps, isFermi, symmetry) #construct dlr basis
 
     Gτ = SemiCircle(dlr, dlr.τ, :τ)
-    Gn = SemiCircle(dlr, dlr.n, :ωn)
+    Gn = SemiCircle(dlr, dlr.n, :n)
 
     n1, n2 = 3, 4
     a = rand(n1, n2)