add sum rule demo

Author: LittleBug

example/demoSumRule.jl

@@ -0,0 +1,36 @@
+using Lehmann
+using Printf
+using Gaston
+β = 40.0
+rtol = 1e-8
+eta = 1e-2
+N = 128
+
+dlr = DLRGrid(β = β, Euv = 1.0, isFermi = true, rtol = rtol, symmetry = :none)
+τgrid = collect(LinRange(0.0, β, N))
+println(τgrid)
+G = Sample.SemiCircle(dlr, :τ, τgrid)
+G += rand(length(G)) / 2.0 * eta
+coeff = tau2dlr(dlr, G, τgrid)
+coeff_sumrule = tau2dlr(dlr, G, τgrid, sumrule = π / 2)
+println("sum rule: ", abs(sum(coeff) - π / 2), " versus ", abs(sum(coeff_sumrule) - π / 2))
+
+
+dlr10 = DLRGrid(β = β, Euv = 100.0, isFermi = true, rtol = rtol, symmetry = :none)
+
+Gtrue = Sample.SemiCircle(dlr, :τ, dlr10.τ)
+Gfit = real.(dlr2tau(dlr, coeff, dlr10.τ))
+Gfit_sumrule = real.(dlr2tau(dlr, coeff_sumrule, dlr10.τ))
+
+@printf("%15s%30s%30s%30s\n", "tau", "true", "no sumrule", "with sumrule")
+for i in 1:dlr10.size
+    @printf("%15.6f%30.15f%30.15e%30.15e\n", dlr10.τ[i], Gtrue[i], abs(Gfit[i] - Gtrue[i]), abs(Gfit_sumrule[i] - Gtrue[i]))
+end
+
+set(term = "qt")
+# p = plot(dlr10.τ, Gtrue, label = "original", axis = "semilogy")
+xrange = (0.0, 40)
+p = plot(dlr10.τ, Gfit - Gtrue, plotstyle = :linespoints, leg = Symbol("DLR no sum rule"), Axes(xrange = xrange))
+plot!(dlr10.τ, Gfit_sumrule - Gtrue, plotstyle = :linespoints, leg = Symbol("DLR with sum rule"))
+display(p)
+readline()

example/test.jl

@@ -122,9 +122,14 @@ function _weightedLeastSqureFit(dlrGrid, Gτ, error, kernel, sumrule)
     Nτ, Nω = size(kernel)
     @assert size(Gτ)[1] == Nτ
     if isnothing(sumrule) == false
+        @assert dlrGrid.symmetry == :none && dlrGrid.isFermi "only unsymmetrized ferminoic sum rule has been implemented!"
         # println(size(Gτ))
-        kernel_m0 = kernel[:, end]
-        kernel = kernel[:, 1:Nω-1] #a copy of kernel submatrix will be created
+        M = Int(floor(dlrGrid.size / 2))
+        # M = dlrGrid.size
+
+        kernel_m0 = kernel[:, M]
+        # kernel = kernel[:, 1:Nω-1] #a copy of kernel submatrix will be created
+        kernel = hcat(kernel[:, 1:M-1], kernel[:, M+1:end])
 
         for i in 1:Nτ
             Gτ[i, :] .-= kernel_m0[i] * sumrule
@@ -161,8 +166,9 @@ function _weightedLeastSqureFit(dlrGrid, Gτ, error, kernel, sumrule)
         #add back the coeff that are fixed by the sum rule
         coeffmore = sumrule' .- sum(coeff, dims = 1)
         cnew = zeros(eltype(coeff), size(coeff)[1] + 1, size(coeff)[2])
-        cnew[1:end-1, :] = coeff
-        cnew[end, :] = coeffmore
+        cnew[1:M-1, :] = coeff[1:M-1, :]
+        cnew[M+1:end, :] = coeff[M:end, :]
+        cnew[M, :] = coeffmore
         return cnew
     else
         return coeff