From 46e6393ad09e5af53f7ca58e47ea4a7e80913e73 Mon Sep 17 00:00:00 2001 From: Tao Wang Date: Tue, 4 Feb 2025 10:05:34 -0500 Subject: [PATCH] cleanup --- test/L2norm_check.jl | 270 -------------------- test/kernel_svd.jl | 557 ----------------------------------------- test/kernel_svd_sym.jl | 442 -------------------------------- test/kernel_svd_wh.jl | 528 -------------------------------------- test/svd.jl | 246 ------------------ test/svd_sym.jl | 499 ------------------------------------ test/svd_wh.jl | 332 ------------------------ 7 files changed, 2874 deletions(-) delete mode 100644 test/L2norm_check.jl delete mode 100644 test/kernel_svd.jl delete mode 100644 test/kernel_svd_sym.jl delete mode 100644 test/kernel_svd_wh.jl delete mode 100644 test/svd.jl delete mode 100644 test/svd_sym.jl delete mode 100644 test/svd_wh.jl diff --git a/test/L2norm_check.jl b/test/L2norm_check.jl deleted file mode 100644 index fe0366a..0000000 --- a/test/L2norm_check.jl +++ /dev/null @@ -1,270 +0,0 @@ -using FastGaussQuadrature, Printf -using CompositeGrids -using Lehmann -using Statistics -import Random -SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr, type, grid, degree=24, regularized=true) -#using DoubleFloats -#rtol(x, y) = maximum(abs.(x - y)) / maximum(abs.(x)) -function fine_τGrid(Λ::Float,degree,ratio::Float) where {Float} - ############## use composite grid ############################################# - # Generating a log densed composite grid with LogDensedGrid() - npo = Int(ceil(log(Λ) / log(ratio))) - 2 # subintervals on [0,1/2] in tau space (# subintervals on [0,1] is 2*npt) - grid = CompositeGrid.LogDensedGrid( - :gauss,# The top layer grid is :gauss, optimized for integration. For interpolation use :cheb - [0.0, 1.0],# The grid is defined on [0.0, β] - [0.0, 1.0],# and is densed at 0.0 and β, as given by 2nd and 3rd parameter. - npo,# N of log grid - 0.5 / ratio^(npo-1), # minimum interval length of log grid - degree, # N of bottom layer - Float - ) - return grid - -end -function fine_ωGrid(Λ::Float, degree, ratio::Float) where {Float} - ############## use composite grid ############################################# - N = Int(floor(log(Λ) / log(ratio) + 1)) - grid = CompositeGrid.LogDensedGrid( - :cheb,# The top layer grid is :gauss, optimized for integration. For interpolation use :cheb - [0.0, Λ],# The grid is defined on [0.0, β] - [0.0,],# and is densed at 0.0 and β, as given by 2nd and 3rd parameter. - N,# N of log grid - Λ / ratio^(N - 1), # minimum interval length of log grid - degree, # N of bottom layer - Float - ) - println("Composite expoential grid size: $(length(grid)), $(grid[1]), $(grid[end])") - return grid -end - -function L2normK(dlr, space, targetspace, case = nothing, pole = nothing, weight = nothing; targetdlr = nothing) - isFermi = dlr.isFermi - symmetry = dlr.symmetry - densewGrid = fine_ωGrid(dlr.Euv, 24, typeof(dlr.Euv)(1.5)) - if space ==:n - denseGrid = dlr.n #collect(-100000:100000) # LinRange(0.0, dlr.β, 100000) - kernel = Spectral.kernelΩ - else - if isnothing(targetdlr) - denseGrid = dlr.τ # LinRange(0.0, dlr.β, 100000) - else - #denseGrid = targetdlr.τ - denseGrid = dlr.τ - end - kernel = Spectral.kernelT - end - - if targetspace ==:n - fineGrid = collect(-1000:1000) - fineGrid_Linf = collect(-1000:1000) - targetkernel = Spectral.kernelΩ - else - fineGrid = fine_τGrid(dlr.Euv, 48, typeof(dlr.Euv)(1.1)) - if isnothing(targetdlr) - fineGrid_Linf = LinRange(0.0, dlr.β, 100000) - else - fineGrid_Linf = targetdlr.τ - #fineGrid_Linf = dlr.τ - end - targetkernel = Spectral.kernelT - end - #print("dlr10$(targetdlr) $(fineGrid_Linf)\n") - if isnothing(case) - value_analy = targetkernel(Float64, Val(isFermi), Val(symmetry), fineGrid, densewGrid, dlr.β) - value_analy_Linf = targetkernel(Float64, Val(isFermi), Val(symmetry), fineGrid_Linf, densewGrid, dlr.β) - value_accurate = kernel(Float64, Val(isFermi), Val(symmetry), denseGrid, densewGrid, dlr.β) - - elseif isnothing(pole) - value_analy = case(dlr, fineGrid, targetspace) - if isnothing(targetdlr) - value_analy_Linf = case(dlr, fineGrid_Linf, targetspace) - else - value_analy_Linf = case(dlr, fineGrid_Linf, targetspace) - end - value_accurate = case(dlr, denseGrid, space) - else - value_analy = case(dlr, fineGrid, targetspace, pole, weight) - value_analy_Linf = case(dlr, fineGrid_Linf, targetspace, pole, weight) - value_accurate = case(dlr, denseGrid, space, pole, weight) - end - #value_Linf2 = 0 - if space == :τ && targetspace == :τ - value = tau2tau(dlr, value_accurate, fineGrid, denseGrid, axis =1) - - value_Linf = tau2tau(dlr, value_accurate, fineGrid_Linf, denseGrid, axis =1) - #value_Linf2 = tau2tau(dlr, value_analy_Linf, denseGrid, fineGrid_Linf ) - print("Ginterp $(value_Linf[1:10])\n") - elseif space == :τ && targetspace == :n - value = tau2matfreq(dlr, value_accurate, fineGrid, denseGrid, axis =1) - value_Linf = tau2matfreq(dlr, value_accurate, fineGrid_Linf, denseGrid, axis =1) - elseif space == :n && targetspace == :τ - value = matfreq2tau(dlr, value_accurate, fineGrid, denseGrid, axis =1) - value_Linf = matfreq2tau(dlr, value_accurate, fineGrid_Linf, denseGrid, axis =1) - elseif space == :n && targetspace == :n - value = matfreq2matfreq(dlr, value_accurate, fineGrid, denseGrid, axis =1) - value_Linf = matfreq2matfreq(dlr, value_accurate, fineGrid_Linf, denseGrid, axis =1) - end - #print("value_analy $(value_analy[1:10])\n" ) - if targetspace == :n - L2_error = maximum(sum( abs.(value-value_analy).^2, dims =1))/dlr.β - else - L2_error = maximum(Interp.integrate1D( real.(value-value_analy).^2, fineGrid)) - end - Linf_error = maximum(abs.(value_Linf-value_analy_Linf).^2) - return L2_error, Linf_error - #print("$(interp_analy) $(interp_analy)\n") -end - - - -function L2normτ(value_dlr, dlr, case, poles=nothing, weight = nothing) - - fineGrid = fine_τGrid(dlr.Euv, 12, typeof(dlr.Euv)(1.5)) - value = real(tau2tau(dlr, value_dlr, fineGrid, dlr.τ )) - if case == MultiPole - value_analy = case(dlr, fineGrid, :τ, poles, weight) - else - value_analy = case(dlr, fineGrid, :τ) - end - #print("value_analy $(value_analy[1:10])\n" ) - interp = Interp.integrate1D( value , fineGrid) - interp_analy = Interp.integrate1D( value_analy , fineGrid) - #print("$(interp_analy) $(interp_analy)\n") - return abs(interp_analy), abs(interp-interp_analy), maximum(abs.(value - value_analy))/ maximum(abs.(value_analy)) -end -#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, type, poles, grid; regularized=true) -#end - -function MultiPole(dlr, grid, type, coeff, weight = nothing) - Euv = dlr.Euv - poles = coeff * Euv - # return Sample.MultiPole(dlr.β, dlr.isFermi, grid, type, poles, dlr.symmetry; regularized = true) - if isnothing(weight) - return Sample.MultiPole(dlr, type, poles, grid; regularized=true) - else - return Sample.MultiPole(dlr, type, poles, grid, weight; regularized=true) - end -end - -function test_dlr_coeff(case, isFermi, symmetry, Euv, β, eps, eff_poles, weight; dtype=Float64, output = false) - para = "fermi=$isFermi, sym=$symmetry, Euv=$Euv, β=$β, rtol=$eps" - dlr = DLRGrid(Euv, β, eps, isFermi, symmetry, dtype=dtype) #construct dlr basis - #dlr10 = DLRGrid(10Euv, β, eps, isFermi, symmetry) #construct denser dlr basis for benchmark purpose - dlr10 = DLRGrid(Euv, β, eps, isFermi, symmetry, dtype=dtype) #construct denser dlr basis for benchmark purpose - - N_poles = 1000 - N = 20 - - Gndlr = case(dlr, dlr.n, :n, eff_poles, weight) - τSample = dlr10.τ - Gsample = case(dlr, τSample, :τ, eff_poles, weight) - - Gfourier = matfreq2tau(dlr, Gndlr, τSample) - dlreff = matfreq2dlr(dlr,Gndlr) - dlreff = imag(dlreff) - print("$(symmetry) $(Euv) $(eps) max $(maximum(abs.(dlreff) )) min $(minimum(abs.(dlreff )))\n") -end - - -function test_errK(isFermi, symmetry, Euv, β, eps, space, targetspace; dtype=Float64, output = false) - # println("Test $case with isFermi=$isFermi, Symmetry = $symmetry, Euv=$Euv, β=$β, rtol=$eps") - #N_poles = 100 - para = "fermi=$isFermi, sym=$symmetry, Euv=$Euv, β=$β, rtol=$eps" - dlr = DLRGrid(Euv, β, eps, isFermi, symmetry, dtype=dtype) #construct dlr basis - L2_error, Linf_error = L2normK(dlr,space, targetspace) - print("K error: L2: $(sqrt(L2_error)), Linf: $(sqrt(Linf_error)), eps: $(eps)\n") -end - - -function test_err(case, isFermi, symmetry, Euv, β, eps, poles, weights, space, targetspace; dtype=Float64, output = false) - # println("Test $case with isFermi=$isFermi, Symmetry = $symmetry, Euv=$Euv, β=$β, rtol=$eps") - #N_poles = 100 - para = "fermi=$isFermi, sym=$symmetry, Euv=$Euv, β=$β, rtol=$eps" - dlr = DLRGrid(Euv, β, eps, isFermi, symmetry, dtype=dtype) #construct dlr basis - dlr10 = DLRGrid(10Euv, β, eps, isFermi, symmetry, dtype=dtype) - #dlr10 = DLRGrid(10Euv, β, eps, isFermi, symmetry) #construct denser dlr basis for benchmark purpose - N_poles = size(poles)[2] - N = size(poles)[1] - L2_sum = 0 - Linf_sum = 0 - for i in 1:N - eff_poles = poles[i,:] - weight = weights[i,:] - L2_error, Linf_error = L2normK(dlr, space, targetspace, case, eff_poles, weights, targetdlr = dlr10 ) - print("G error: L2: $(sqrt(L2_error)), Linf:$(sqrt(Linf_error)), eps: $(eps)\n") - L2_sum += L2_error - Linf_sum += Linf_error - end - - if output - file = open("./accuracy_test3.dat", "a") - #@printf(file, "%48.40g %48.40g %48.40g\n", eps, abs(b-c), ) - #@printf(file, "%24.20g %24.20g %24.20g %24.20g %24.20g %24.20g\n", eps, value_sum/N, err_sum/N, err_sum /N/eps*Euv, max_err_sum/N, eta/N) - @printf(file, "%24.20g %24.20g %24.20g\n", eps, sqrt(L2_sum/N)/eps, sqrt(Linf_sum/N)/eps ) - close(file) - end -end - -function test_err(case, isFermi, symmetry, Euv, β, eps, space, targetspace; dtype=Float64, output = false) - # println("Test $case with isFermi=$isFermi, Symmetry = $symmetry, Euv=$Euv, β=$β, rtol=$eps") - #N_poles = 100 - para = "fermi=$isFermi, sym=$symmetry, Euv=$Euv, β=$β, rtol=$eps" - dlr = DLRGrid(Euv, β, eps, isFermi, symmetry, dtype=dtype) #construct dlr basis - dlr10 = DLRGrid(10Euv, β, eps, isFermi, symmetry, dtype=dtype) - #dlr10 = DLRGrid(10Euv, β, eps, isFermi, symmetry) #construct denser dlr basis for benchmark purpose - L2_error, Linf_error = L2normK(dlr, space, targetspace, case, targetdlr = dlr10) - print("G error: L2: $(sqrt(L2_error)), Linf:$(sqrt(Linf_error)), eps: $(eps)\n") - - if output - file = open("./accuracy_test2.dat", "a") - #@printf(file, "%48.40g %48.40g %48.40g\n", eps, abs(b-c), ) - @printf(file, "%24.20g %24.20g %24.20g\n", eps, sqrt(L2_error)/eps, sqrt(Linf_error)/eps ) - close(file) - end -end - - -Λ = [1e4,1e5, 1e6] -#rtol = [ 1e-12] -rtol = [1e-4, 1e-6, 1e-8, 1e-10, 1e-12] -case = MultiPole -N_poles = 1000 -N = 10 -#setprecision(128) -dtype = Float64 -#dtype = BigFloat -poles = zeros(dtype,(N,N_poles) ) -weights = zeros(dtype, (N,N_poles)) -Random.seed!(1234) -for i in 1:N - poles[i,:] = 2.0*rand(dtype, N_poles) .- 1.0 - weights[i,:] = 2.0 *rand(dtype, N_poles) .- 1.0 -end -weights = weights/abs(sum(weights)) -for l in Λ - for r in rtol - #test_errK(true, :none, l, 1.0, r) - space = :τ - #targetspace = :τ - #space = :n - targetspace = :n - #test_err(case, true, :none, l, 1.0, r, poles , weights, dtype = Float64, output = true) - test_errK(true, :none, l, 1.0, r,space, targetspace) - if case == MultiPole - test_err(case, true, :none, l, 1.0, r, poles , weights, space,targetspace, dtype = Float64, output = true) - #test_err(case, false, :none, l, 1.0, r, poles , weights, space,targetspace, dtype = Float64, output = true) - else - test_err(case, true, :none, l, 1.0, r, space,targetspace, dtype = Float64, output = true) - #test_err(case, false, :none, l, 1.0, r, space,targetspace, dtype = Float64, output = true) - end - #test_err(case, true, :sym, l, 1.0, r, poles, weights, dtype = Float64, output = true) - end -end - - - diff --git a/test/kernel_svd.jl b/test/kernel_svd.jl deleted file mode 100644 index ac05a55..0000000 --- a/test/kernel_svd.jl +++ /dev/null @@ -1,557 +0,0 @@ -using Lehmann -using Printf -using CompositeGrids -using LinearAlgebra -using GenericLinearAlgebra -using Random -using Plots -using LaTeXStrings -include("grid.jl") -""" -composite expoential grid -""" -function Htau(tau::Vector{T}, weight::Vector{T}, gamma) where {T} - result = zeros(T, (length(tau), length(tau))) - for i in eachindex(tau) - for j in eachindex(tau) - result[i, j] = sqrt(weight[j] * weight[i]) * (exp(-(tau[i] + tau[j])) - exp(-(tau[i] + tau[j]) * gamma)) / (tau[i] + tau[j]) - end - end - #print(result) - return result -end - - -@inline function F1(a::T, b::T) where {T} - if abs(a + b) > EPS - return (1 - exp(-(a + b))) / (a + b) - else - return T(1-(a+b)/2 + (a+b)^2/6 - (a+b)^3/24) - end -end - -""" -``G(x, y) = (exp(-x)-exp(-y))/(x-y)`` -``G(x, x) = -exp(-x)`` -""" -@inline function G1(a::T, b::T) where {T} - if abs(a - b) > EPS - return (exp(-a) - exp(-b)) / (b - a) - else - return (exp(-a) + exp(-b)) / 2 - end -end - -function Homega(omega::Vector{T}, weight::Vector{T}) where {T} - result = zeros(T, (length(omega), length(omega))) - for i in eachindex(omega) - for j in eachindex(omega) - if omega[i]*omega[j]>0 - #result[i, j] = sqrt(weight[j] * weight[i]) * F1(abs(omega[i]), abs(omega[j])) - result[i, j] = F1(abs(omega[i]), abs(omega[j])) - else - result[i, j] = G1(abs(omega[i]), abs(omega[j])) - end - end - end - #print(result) - return result -end - -function IntFermiT(omega::T) where {T} - omega_new = omega/2 - if omega_new < 1e-6 - return 0.5*(1 - omega_new^2/3 + 2omega_new^4/15 - 17omega_new^6/315) - else - return tanh(omega_new)/omega - end -end - - -function Kfunc(omega::Vector{T}, tau::Vector{T}, weight_omega::Vector{T}, weight_tau::Vector{T} , ifregular::Bool=false, omega0::T = 1e-4) where {T} - result = zeros(T, (length(tau), length(omega))) - omega0 = T(omega0) - for i in eachindex(tau) - for j in eachindex(omega) - #result[i, j] = sqrt(weight_tau[i] * weight_omega[j]) * Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - #result[i, j] = sqrt(weight_tau[i] * weight_omega[j]) * Spectral.kernelFermiT(tau[i], omega[j], T(1.0)) - if ifregular - #result[i, j] = Spectral.kernelFermiT(tau[i], omega[j], T(1.0)) - 1/2.0 # - IntFermiT(omega[j]) - result[i, j] =sqrt(weight_tau[i])*(Spectral.kernelFermiT(tau[i], omega[j], T(1.0)) - 1.0/2.0 - (1- 2*tau[i])*omega[j]/4.0 - (tau[i]-1)*tau[i]*omega[j]^2/4.0) - - #(Spectral.kernelFermiT(tau[i], omega0, T(1.0)) - #+ Spectral.kernelFermiT(1.0 - tau[i], omega0, T(1.0)))/2.0 - else - #result[i, j] = Spectral.kernelFermiT(tau[i], omega[j], T(1.0)) - result[i, j] = sqrt(weight_tau[i])*Spectral.kernelFermiT(tau[i], omega[j], T(1.0)) - #result[i, j] = Spectral.kernelFermiT(tau[i], omega[j], T(1.0)) - end - #result[i, j] = Spectral.kernelFermiT_PH(tau[i], omega[j], T(1.0)) - #result[i, j] = weight_tau[i] * weight_omega[j]*kernelSymT_test(tau[i], omega[j], T(1.0)) - #result[i, j] = sqrt(weight_omega[i] * weight_tau[j]) * - end - end - #print(result) - return result -end - -function Kfunc(omega::Vector{T}, tau::Vector{T}, ifregular::Bool=false, omega0::T = 1e-4) where {T} - result = zeros(T, (length(tau), length(omega))) - omega0 = T(omega0) - for i in eachindex(tau) - for j in eachindex(omega) - if ifregular - - result[i, j] =Spectral.kernelFermiT(tau[i], omega[j], T(1.0)) - 1.0/2.0 - (1- 2*tau[i])*omega[j]/4.0 - (tau[i]-1)*tau[i]*omega[j]^2/4.0 - - else - result[i, j] = Spectral.kernelFermiT(tau[i], omega[j], T(1.0)) - - end - - end - end - - return result -end - - -function Kfunc_freq(omega::Vector{T}, n::Vector{Int}, weight_omega::Vector{T}, regular::Bool=false ,omega0::T=1e-4 ) where {T} - result = zeros(Complex{T}, (length(n), length(omega))) - omega0 = T(omega0) - for i in eachindex(n) - omega_n = (2*n[i]+1)* π - for j in eachindex(omega) - if regular - result[i, j] =(Spectral.kernelFermiΩ(n[i], omega[j], T(1.0)) + 1/(im*omega_n) + omega[j]/(im*omega_n)^2 + omega[j]^2/(im*omega_n)^3) - else - result[i, j] =Spectral.kernelFermiΩ(n[i], omega[j], T(1.0)) - end - end - end - return result -end - - -function Kfunc_freq(omega::Vector{T}, n::Vector{Int}, regular::Bool=false ,omega0::T=1e-4 ) where {T} - result = zeros(Complex{T}, (length(n), length(omega))) - omega0 = T(omega0) - for i in eachindex(n) - omega_n = (2*n[i]+1)* π - for j in eachindex(omega) - if regular - result[i, j] =(Spectral.kernelFermiΩ(n[i], omega[j], T(1.0)) + 1/(im*omega_n) + omega[j]/(im*omega_n)^2 + omega[j]^2/(im*omega_n)^3) - else - result[i, j] =Spectral.kernelFermiΩ(n[i], omega[j], T(1.0)) - end - end - end - return result -end - -function Kfunc_freq!(result::Matrix{Complex{T}}, omega::Vector{T}, n::Vector{Int}, regular::Bool=false ,omega0::T=1e-4 ) where {T} - omega0 = T(omega0) - for i in eachindex(n) - omega_n = (2*n[i]+1)* π - for j in eachindex(omega) - if regular - result[i, j] =(Spectral.kernelFermiΩ(n[i], omega[j], T(1.0)) + 1/(im*omega_n) + omega[j]/(im*omega_n)^2 + omega[j]^2/(im*omega_n)^3) - else - result[i, j] =Spectral.kernelFermiΩ(n[i], omega[j], T(1.0)) - end - end - end -end - -function Kfunc_expan(omega::Vector{T}, n::Vector{Int}, weight_omega::Vector{T}, Lambda::T) where {T} - result = zeros(T, (length(n), length(omega))) - for i in eachindex(n) - for j in eachindex(omega) - result[i, j] = (omega[j]/Lambda)^(i-1)*weight_omega[j] - end - end - return result -end - - -function IR(grid, U, idx, name, Lambda=100, ifplot = false) - Un = U[:, 1:idx] - qr_Un = qr(transpose(Un), Val(true)) - qr_nidx = qr_Un.p - left = searchsortedfirst(grid, -Lambda) - right = searchsortedfirst(grid, Lambda) - #print("selected: $(length(qr_nidx)) $(minimum(qr_nidx[1:idx])) $(maximum(qr_nidx[1:idx]))\n" ) - if name == "omega_n" - shift = 0#20 - idx - else - shift = 0 - end - ir_grid = sort(grid[qr_nidx[1:idx+shift]]) - Unew = Un[sort(qr_nidx[1:idx+shift]),:] - Hnew = Unew*transpose(conj(Unew)) - Hfull = Un*transpose(conj(Un)) - #print("eigenvalue: $(eigvals(Hfull))\n") - if ifplot - #for plidx in [1,5,10, 15, 20, 25, 30] - #plot!(pl, grid[qr_nidx[1:idx]] , abs.(Un[qr_nidx[1:idx],plidx]), seriestype=:scatter, markershape=:circle) - - #end - pl = plot() - if name == "omega_n" - heatmap!(pl, abs.(Unew), title=L"U_{IR}", xlabel=L"s", ylabel=L"\omega_n", color=:viridis) - else - heatmap!(pl, abs.(Unew), title=L"U_{IR}", xlabel=L"s", ylabel=L"\tau", color=:viridis) - end - #xlabel!(L"\omega") - #legend() - #pl = plot(wgrid , abs.(Kn[1,:]) ) - savefig(pl, name*"UIR.pdf") - if name == "omega_n" - pl = heatmap(abs.(Un[left:right, :]), title=L"U_{full}", xlabel=L"s", ylabel=L"\tau", color=:viridis) - else - pl = heatmap(abs.(Un), title=L"U_{full}", xlabel=L"s", ylabel=L"\tau", color=:viridis) - end - diag_full = diag(abs.(Hfull)) - # print("max Ufull:$(maximum(diag_full))\n") - # print("max Ufull:$(minimum(diag_full))\n") - diag_IR = diag(abs.(Hnew)) - # print("max UIR:$(maximum(diag_IR))\n") - # print("max UIR:$(minimum(diag_IR))\n") - savefig(pl, name*"Ufull.pdf") - end - #print(name*" R cond :$(cond(qr_Un.R[:, qr_nidx[1:idx]]))\n") - #print(name*" Q cond :$(cond(qr_Un.Q))\n") - #print(name*" IR cond :$(cond(Un[qr_nidx[1:idx] , 1:idx]))\n") - #print(name*" Full cond:$(cond(Un[:, 1:idx]))\n") - #print(name*" IR cond :$(cond(Unew))\n") - #print(name*" Full cond:$(cond(Un))\n") - # print(name*" IR cond UUT:$(cond(abs.(Hnew)))\n") - # print(name*" Full cond UUT:$(cond(abs.(Hfull)))\n") - return qr_nidx, ir_grid -end - -# function generate_grid(eps::T, Lambda::T, n_trunc::T, space::Symbol=:τ, regular = false, -# omega0::T = Lambda,) where {T} -# # generate frequency finegrid -# w_grid = fine_ωGrid_test(T(Lambda), 24, T(1.5)) -# weight_w = zeros(T,length(w_grid)) -# #calculate grid weights -# for i in 1:length(w_grid) -# data = zeros(T,length(w_grid)) -# data[i] = 1.0 -# weight_w[i] = Interp.integrate1D(data, w_grid) -# end - -# #symmetrize the grid -# wgrid = vcat(-w_grid.grid[end:-1:1], w_grid.grid) -# weight_w = vcat(weight_w[end:-1:1], weight_w) - - -# #generate tau fine grid -# t_grid = fine_τGrid_test(T(Lambda),128, T(1.5)) -# weight_t = zeros(T,length(t_grid)) -# for i in 1:length(t_grid) -# data = zeros(T,length(t_grid)) -# data[i] = 1.0 -# weight_t[i] = Interp.integrate1D(data, t_grid) -# end -# tgrid = t_grid.grid - -# # generate fine n grid - -# #ngrid = nGrid_test(true, T(Lambda), 12, T(1.5)) -# ngrid = uni_ngrid(true, T(n_trunc*Lambda)) -# omega = (2*ngrid.+1)* π - -# #ngrid = vcat(ngrid, dlr.n) -# #unique!(ngrid) -# #ngrid = sort(ngrid) - -# #regular controls if we add 1/(iwn - Lambda) term to compensate the tail -# expangrid = collect(1:10000) - -# F = Fourier(ngrid, tgrid, weight_t) -# Kn = Kfunc_freq(wgrid, Int.(ngrid), weight_w, regular, omega0) -# Ktau = Kfunc(wgrid, tgrid, weight_w, weight_t, regular, omega0) -# Kexpan = Kfunc_expan(wgrid, expangrid, weight_w, Lambda) -# # Kn_fourier = F*Ktau -# # print("fourier error: $(maximum(abs.(Kn- Kn_fourier)))\n") -# left = searchsortedfirst(omega, -n_trunc*Lambda/10) -# right = searchsortedfirst(omega, n_trunc*Lambda/10) -# # print("$(left) $(right)\n") -# # Kn_new = copy(Kn[left:right, :]) -# # Kn_new[1, :] = sum(Kn[1:left, :], dims=1) -# # Kn_new[end, :] = sum(Kn[right:end, :], dims=1) - -# # #print("$(maximum(abs.(Kn), dims=1))\n $(abs.(Kn[1,:]))\n $(abs.(Kn[end,:])) \n") - -# # Kn = Kn_new - -# if space == :n -# eig = svd(Kn, full = true) -# elseif space == :τ -# eig = svd(Ktau, full = true) -# elseif space == :ex -# eig = svd(Kexpan, full = true) -# print("eig highfreq expansion: $(eig.S[1:10])\n") -# pl = plot( collect(1:70), eig.S[1:70], linewidth = 1,label = L"max(\left|K_n\right|)", yaxis=:log) -# # plot!(pl,wgrid , abs.(Kn_new[1,:]), label=L"\left|\sum_{\omega_n > \Lambda} K_n \right|" ) -# xlabel!(L"s") -# #legend() -# #pl = plot(wgrid , abs.(Kn[1,:]) ) -# savefig(pl, "expan_eig.pdf") -# end - -# idx = searchsortedfirst(eig.S./eig.S[1], eps, rev=true) - -# print("rank: $(idx)\n") -# if space == :n -# n_grid = IR(ngrid, eig.U, idx, "omega_n") -# #print("tail selected: left $(ngrid[left] in n_grid) right $(ngrid[right] in n_grid)\n") -# elseif space == :τ -# tau_grid = IR(tgrid, eig.U, idx, "tau", Lambda, true) -# elseif space==:ex -# expan_grid = IR(expangrid, eig.U, idx, "expan") -# end - -# omega_grid = IR(wgrid, eig.V, idx, "omega") - -# #Use implicit fourier to get U in the other space - -# if space == :n -# U = (Ktau * eig.V)[:, 1:idx] -# elseif space == :τ -# U = (Kn * eig.V)[:, 1:idx] -# U_compare = F*eig.U[:, 1:idx] - -# end - -# for i in 1:idx -# U[:, i] = U[:, i] ./ eig.S[i] -# end -# print("fourier error: $(maximum(abs.(U- U_compare)))\n") -# #Un = U -# Un = U_compare -# Un_new = copy(Un[left:right, :]) -# Un_new[1, :] = sum(Un[1:left, :], dims=1) -# Un_new[end, :] = sum(Un[right:end, :], dims=1) -# pl1 = plot(omega , abs.(Un[: , 15]).*omega.^2, linewidth = 1) -# #plot!(pl1, omega[left:right] , abs.(Un_new[1 , 15])* ones(length(omega)), linewidth = 1) -# ylabel!(L"\omega_n^2 U") -# xlabel!(L"\omega") -# savefig(pl1, "U_omega_n_2.pdf") - -# pl2 = plot(omega , abs.(Un[: , 15]), linewidth = 1) -# plot!(pl2, omega , abs.(Un_new[1 , 15])* ones(length(omega)), linewidth = 1, label="sum tail") -# ylabel!(L"U") -# xlabel!(L"\omega") -# savefig(pl2, "U_tail.pdf") - -# pl = plot( collect(1:idx) , maximum(abs.(Un), dims=1)[1,:], linewidth = 1,label = L"max(\left|U_n\right|)") -# plot!(pl,collect(1:idx), abs.(Un_new[1,:]), label=L"\left|\sum_{\omega_n > \Lambda} U_n \right|" ) -# xlabel!(L"s") -# #legend() -# #pl = plot(wgrid , abs.(Kn[1,:]) ) -# savefig(pl, "U.pdf") -# #print("U diff: $(maximum(abs.(U - U_compare)))\n") - -# if space == :n -# tau_grid = IR(tgrid, U, idx, "tau") -# elseif space == :τ -# n_grid = IR(ngrid, U, idx, "omega_n", Lambda, true) -# end -# dlr = DLRGrid(Lambda, beta, eps, true, :none, dtype=T) - -# test_err(dlr, tau_grid, tgrid, t_grid, omega_grid, :τ, regular, omega0) -# return omega_grid, tau_grid, n_grid -# end - - -function generate_grid_expan(eps::T, Lambda::T, n_trunc::Int, space::Symbol=:τ, regular = false, - omega0::T = Lambda,) where {T} - # generate frequency finegrid - w_grid = fine_ωGrid_test(T(Lambda), 24, T(1.5)) - weight_w = zeros(T,length(w_grid)) - #calculate grid weights - for i in 1:length(w_grid) - data = zeros(T,length(w_grid)) - data[i] = 1.0 - weight_w[i] = Interp.integrate1D(data, w_grid) - end - - #symmetrize the grid - wgrid = vcat(-w_grid.grid[end:-1:1], w_grid.grid) - weight_w = vcat(weight_w[end:-1:1], weight_w) - - - #generate tau fine grid - t_grid = fine_τGrid_test(T(Lambda),128, T(1.5)) - weight_t = zeros(T,length(t_grid)) - for i in 1:length(t_grid) - data = zeros(T,length(t_grid)) - data[i] = 1.0 - weight_t[i] = Interp.integrate1D(data, t_grid) - end - tgrid = t_grid.grid - - # generate fine n grid - - #ngrid = nGrid_test(true, T(Lambda), 12, T(1.5)) - ngrid = uni_ngrid(true, T(n_trunc*Lambda)) - omega = (2*ngrid.+1)* π - - expangrid = collect(1:n_trunc) - - #F = Fourier(ngrid, tgrid, weight_t) - #Kn = Kfunc_freq(wgrid, Int.(ngrid), weight_w, regular, omega0) - Ktau = Kfunc(wgrid, tgrid, weight_w, weight_t, regular, omega0) - Kexpan = Kfunc_expan(wgrid, expangrid, weight_w, Lambda) - # Kn_fourier = F*Ktau - # print("fourier error: $(maximum(abs.(Kn- Kn_fourier)))\n") - #left = searchsortedfirst(omega, -n_trunc*Lambda/10) - #right = searchsortedfirst(omega, n_trunc*Lambda/10) - - if space == :τ - eig = svd(Ktau, full = true) - elseif space == :ex - eig = svd(Kexpan, full = true) - print("eig highfreq expansion: $(eig.S[1:10])\n") - pl = plot( collect(1:70), eig.S[1:70], linewidth = 1,label = L"max(\left|K_n\right|)", yaxis=:log) - # plot!(pl,wgrid , abs.(Kn_new[1,:]), label=L"\left|\sum_{\omega_n > \Lambda} K_n \right|" ) - xlabel!(L"s") - #legend() - #pl = plot(wgrid , abs.(Kn[1,:]) ) - savefig(pl, "expan_eig.pdf") - end - - idx = searchsortedfirst(eig.S./eig.S[1], eps, rev=true) - - print("rank: $(idx)\n") - if space == :τ - tau_grid = IR(tgrid, eig.U, idx, "tau", Lambda, true) - elseif space==:ex - expan_grid = IR(expangrid, eig.U, idx, "expan") - end - - omega_grid = IR(wgrid, eig.V, idx, "omega") - - #Use implicit fourier to get U in the other space - - if space == :ex - U = (Ktau * eig.V)[:, 1:idx] - elseif space == :τ - U = (Kexpan * eig.V)[:, 1:idx] - end - - for i in 1:idx - U[:, i] = U[:, i] ./ eig.S[i] - end - - if space == :ex - tau_grid = IR(tgrid, U, idx, "tau") - elseif space == :τ - expan_grid = IR(expangrid, U, idx, "expan", Lambda, true) - end - return omega_grid, tau_grid, expan_grid -end - - -function Fourier(ngrid, taugrid::Vector{T}, tauweight::Vector{T}) where {T} - result = zeros(Complex{T}, (length(ngrid), length(taugrid))) - for i in eachindex(ngrid) - for j in eachindex(taugrid) - result[i, j] = exp(im *(2ngrid[i]+1)* π *taugrid[j]) *(tauweight[j]) - end - end - return result -end - -SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr, type, grid, degree=48, regularized=true) -function MultiPole(dlr, grid, type, coeff, weight=nothing) - Euv = dlr.Euv - poles = coeff * Euv - # return Sample.MultiPole(dlr.β, dlr.isFermi, grid, type, poles, dlr.symmetry; regularized = true) - if isnothing(weight) - return Sample.MultiPole(dlr, type, poles, grid; regularized=true) - else - return Sample.MultiPole(dlr, type, poles, grid, weight; regularized=true) - end -end - -function test_err_dlr(dlr, ir_grid, target_fine_grid, ir_omega_grid, space, target_space, regular, omega0, hasweight=false, weight=nothing) - - - #generate_grid_expan(eps, Lambda, expan_trunc, :ex, false, datatype(Lambda)) - Gsample = SemiCircle(dlr, ir_grid, space) - T = typeof(Gsample[1]) - if space == :n - K = Kfunc_freq(ir_omega_grid , (ir_grid), regular, omega0) - noise = (rand(T, length(Gsample))) - noise = 1e-6 * noise ./ norm(noise) - Gsample += noise - elseif space == :τ - K = Kfunc(ir_omega_grid, ir_grid, regular, omega0) - noise = (rand(T, length(Gsample))) - noise = 1e-6 * noise ./ norm(noise) - Gsample += noise - end - if target_space==:τ - Kfull = Kfunc( ir_omega_grid , target_fine_grid.grid, regular, omega0) - G_analy = SemiCircle(dlr, target_fine_grid.grid, target_space) - else - Kfull = Kfunc_freq( ir_omega_grid , target_fine_grid, regular, omega0) - G_analy = SemiCircle(dlr, target_fine_grid, target_space) - end - # if hasweight - # if space == :n - # Gsample = Gsample .* weight[] - # end - rho = K \ Gsample - G = Kfull * rho - - #interp_err = sqrt(Interp.integrate1D((G - G_analy) .^ 2, target_fine_grid )) - if target_space==:n - interp_err = norm(G - G_analy) - else - interp_err = sqrt(Interp.integrate1D( (abs.(G - G_analy)).^2, target_fine_grid )) - end - #print("condition KIR: $(cond(K))\n") - print("Exact Green err: $(interp_err)\n") - - -end -# function test_err(dlr, ir_grid, fine_grid, target_fine_grid, ir_omega_grid, space, regular, omega0) - -# #generate_grid_expan(eps, Lambda, expan_trunc, :ex, false, datatype(Lambda)) -# Gsample = SemiCircle(dlr, ir_grid, space) -# if space == :n -# K = Kfunc_freq(ir_omega_grid , ir_grid, regular, omega0) -# elseif space == :τ -# K = Kfunc(ir_omega_grid, ir_grid, regular, omega0) -# end -# Ktau = Kfunc( ir_omega_grid , target_fine_grid.grid, regular, omega0) -# rho = K \ Gsample -# G = Ktau * rho -# G_analy = SemiCircle(dlr, target_fine_grid, :τ) -# interp_err = sqrt(Interp.integrate1D((G - G_analy) .^ 2, target_fine_grid )) -# print("Exact Green err: $(interp_err)\n") - - -# end -# if abspath(PROGRAM_FILE) == @__FILE__ -# # dlr = DLRGrid(Euv=lambda, β=beta, isFermi=true, rtol=1e-12, symmetry=:sym) - -# datatype = Float64 -# #setprecision(128) -# #atatype = BigFloat -# isFermi = true -# symmetry = :none -# beta = datatype(1.0) -# Lambda = datatype(1000) -# eps = datatype(1e-10) -# n_trunc = datatype(10) #omega_n is truncated at n_trunc * Lambda -# expan_trunc = 1000 -# omega_grid, tau_grid, n_grid = generate_grid(eps, Lambda, n_trunc, :τ, false, datatype(Lambda)) - -# #generate_grid_resum(eps, Lambda, n_trunc, false, datatype(Lambda)) -# end diff --git a/test/kernel_svd_sym.jl b/test/kernel_svd_sym.jl deleted file mode 100644 index a63073d..0000000 --- a/test/kernel_svd_sym.jl +++ /dev/null @@ -1,442 +0,0 @@ -using Lehmann -using Printf -using CompositeGrids -using LinearAlgebra -using GenericLinearAlgebra -using Random -using Plots -using LaTeXStrings -include("grid.jl") -include("symQR.jl") -""" -composite expoential grid -""" -function Htau(tau::Vector{T}, weight::Vector{T}, gamma) where {T} - result = zeros(T, (length(tau), length(tau))) - for i in eachindex(tau) - for j in eachindex(tau) - result[i, j] = sqrt(weight[j] * weight[i]) * (exp(-(tau[i] + tau[j])) - exp(-(tau[i] + tau[j]) * gamma)) / (tau[i] + tau[j]) - end - end - #print(result) - return result -end - - -@inline function F1(a::T, b::T) where {T} - if abs(a + b) > EPS - return (1 - exp(-(a + b))) / (a + b) - else - return T(1-(a+b)/2 + (a+b)^2/6 - (a+b)^3/24) - end -end - -""" -``G(x, y) = (exp(-x)-exp(-y))/(x-y)`` -``G(x, x) = -exp(-x)`` -""" -@inline function G1(a::T, b::T) where {T} - if abs(a - b) > EPS - return (exp(-a) - exp(-b)) / (b - a) - else - return (exp(-a) + exp(-b)) / 2 - end -end - -function Homega(omega::Vector{T}, weight::Vector{T}) where {T} - result = zeros(T, (length(omega), length(omega))) - for i in eachindex(omega) - for j in eachindex(omega) - if omega[i]*omega[j]>0 - #result[i, j] = sqrt(weight[j] * weight[i]) * F1(abs(omega[i]), abs(omega[j])) - result[i, j] = F1(abs(omega[i]), abs(omega[j])) - else - result[i, j] = G1(abs(omega[i]), abs(omega[j])) - end - end - end - #print(result) - return result -end - -function IntFermiT(omega::T) where {T} - omega_new = omega/2 - if omega_new < 1e-6 - return 0.5*(1 - omega_new^2/3 + 2omega_new^4/15 - 17omega_new^6/315) - else - return tanh(omega_new)/omega - end -end - - -function Kfunc(omega::Vector{T}, tau::Vector{T}, weight_omega::Vector{T}, weight_tau::Vector{T} , ifregular::Bool=false, omega0::T = 1e-4) where {T} - result = zeros(T, (length(tau), length(omega))) - omega0 = T(omega0) - for i in eachindex(tau) - for j in eachindex(omega) - #result[i, j] = sqrt(weight_tau[i] * weight_omega[j]) * Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - #result[i, j] = sqrt(weight_tau[i] * weight_omega[j]) * Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - if ifregular - #result[i, j] = Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - 1/2.0 # - IntFermiT(omega[j]) - result[i, j] =sqrt(weight_tau[i])*(Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - 1.0/2.0 - (1- 2*tau[i])*omega[j]/4.0 - (tau[i]-1)*tau[i]*omega[j]^2/4.0) - - #(Spectral.kernelSymT(tau[i], omega0, T(1.0)) - #+ Spectral.kernelSymT(1.0 - tau[i], omega0, T(1.0)))/2.0 - else - #result[i, j] = Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - result[i, j] = sqrt(weight_tau[i])*sqrt(weight_omega[j])*Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - #result[i, j] = sqrt(weight_tau[i])*Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - #result[i, j] = sqrt(weight_omega[j])*Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - #result[i, j] = Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - end - #result[i, j] = Spectral.kernelSymT_PH(tau[i], omega[j], T(1.0)) - #result[i, j] = weight_tau[i] * weight_omega[j]*kernelSymT_test(tau[i], omega[j], T(1.0)) - #result[i, j] = sqrt(weight_omega[i] * weight_tau[j]) * - end - end - #print(result) - return result -end - -function Kfunc(omega::Vector{T}, tau::Vector{T}, ifregular::Bool=false, omega0::T = 1e-4) where {T} - result = zeros(T, (length(tau), length(omega))) - omega0 = T(omega0) - for i in eachindex(tau) - for j in eachindex(omega) - if ifregular - - result[i, j] =Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - 1.0/2.0 - (1- 2*tau[i])*omega[j]/4.0 - (tau[i]-1)*tau[i]*omega[j]^2/4.0 - - else - #result[i, j] = Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - result[i, j] = Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - - end - - end - end - - return result -end - - -function Kfunc_freq(omega::Vector{T}, n::Vector{Int}, weight_omega::Vector{T}, regular::Bool=false ,omega0::T=1e-4 ) where {T} - result = zeros(Complex{T}, (length(n), length(omega))) - omega0 = T(omega0) - for i in eachindex(n) - omega_n = (2*n[i]+1)* π - for j in eachindex(omega) - if regular - result[i, j] =(Spectral.kernelFermiSymΩ(n[i], omega[j], T(1.0)) + 1/(im*omega_n) + omega[j]/(im*omega_n)^2 + omega[j]^2/(im*omega_n)^3) - else - result[i, j] = sqrt(weight_omega[j])*Spectral.kernelFermiSymΩ(n[i], omega[j], T(1.0)) - #result[i, j] = Spectral.kernelFermiSymΩ(n[i], omega[j], T(1.0)) - end - end - end - return result -end - - -function Kfunc_freq(omega::Vector{T}, n::Vector{Int}, regular::Bool=false ,omega0::T=1e-4 ) where {T} - result = zeros(Complex{T}, (length(n), length(omega))) - omega0 = T(omega0) - for i in eachindex(n) - omega_n = (2*n[i]+1)* π - for j in eachindex(omega) - if regular - result[i, j] =(Spectral.kernelFermiSymΩ(n[i], omega[j], T(1.0)) + 1/(im*omega_n) + omega[j]/(im*omega_n)^2 + omega[j]^2/(im*omega_n)^3) - else - result[i, j] =Spectral.kernelFermiSymΩ(n[i], omega[j], T(1.0)) - end - end - end - return result -end - -function Kfunc_freq!(result::Matrix{Complex{T}}, omega::Vector{T}, n::Vector{Int}, regular::Bool=false ,omega0::T=1e-4 ) where {T} - omega0 = T(omega0) - for i in eachindex(n) - omega_n = (2*n[i]+1)* π - for j in eachindex(omega) - if regular - result[i, j] =(Spectral.kernelFermiSymΩ(n[i], omega[j], T(1.0)) + 1/(im*omega_n) + omega[j]/(im*omega_n)^2 + omega[j]^2/(im*omega_n)^3) - else - result[i, j] =Spectral.kernelFermiSymΩ(n[i], omega[j], T(1.0)) - end - end - end -end - -function Kfunc_expan(omega::Vector{T}, n::Vector{Int}, weight_omega::Vector{T}, Lambda::T) where {T} - result = zeros(T, (length(n), length(omega))) - for i in eachindex(n) - for j in eachindex(omega) - result[i, j] = (omega[j]/Lambda)^(i-1)*weight_omega[j] - end - end - return result -end - -function symmetrize_idx(idx_list) - result = Int[] - for idx in idx_list - if !(idx in result) - append!(result, idx) - append!(result, length(idx_list) - idx +1) - end - end - @assert length(idx_list) == length(result) - return result -end - - -function IR(grid, U, idx, name, Lambda=100, ifplot = false) - Un = U[:, 1:idx] - qr_Un = qr(transpose(Un), Val(true)) - qr_nidx = qr_Un.p - - left = searchsortedfirst(grid, -Lambda) - right = searchsortedfirst(grid, Lambda) - #print("selected: $(length(qr_nidx)) $(minimum(qr_nidx[1:idx])) $(maximum(qr_nidx[1:idx]))\n" ) - if name == "omega_n" - shift = 0#20 - idx - else - shift = 0 - end - ir_grid = sort(grid[qr_nidx[1:idx+shift]]) - # Unew = Un[sort(qr_nidx[1:idx+shift]),:] - # Hnew = Unew*transpose(conj(Unew)) - # Hfull = Un*transpose(conj(Un)) - # #print("eigenvalue: $(eigvals(Hfull))\n") - # if ifplot - # #for plidx in [1,5,10, 15, 20, 25, 30] - # #plot!(pl, grid[qr_nidx[1:idx]] , abs.(Un[qr_nidx[1:idx],plidx]), seriestype=:scatter, markershape=:circle) - - # #end - # pl = plot() - # if name == "omega_n" - # heatmap!(pl, abs.(Unew), title=L"U_{IR}", xlabel=L"s", ylabel=L"\omega_n", color=:viridis) - # else - # heatmap!(pl, abs.(Unew), title=L"U_{IR}", xlabel=L"s", ylabel=L"\tau", color=:viridis) - # end - # #xlabel!(L"\omega") - # #legend() - # #pl = plot(wgrid , abs.(Kn[1,:]) ) - # savefig(pl, name*"UIR.pdf") - # if name == "omega_n" - # pl = heatmap(abs.(Un[left:right, :]), title=L"U_{full}", xlabel=L"s", ylabel=L"\tau", color=:viridis) - # else - # pl = heatmap(abs.(Un), title=L"U_{full}", xlabel=L"s", ylabel=L"\tau", color=:viridis) - # end - # diag_full = diag(abs.(Hfull)) - # # print("max Ufull:$(maximum(diag_full))\n") - # # print("max Ufull:$(minimum(diag_full))\n") - # diag_IR = diag(abs.(Hnew)) - # # print("max UIR:$(maximum(diag_IR))\n") - # # print("max UIR:$(minimum(diag_IR))\n") - # savefig(pl, name*"Ufull.pdf") - # end - #print(name*" R cond :$(cond(qr_Un.R[:, qr_nidx[1:idx]]))\n") - #print(name*" Q cond :$(cond(qr_Un.Q))\n") - #print(name*" IR cond :$(cond(Un[qr_nidx[1:idx] , 1:idx]))\n") - #print(name*" Full cond:$(cond(Un[:, 1:idx]))\n") - #print(name*" IR cond :$(cond(Unew))\n") - #print(name*" Full cond:$(cond(Un))\n") - # print(name*" IR cond UUT:$(cond(abs.(Hnew)))\n") - # print(name*" Full cond UUT:$(cond(abs.(Hfull)))\n") - return qr_nidx, ir_grid -end - -function generate_grid(eps::T, Lambda::T, n_trunc::T, space::Symbol=:τ, regular = false, - omega0::T = Lambda,) where {T} - # generate frequency finegrid - w_grid = fine_ωGrid_test(T(Lambda), 24, T(1.5)) - weight_w = zeros(T,length(w_grid)) - #calculate grid weights - for i in 1:length(w_grid) - data = zeros(T,length(w_grid)) - data[i] = 1.0 - weight_w[i] = Interp.integrate1D(data, w_grid) - end - - #symmetrize the grid - wgrid = vcat(-w_grid.grid[end:-1:1], w_grid.grid) - weight_w = vcat(weight_w[end:-1:1], weight_w) - - - #generate tau fine grid - t_grid = fine_τGrid_test(T(Lambda),128, T(1.5)) - weight_t = zeros(T,length(t_grid)) - for i in 1:length(t_grid) - data = zeros(T,length(t_grid)) - data[i] = 1.0 - weight_t[i] = Interp.integrate1D(data, t_grid) - end - tgrid = t_grid.grid - - # generate fine n grid - - ngrid = nGrid_test(true, T(Lambda), 12, T(1.5)) - #ngrid = uni_ngrid(true, T(n_trunc*Lambda)) - omega = (2*ngrid.+1)* π - - #ngrid = vcat(ngrid, dlr.n) - #unique!(ngrid) - #ngrid = sort(ngrid) - - #regular controls if we add 1/(iwn - Lambda) term to compensate the tail - Kn = Kfunc_freq(wgrid, Int.(ngrid), weight_w, regular, omega0) - Ktau = Kfunc(wgrid, tgrid, weight_w, weight_t, regular, omega0) - # Kn_fourier = F*Ktau - # print("fourier error: $(maximum(abs.(Kn- Kn_fourier)))\n") - left = searchsortedfirst(omega, -n_trunc*Lambda/10) - right = searchsortedfirst(omega, n_trunc*Lambda/10) - # print("$(left) $(right)\n") - # Kn_new = copy(Kn[left:right, :]) - # Kn_new[1, :] = sum(Kn[1:left, :], dims=1) - # Kn_new[end, :] = sum(Kn[right:end, :], dims=1) - - # #print("$(maximum(abs.(Kn), dims=1))\n $(abs.(Kn[1,:]))\n $(abs.(Kn[end,:])) \n") - - # Kn = Kn_new - - if space == :n - eig = svd(Kn, full = true) - elseif space == :τ - eig = svd(Ktau, full = true) - end - - idx = searchsortedfirst(eig.S./eig.S[1], eps, rev=true) - #idx = idx-5 - print("rank: $(idx)\n") - if space == :n - pivoted_idx, n_grid = IR(ngrid, eig.U, idx, "omega_n") - #print("tail selected: left $(ngrid[left] in n_grid) right $(ngrid[right] in n_grid)\n") - elseif space == :τ - pivoted_idx, tau_grid = IR(tgrid, eig.U, idx, "tau", Lambda, true) - end - - pivoted_idx, omega_grid = IR(wgrid, eig.V, idx, "omega") - - #Use implicit fourier to get U in the other space - - if space == :n - U = (Ktau * eig.V)[:, 1:idx] - elseif space == :τ - U = (Kn * eig.V)[:, 1:idx] - end - - for i in 1:idx - U[:, i] = U[:, i] ./ eig.S[i] - end - - if space == :n - pivoted_idx, tau_grid = IR(tgrid, U, idx, "tau") - elseif space == :τ - pivoted_idx, n_grid = IR(ngrid, U, idx, "omega_n", Lambda, true) - end - return omega_grid, tau_grid, n_grid -end - - - -function Fourier(ngrid, taugrid::Vector{T}, tauweight::Vector{T}) where {T} - result = zeros(Complex{T}, (length(ngrid), length(taugrid))) - for i in eachindex(ngrid) - for j in eachindex(taugrid) - result[i, j] = exp(im *(2ngrid[i]+1)* π *taugrid[j]) *(tauweight[j]) - end - end - return result -end - -SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr, type, grid, degree=48, regularized=true) -function MultiPole(dlr, grid, type, coeff, weight=nothing) - Euv = dlr.Euv - poles = coeff * Euv - # return Sample.MultiPole(dlr.β, dlr.isFermi, grid, type, poles, dlr.symmetry; regularized = true) - if isnothing(weight) - return Sample.MultiPole(dlr, type, poles, grid; regularized=true) - else - return Sample.MultiPole(dlr, type, poles, grid, weight; regularized=true) - end -end - -function test_err_dlr(dlr, ir_grid, target_fine_grid, ir_omega_grid, space, target_space, regular, omega0, hasweight=false, weight=nothing) - - - #generate_grid_expan(eps, Lambda, expan_trunc, :ex, false, datatype(Lambda)) - Gsample = SemiCircle(dlr, ir_grid, space) - T = typeof(Gsample[1]) - if space == :n - K = Kfunc_freq(ir_omega_grid , (ir_grid), regular, omega0) - noise = (rand(T, length(Gsample))) - noise = 1e-6 * noise ./ norm(noise) - Gsample += noise - elseif space == :τ - K = Kfunc(ir_omega_grid, ir_grid, regular, omega0) - noise = (rand(T, length(Gsample))) - noise = 1e-6 * noise ./ norm(noise) - Gsample += noise - end - if target_space==:τ - Kfull = Kfunc( ir_omega_grid , target_fine_grid.grid, regular, omega0) - G_analy = SemiCircle(dlr, target_fine_grid.grid, target_space) - else - Kfull = Kfunc_freq( ir_omega_grid , target_fine_grid, regular, omega0) - G_analy = SemiCircle(dlr, target_fine_grid, target_space) - end - # if hasweight - # if space == :n - # Gsample = Gsample .* weight[] - # end - rho = K \ Gsample - G = Kfull * rho - - #interp_err = sqrt(Interp.integrate1D((G - G_analy) .^ 2, target_fine_grid )) - if target_space==:n - interp_err = norm(G - G_analy) - else - interp_err = sqrt(Interp.integrate1D( (abs.(G - G_analy)).^2, target_fine_grid )) - end - #print("condition KIR: $(cond(K))\n") - print("Exact Green err: $(interp_err)\n") - - -end -# function test_err(dlr, ir_grid, fine_grid, target_fine_grid, ir_omega_grid, space, regular, omega0) - -# #generate_grid_expan(eps, Lambda, expan_trunc, :ex, false, datatype(Lambda)) -# Gsample = SemiCircle(dlr, ir_grid, space) -# if space == :n -# K = Kfunc_freq(ir_omega_grid , ir_grid, regular, omega0) -# elseif space == :τ -# K = Kfunc(ir_omega_grid, ir_grid, regular, omega0) -# end -# Ktau = Kfunc( ir_omega_grid , target_fine_grid.grid, regular, omega0) -# rho = K \ Gsample -# G = Ktau * rho -# G_analy = SemiCircle(dlr, target_fine_grid, :τ) -# interp_err = sqrt(Interp.integrate1D((G - G_analy) .^ 2, target_fine_grid )) -# print("Exact Green err: $(interp_err)\n") - - -# end -# if abspath(PROGRAM_FILE) == @__FILE__ -# # dlr = DLRGrid(Euv=lambda, β=beta, isFermi=true, rtol=1e-12, symmetry=:sym) - -# datatype = Float64 -# #setprecision(128) -# #atatype = BigFloat -# isFermi = true -# symmetry = :none -# beta = datatype(1.0) -# Lambda = datatype(1000) -# eps = datatype(1e-10) -# n_trunc = datatype(10) #omega_n is truncated at n_trunc * Lambda -# expan_trunc = 1000 -# omega_grid, tau_grid, n_grid = generate_grid(eps, Lambda, n_trunc, :τ, false, datatype(Lambda)) - -# #generate_grid_resum(eps, Lambda, n_trunc, false, datatype(Lambda)) -# end diff --git a/test/kernel_svd_wh.jl b/test/kernel_svd_wh.jl deleted file mode 100644 index e161075..0000000 --- a/test/kernel_svd_wh.jl +++ /dev/null @@ -1,528 +0,0 @@ -using Lehmann -using Printf -using CompositeGrids -using LinearAlgebra -using GenericLinearAlgebra -using Random -using Plots -using LaTeXStrings -include("grid.jl") -""" -composite expoential grid -""" -function Htau(tau::Vector{T}, weight::Vector{T}, gamma) where {T} - result = zeros(T, (length(tau), length(tau))) - for i in eachindex(tau) - for j in eachindex(tau) - result[i, j] = sqrt(weight[j] * weight[i]) * (exp(-(tau[i] + tau[j])) - exp(-(tau[i] + tau[j]) * gamma)) / (tau[i] + tau[j]) - end - end - #print(result) - return result -end - - -@inline function F1(a::T, b::T) where {T} - if abs(a + b) > EPS - return (1 - exp(-(a + b))) / (a + b) - else - return T(1-(a+b)/2 + (a+b)^2/6 - (a+b)^3/24) - end -end - -""" -``G(x, y) = (exp(-x)-exp(-y))/(x-y)`` -``G(x, x) = -exp(-x)`` -""" -@inline function G1(a::T, b::T) where {T} - if abs(a - b) > EPS - return (exp(-a) - exp(-b)) / (b - a) - else - return (exp(-a) + exp(-b)) / 2 - end -end - -function Homega(omega::Vector{T}, weight::Vector{T}) where {T} - result = zeros(T, (length(omega), length(omega))) - for i in eachindex(omega) - for j in eachindex(omega) - if omega[i]*omega[j]>0 - #result[i, j] = sqrt(weight[j] * weight[i]) * F1(abs(omega[i]), abs(omega[j])) - result[i, j] = F1(abs(omega[i]), abs(omega[j])) - else - result[i, j] = G1(abs(omega[i]), abs(omega[j])) - end - end - end - #print(result) - return result -end - -function IntFermiT(omega::T) where {T} - omega_new = omega/2 - if omega_new < 1e-6 - return 0.5*(1 - omega_new^2/3 + 2omega_new^4/15 - 17omega_new^6/315) - else - return tanh(omega_new)/omega - end -end - - -function Kfunc(omega::Vector{T}, tau::Vector{T}, weight_omega::Vector{T}, weight_tau::Vector{T}) where {T} - result = zeros(T, (length(tau), length(omega))) - for i in eachindex(tau) - for j in eachindex(omega) - result[i, j] = weight_omega[j]*weight_tau[i]*Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - end - end - return result -end - -function Kfunc(omega::Vector{T}, tau::Vector{T}) where {T} - result = zeros(T, (length(tau), length(omega))) - for i in eachindex(tau) - for j in eachindex(omega) - result[i, j] = Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - end - end - return result -end - -function Kfunc(omega::Vector{T}, tau::Vector{T}, weight_tau::Vector{T}) where {T} - result = zeros(T, (length(tau), length(omega))) - for i in eachindex(tau) - for j in eachindex(omega) - result[i, j] = weight_tau[i]*Spectral.kernelSymT(tau[i], omega[j], T(1.0)) - end - end - return result -end - - -function Kfunc_freq(omega::Vector{T}, n::Vector{Int}, weight_omega::Vector{T}) where {T} - result = zeros(Complex{T}, (length(n), length(omega))) - for i in eachindex(n) - for j in eachindex(omega) - result[i, j] =weight_omega[j]*Spectral.kernelFermiSymΩ(n[i], omega[j], T(1.0)) - end - end - return result -end - - -function Kfunc_freq(omega::Vector{T}, n::Vector{Int}) where {T} - result = zeros(Complex{T}, (length(n), length(omega))) - for i in eachindex(n) - for j in eachindex(omega) - result[i, j] =Spectral.kernelFermiSymΩ(n[i], omega[j], T(1.0)) - end - end - return result -end - -function Kfunc_freq!(result::Matrix{Complex{T}}, omega::Vector{T}, n::Vector{Int}, regular::Bool=false ,omega0::T=1e-4 ) where {T} - omega0 = T(omega0) - for i in eachindex(n) - omega_n = (2*n[i]+1)* π - for j in eachindex(omega) - if regular - result[i, j] =(Spectral.kernelFermiSymΩ(n[i], omega[j], T(1.0)) + 1/(im*omega_n) + omega[j]/(im*omega_n)^2 + omega[j]^2/(im*omega_n)^3) - else - result[i, j] =Spectral.kernelFermiSymΩ(n[i], omega[j], T(1.0)) - end - end - end -end - -function Kfunc_expan(omega::Vector{T}, n::Vector{Int}, weight_omega::Vector{T}, Lambda::T) where {T} - result = zeros(T, (length(n), length(omega))) - for i in eachindex(n) - for j in eachindex(omega) - result[i, j] = (omega[j]/Lambda)^(i-1)*weight_omega[j] - end - end - return result -end - - -function IR(grid, U, idx, name, Lambda=100, ifplot = false) - Un = U[:, 1:idx] - qr_Un = qr(transpose(Un), Val(true)) - qr_nidx = qr_Un.p - left = searchsortedfirst(grid, -Lambda) - right = searchsortedfirst(grid, Lambda) - #print("selected: $(length(qr_nidx)) $(minimum(qr_nidx[1:idx])) $(maximum(qr_nidx[1:idx]))\n" ) - if name == "omega_n" - shift = 0#20 - idx - else - shift = 0 - end - ir_grid = sort(grid[qr_nidx[1:idx+shift]]) - Unew = Un[sort(qr_nidx[1:idx+shift]),:] - Hnew = Unew*transpose(conj(Unew)) - Hfull = Un*transpose(conj(Un)) - #print("eigenvalue: $(eigvals(Hfull))\n") - if ifplot - #for plidx in [1,5,10, 15, 20, 25, 30] - #plot!(pl, grid[qr_nidx[1:idx]] , abs.(Un[qr_nidx[1:idx],plidx]), seriestype=:scatter, markershape=:circle) - - #end - pl = plot() - if name == "omega_n" - heatmap!(pl, abs.(Unew), title=L"U_{IR}", xlabel=L"s", ylabel=L"\omega_n", color=:viridis) - else - heatmap!(pl, abs.(Unew), title=L"U_{IR}", xlabel=L"s", ylabel=L"\tau", color=:viridis) - end - #xlabel!(L"\omega") - #legend() - #pl = plot(wgrid , abs.(Kn[1,:]) ) - savefig(pl, name*"UIR.pdf") - if name == "omega_n" - pl = heatmap(abs.(Un[left:right, :]), title=L"U_{full}", xlabel=L"s", ylabel=L"\tau", color=:viridis) - else - pl = heatmap(abs.(Un), title=L"U_{full}", xlabel=L"s", ylabel=L"\tau", color=:viridis) - end - diag_full = diag(abs.(Hfull)) - # print("max Ufull:$(maximum(diag_full))\n") - # print("max Ufull:$(minimum(diag_full))\n") - diag_IR = diag(abs.(Hnew)) - # print("max UIR:$(maximum(diag_IR))\n") - # print("max UIR:$(minimum(diag_IR))\n") - savefig(pl, name*"Ufull.pdf") - end - #print(name*" R cond :$(cond(qr_Un.R[:, qr_nidx[1:idx]]))\n") - #print(name*" Q cond :$(cond(qr_Un.Q))\n") - #print(name*" IR cond :$(cond(Un[qr_nidx[1:idx] , 1:idx]))\n") - #print(name*" Full cond:$(cond(Un[:, 1:idx]))\n") - #print(name*" IR cond :$(cond(Unew))\n") - #print(name*" Full cond:$(cond(Un))\n") - # print(name*" IR cond UUT:$(cond(abs.(Hnew)))\n") - # print(name*" Full cond UUT:$(cond(abs.(Hfull)))\n") - return qr_nidx, ir_grid -end - -# function generate_grid(eps::T, Lambda::T, n_trunc::T, space::Symbol=:τ, regular = false, -# omega0::T = Lambda,) where {T} -# # generate frequency finegrid -# w_grid = fine_ωGrid_test(T(Lambda), 24, T(1.5)) -# weight_w = zeros(T,length(w_grid)) -# #calculate grid weights -# for i in 1:length(w_grid) -# data = zeros(T,length(w_grid)) -# data[i] = 1.0 -# weight_w[i] = Interp.integrate1D(data, w_grid) -# end - -# #symmetrize the grid -# wgrid = vcat(-w_grid.grid[end:-1:1], w_grid.grid) -# weight_w = vcat(weight_w[end:-1:1], weight_w) - - -# #generate tau fine grid -# t_grid = fine_τGrid_test(T(Lambda),128, T(1.5)) -# weight_t = zeros(T,length(t_grid)) -# for i in 1:length(t_grid) -# data = zeros(T,length(t_grid)) -# data[i] = 1.0 -# weight_t[i] = Interp.integrate1D(data, t_grid) -# end -# tgrid = t_grid.grid - -# # generate fine n grid - -# #ngrid = nGrid_test(true, T(Lambda), 12, T(1.5)) -# ngrid = uni_ngrid(true, T(n_trunc*Lambda)) -# omega = (2*ngrid.+1)* π - -# #ngrid = vcat(ngrid, dlr.n) -# #unique!(ngrid) -# #ngrid = sort(ngrid) - -# #regular controls if we add 1/(iwn - Lambda) term to compensate the tail -# expangrid = collect(1:10000) - -# F = Fourier(ngrid, tgrid, weight_t) -# Kn = Kfunc_freq(wgrid, Int.(ngrid), weight_w, regular, omega0) -# Ktau = Kfunc(wgrid, tgrid, weight_w, weight_t, regular, omega0) -# Kexpan = Kfunc_expan(wgrid, expangrid, weight_w, Lambda) -# # Kn_fourier = F*Ktau -# # print("fourier error: $(maximum(abs.(Kn- Kn_fourier)))\n") -# left = searchsortedfirst(omega, -n_trunc*Lambda/10) -# right = searchsortedfirst(omega, n_trunc*Lambda/10) -# # print("$(left) $(right)\n") -# # Kn_new = copy(Kn[left:right, :]) -# # Kn_new[1, :] = sum(Kn[1:left, :], dims=1) -# # Kn_new[end, :] = sum(Kn[right:end, :], dims=1) - -# # #print("$(maximum(abs.(Kn), dims=1))\n $(abs.(Kn[1,:]))\n $(abs.(Kn[end,:])) \n") - -# # Kn = Kn_new - -# if space == :n -# eig = svd(Kn, full = true) -# elseif space == :τ -# eig = svd(Ktau, full = true) -# elseif space == :ex -# eig = svd(Kexpan, full = true) -# print("eig highfreq expansion: $(eig.S[1:10])\n") -# pl = plot( collect(1:70), eig.S[1:70], linewidth = 1,label = L"max(\left|K_n\right|)", yaxis=:log) -# # plot!(pl,wgrid , abs.(Kn_new[1,:]), label=L"\left|\sum_{\omega_n > \Lambda} K_n \right|" ) -# xlabel!(L"s") -# #legend() -# #pl = plot(wgrid , abs.(Kn[1,:]) ) -# savefig(pl, "expan_eig.pdf") -# end - -# idx = searchsortedfirst(eig.S./eig.S[1], eps, rev=true) - -# print("rank: $(idx)\n") -# if space == :n -# n_grid = IR(ngrid, eig.U, idx, "omega_n") -# #print("tail selected: left $(ngrid[left] in n_grid) right $(ngrid[right] in n_grid)\n") -# elseif space == :τ -# tau_grid = IR(tgrid, eig.U, idx, "tau", Lambda, true) -# elseif space==:ex -# expan_grid = IR(expangrid, eig.U, idx, "expan") -# end - -# omega_grid = IR(wgrid, eig.V, idx, "omega") - -# #Use implicit fourier to get U in the other space - -# if space == :n -# U = (Ktau * eig.V)[:, 1:idx] -# elseif space == :τ -# U = (Kn * eig.V)[:, 1:idx] -# U_compare = F*eig.U[:, 1:idx] - -# end - -# for i in 1:idx -# U[:, i] = U[:, i] ./ eig.S[i] -# end -# print("fourier error: $(maximum(abs.(U- U_compare)))\n") -# #Un = U -# Un = U_compare -# Un_new = copy(Un[left:right, :]) -# Un_new[1, :] = sum(Un[1:left, :], dims=1) -# Un_new[end, :] = sum(Un[right:end, :], dims=1) -# pl1 = plot(omega , abs.(Un[: , 15]).*omega.^2, linewidth = 1) -# #plot!(pl1, omega[left:right] , abs.(Un_new[1 , 15])* ones(length(omega)), linewidth = 1) -# ylabel!(L"\omega_n^2 U") -# xlabel!(L"\omega") -# savefig(pl1, "U_omega_n_2.pdf") - -# pl2 = plot(omega , abs.(Un[: , 15]), linewidth = 1) -# plot!(pl2, omega , abs.(Un_new[1 , 15])* ones(length(omega)), linewidth = 1, label="sum tail") -# ylabel!(L"U") -# xlabel!(L"\omega") -# savefig(pl2, "U_tail.pdf") - -# pl = plot( collect(1:idx) , maximum(abs.(Un), dims=1)[1,:], linewidth = 1,label = L"max(\left|U_n\right|)") -# plot!(pl,collect(1:idx), abs.(Un_new[1,:]), label=L"\left|\sum_{\omega_n > \Lambda} U_n \right|" ) -# xlabel!(L"s") -# #legend() -# #pl = plot(wgrid , abs.(Kn[1,:]) ) -# savefig(pl, "U.pdf") -# #print("U diff: $(maximum(abs.(U - U_compare)))\n") - -# if space == :n -# tau_grid = IR(tgrid, U, idx, "tau") -# elseif space == :τ -# n_grid = IR(ngrid, U, idx, "omega_n", Lambda, true) -# end -# dlr = DLRGrid(Lambda, beta, eps, true, :none, dtype=T) - -# test_err(dlr, tau_grid, tgrid, t_grid, omega_grid, :τ, regular, omega0) -# return omega_grid, tau_grid, n_grid -# end - - -function generate_grid_expan(eps::T, Lambda::T, n_trunc::Int, space::Symbol=:τ, regular = false, - omega0::T = Lambda,) where {T} - # generate frequency finegrid - w_grid = fine_ωGrid_test(T(Lambda), 24, T(1.5)) - weight_w = zeros(T,length(w_grid)) - #calculate grid weights - for i in 1:length(w_grid) - data = zeros(T,length(w_grid)) - data[i] = 1.0 - weight_w[i] = Interp.integrate1D(data, w_grid) - end - - #symmetrize the grid - wgrid = vcat(-w_grid.grid[end:-1:1], w_grid.grid) - weight_w = vcat(weight_w[end:-1:1], weight_w) - - - #generate tau fine grid - t_grid = fine_τGrid_test(T(Lambda),128, T(1.5)) - weight_t = zeros(T,length(t_grid)) - for i in 1:length(t_grid) - data = zeros(T,length(t_grid)) - data[i] = 1.0 - weight_t[i] = Interp.integrate1D(data, t_grid) - end - tgrid = t_grid.grid - - # generate fine n grid - - #ngrid = nGrid_test(true, T(Lambda), 12, T(1.5)) - ngrid = uni_ngrid(true, T(n_trunc*Lambda)) - omega = (2*ngrid.+1)* π - - expangrid = collect(1:n_trunc) - - #F = Fourier(ngrid, tgrid, weight_t) - #Kn = Kfunc_freq(wgrid, Int.(ngrid), weight_w, regular, omega0) - Ktau = Kfunc(wgrid, tgrid, weight_w, weight_t, regular, omega0) - Kexpan = Kfunc_expan(wgrid, expangrid, weight_w, Lambda) - # Kn_fourier = F*Ktau - # print("fourier error: $(maximum(abs.(Kn- Kn_fourier)))\n") - #left = searchsortedfirst(omega, -n_trunc*Lambda/10) - #right = searchsortedfirst(omega, n_trunc*Lambda/10) - - if space == :τ - eig = svd(Ktau, full = true) - elseif space == :ex - eig = svd(Kexpan, full = true) - print("eig highfreq expansion: $(eig.S[1:10])\n") - pl = plot( collect(1:70), eig.S[1:70], linewidth = 1,label = L"max(\left|K_n\right|)", yaxis=:log) - # plot!(pl,wgrid , abs.(Kn_new[1,:]), label=L"\left|\sum_{\omega_n > \Lambda} K_n \right|" ) - xlabel!(L"s") - #legend() - #pl = plot(wgrid , abs.(Kn[1,:]) ) - savefig(pl, "expan_eig.pdf") - end - - idx = searchsortedfirst(eig.S./eig.S[1], eps, rev=true) - - print("rank: $(idx)\n") - if space == :τ - tau_grid = IR(tgrid, eig.U, idx, "tau", Lambda, true) - elseif space==:ex - expan_grid = IR(expangrid, eig.U, idx, "expan") - end - - omega_grid = IR(wgrid, eig.V, idx, "omega") - - #Use implicit fourier to get U in the other space - - if space == :ex - U = (Ktau * eig.V)[:, 1:idx] - elseif space == :τ - U = (Kexpan * eig.V)[:, 1:idx] - end - - for i in 1:idx - U[:, i] = U[:, i] ./ eig.S[i] - end - - if space == :ex - tau_grid = IR(tgrid, U, idx, "tau") - elseif space == :τ - expan_grid = IR(expangrid, U, idx, "expan", Lambda, true) - end - return omega_grid, tau_grid, expan_grid -end - - -function Fourier(ngrid, taugrid::Vector{T}, tauweight::Vector{T}) where {T} - result = zeros(Complex{T}, (length(ngrid), length(taugrid))) - for i in eachindex(ngrid) - for j in eachindex(taugrid) - result[i, j] = exp(im *(2ngrid[i]+1)* π *taugrid[j]) *(tauweight[j]) - end - end - return result -end - -SemiCircle(dlr, grid, type) = Sample.SemiCircle(dlr, type, grid, degree=48, regularized=true) -function MultiPole(dlr, grid, type, coeff, weight=nothing) - Euv = dlr.Euv - poles = coeff * Euv - # return Sample.MultiPole(dlr.β, dlr.isFermi, grid, type, poles, dlr.symmetry; regularized = true) - if isnothing(weight) - return Sample.MultiPole(dlr, type, poles, grid; regularized=true) - else - return Sample.MultiPole(dlr, type, poles, grid, weight; regularized=true) - end -end - -function test_err_dlr(dlr, ir_grid, target_fine_grid, ir_omega_grid, space, target_space, regular, omega0, hasweight=false, weight=nothing) - - - #generate_grid_expan(eps, Lambda, expan_trunc, :ex, false, datatype(Lambda)) - Gsample = SemiCircle(dlr, ir_grid, space) - T = typeof(Gsample[1]) - if space == :n - K = Kfunc_freq(ir_omega_grid , (ir_grid), regular, omega0) - noise = (rand(T, length(Gsample))) - noise = 1e-6 * noise ./ norm(noise) - Gsample += noise - elseif space == :τ - K = Kfunc(ir_omega_grid, ir_grid, regular, omega0) - noise = (rand(T, length(Gsample))) - noise = 1e-6 * noise ./ norm(noise) - Gsample += noise - end - if target_space==:τ - Kfull = Kfunc( ir_omega_grid , target_fine_grid.grid, regular, omega0) - G_analy = SemiCircle(dlr, target_fine_grid.grid, target_space) - else - Kfull = Kfunc_freq( ir_omega_grid , target_fine_grid, regular, omega0) - G_analy = SemiCircle(dlr, target_fine_grid, target_space) - end - # if hasweight - # if space == :n - # Gsample = Gsample .* weight[] - # end - rho = K \ Gsample - G = Kfull * rho - - #interp_err = sqrt(Interp.integrate1D((G - G_analy) .^ 2, target_fine_grid )) - if target_space==:n - interp_err = norm(G - G_analy) - else - interp_err = sqrt(Interp.integrate1D( (abs.(G - G_analy)).^2, target_fine_grid )) - end - #print("condition KIR: $(cond(K))\n") - print("Exact Green err: $(interp_err)\n") - - -end -# function test_err(dlr, ir_grid, fine_grid, target_fine_grid, ir_omega_grid, space, regular, omega0) - -# #generate_grid_expan(eps, Lambda, expan_trunc, :ex, false, datatype(Lambda)) -# Gsample = SemiCircle(dlr, ir_grid, space) -# if space == :n -# K = Kfunc_freq(ir_omega_grid , ir_grid, regular, omega0) -# elseif space == :τ -# K = Kfunc(ir_omega_grid, ir_grid, regular, omega0) -# end -# Ktau = Kfunc( ir_omega_grid , target_fine_grid.grid, regular, omega0) -# rho = K \ Gsample -# G = Ktau * rho -# G_analy = SemiCircle(dlr, target_fine_grid, :τ) -# interp_err = sqrt(Interp.integrate1D((G - G_analy) .^ 2, target_fine_grid )) -# print("Exact Green err: $(interp_err)\n") - - -# end -# if abspath(PROGRAM_FILE) == @__FILE__ -# # dlr = DLRGrid(Euv=lambda, β=beta, isFermi=true, rtol=1e-12, symmetry=:sym) - -# datatype = Float64 -# #setprecision(128) -# #atatype = BigFloat -# isFermi = true -# symmetry = :none -# beta = datatype(1.0) -# Lambda = datatype(1000) -# eps = datatype(1e-10) -# n_trunc = datatype(10) #omega_n is truncated at n_trunc * Lambda -# expan_trunc = 1000 -# omega_grid, tau_grid, n_grid = generate_grid(eps, Lambda, n_trunc, :τ, false, datatype(Lambda)) - -# #generate_grid_resum(eps, Lambda, n_trunc, false, datatype(Lambda)) -# end diff --git a/test/svd.jl b/test/svd.jl deleted file mode 100644 index e79074b..0000000 --- a/test/svd.jl +++ /dev/null @@ -1,246 +0,0 @@ -include("kernel_svd.jl") - -function generate_grid(eps::T, Lambda::T, n_trunc::T, space::Symbol=:τ, regular = false, - omega0::T = Lambda, hasweight =false) where {T} - # generate frequency finegrid - w_grid = fine_ωGrid_test(T(Lambda), 24, T(1.5)) - weight_w = zeros(T,length(w_grid)) - #calculate grid weights - for i in 1:length(w_grid) - data = zeros(T,length(w_grid)) - data[i] = 1.0 - weight_w[i] = Interp.integrate1D(data, w_grid) - end - - #symmetrize the grid - wgrid = vcat(-w_grid.grid[end:-1:1], w_grid.grid) - weight_w = vcat(weight_w[end:-1:1], weight_w) - - - #generate tau fine grid - t_grid = fine_τGrid_test(T(Lambda),128, T(1.5)) - weight_t = zeros(T,length(t_grid)) - for i in 1:length(t_grid) - data = zeros(T,length(t_grid)) - data[i] = 1.0 - weight_t[i] = Interp.integrate1D(data, t_grid) - end - tgrid = t_grid.grid - - # generate fine n grid - - ngrid = nGrid_test(true, T(Lambda), 12, T(1.5)) - #ngrid = uni_ngrid(true, T(n_trunc*Lambda)) - fine_ngrid = uni_ngrid(true, T(n_trunc*Lambda)) - omega = (2*ngrid.+1)* π - dlr = DLRGrid(Lambda, beta, eps, true, :none, dtype=T) - - Kn = Kfunc_freq(wgrid, Int.(ngrid), weight_w, regular, omega0) - if hasweight - Ktau = Kfunc(wgrid, tgrid, weight_w, weight_t, regular, omega0) - else - Ktau = Kfunc(wgrid, tgrid, regular, omega0) - end - - left = searchsortedfirst(omega, -n_trunc*Lambda/10) - right = searchsortedfirst(omega, n_trunc*Lambda/10) - - if space == :n - eig = svd(Kn, full = true) - elseif space == :τ - eig = svd(Ktau, full = true) - end - - idx = searchsortedfirst(eig.S./eig.S[1], eps, rev=true) - - print("rank: $(idx)\n") - if space == :n - pivoted_idx, n_grid = IR(ngrid, eig.U, idx, "omega_n") - Un_full = eig.U - n_idx = pivoted_idx - #print("tail selected: left $(ngrid[left] in n_grid) right $(ngrid[right] in n_grid)\n") - elseif space == :τ - pivoted_idx, tau_grid = IR(tgrid, eig.U, idx, "tau", Lambda) - Utau_full = eig.U - tau_idx = pivoted_idx - end - - omega_idx, omega_grid = IR(wgrid, eig.V, idx, "omega") - - #Use implicit fourier to get U in the other space - if space == :n - U = (Ktau * eig.V) - elseif space == :τ - U = (Kn * eig.V) - #U_compare = F*eig.U - end - - for i in 1:idx - U[:, i] = U[:, i] ./ eig.S[i] - end - - if space == :n - pivoted_idx, tau_grid = IR(tgrid, U, idx, "tau") - Utau_full = U - tau_idx = pivoted_idx - elseif space == :τ - pivoted_idx, n_grid = IR(ngrid, U, idx, "omega_n", Lambda) - Un_full = U - n_idx = pivoted_idx - end - - maxidx = searchsortedfirst(eig.S./eig.S[1], 1e-16, rev=true) - fine_n_idx = zeros(Int, length(n_grid)) - fidx = 1 - for i in eachindex(n_grid) - while Int(n_grid[i]) != Int(fine_ngrid[fidx]) - fidx += 1 - end - fine_n_idx[i] = fidx - end - print("test idx: $(fine_ngrid[fine_n_idx]) $(n_grid)\n") - Kn_fine = Kfunc_freq(wgrid, Int.(fine_ngrid), weight_w, regular, omega0) - # This test with U matrix - # n space sparse grid: n_grid - # n space fine grid: ngrid - # τ space sparse grid: tau_grid - # τ space fine grid: tgrid - - #test_err(dlr, tau_grid, tgrid, tgrid, :τ, :τ, Un_full[:, 1:maxidx], n_idx, Utau_full[:, 1:maxidx], tau_idx, collect(1:maxidx), idx; - #test_err(dlr, Int.(n_grid), Int.(ngrid), tgrid, :n, :τ, Un_full[:, 1:maxidx], n_idx, Utau_full[:, 1:maxidx], tau_idx, collect(1:maxidx), idx; - - # This test with K matrix - #test_err(dlr, Int.(n_grid), Int.(ngrid), t_grid, :n, :τ, Kn, n_idx, Ktau, tau_idx, omega_idx, idx; - test_err(dlr, Int.(n_grid), Int.(fine_ngrid), t_grid, :n, :τ, Kn_fine, fine_n_idx, Ktau, tau_idx, omega_idx, idx; - case = "SemiCircle", hasnoise = true, hasweight=hasweight, weight_tau = sqrt.(weight_t)) - - filename = "newDLReig.txt" - folder="./" - file = open(joinpath(folder, filename), "a") - #max_res = maximum((res[:])) - for i in 1:idx - @printf(file, "%32.30g\n", eig.S[i]) - end - close(file) - - return omega_grid, tau_grid, n_grid -end - - - -function test_err(dlr, ir_grid, fine_grid, target_fine_grid, space, target_space, Un_full, n_idx, Utau_full, tau_idx, omega_idx, rank; - case = "SemiCircle", hasnoise=false, hasweight=false, weight_tau =nothing) - print("size: $(size(Un_full)) $(size(Utau_full))\n") - if space == :n - U11 = Un_full[sort(n_idx[1:rank]), sort(omega_idx[1:rank])] - U12 = Un_full[sort(n_idx[1:rank]), sort(omega_idx[rank+1:end])] - Uinit_full = Un_full - elseif space== :τ - U11 = Utau_full[sort(tau_idx[1:rank]), sort(omega_idx[1:rank])] - U12 = Utau_full[sort(tau_idx[1:rank]), sort(omega_idx[rank+1:end])] - Uinit_full = Utau_full - end - if target_space == :τ - U21 = Utau_full[sort(tau_idx[rank+1 : end]),sort(omega_idx[1:rank])] - U22 = Utau_full[sort(tau_idx[rank+1 : end]),sort(omega_idx[rank+1:end])] - U_full = Utau_full - elseif target_space== :n - U21 = Un_full[sort(n_idx[rank+1 : end]), sort(omega_idx[1:rank])] - U22 = Un_full[sort(n_idx[rank+1 : end]), sort(omega_idx[rank+1:end])] - U_full = Un_full - end - print("noise err bound: $(opnorm(U21*inv(U11)))\n") - print("Un * inv(U11): $(opnorm(Un_full[:, sort(omega_idx[1:rank])]*inv(U11)))\n") - print("Utau * inv(U11): $(opnorm(Utau_full[:, sort(omega_idx[1:rank])]*inv(U11)))\n") - print("epsilon err bound: $(opnorm(U21*inv(U11)*U12 -U22))\n") - #generate_grid_expan(eps, Lambda, expan_trunc, :ex, false, datatype(Lambda)) - - N = 1 - if case == "MultiPole" - N = 5 - N_poles = 100 - dtype = typeof(Utau_full[1,1]) - poles = zeros(dtype, (N, N_poles)) - weights = zeros(dtype, (N, N_poles)) - Random.seed!(8) - - for i in 1:N - #poles[i, :] = dlr.ω - poles[i, :] = 2.0 * rand(dtype, N_poles) .- 1.0 - weights[i, :] = rand(dtype, N_poles)#2.0 * rand(dtype, N_poles) .- 1.0 - weights[i, :] = weights[i, :] / sum(abs.(weights[i, :])) - end - end - - for i in 1:N - if case == "SemiCircle" - Gsample = SemiCircle(dlr, ir_grid, space) - Gsample_full = SemiCircle(dlr, fine_grid, space) - G_analy = SemiCircle(dlr, target_fine_grid, target_space) - else - Gsample = MultiPole(dlr, ir_grid, space, poles[i], weights[i]) - Gsample_full = MultiPole(dlr, fine_grid, space, poles[i], weights[i]) - G_analy = MultiPole(dlr, target_fine_grid, target_space, poles[i], weights[i]) - end - - if hasnoise - T = typeof(Gsample[1]) - noise = (rand(T, length(Gsample))) - noise = dlr.rtol * noise ./ norm(noise) - print("err norm: $(norm(noise))\n") - Gsample += noise - end - - - if hasweight - if target_space == :τ - G_analy .*= weight_tau - end - if space == :τ - Gsample .*= weight_tau[sort(tau_idx[1:rank])] - Gsample_full .*= weight_tau - end - end - #Gsample += 1e-11 * (randn(length(Gsample)) + im * randn(length(Gsample))) - #print("U21U11^-1: $(norm(Ures/U_IR))\n") - rho = U11 \ Gsample - G = U_full[:, sort(omega_idx[1:rank])] * rho - GIR = U11 *rho - - print("norm C_IR: $(norm(rho))\n") - init_rho_full = Uinit_full \ Gsample_full - print("norm in $(space): C1: $(norm(init_rho_full[ sort(omega_idx[1:rank])])) C2:$(norm(init_rho_full[ sort(omega_idx[rank+1:end])]))\n") - rho_full = U_full \ G_analy - print("norm in $(target_space): C1: $(norm(rho_full[ sort(omega_idx[1:rank])])) C2:$(norm(rho_full[ sort(omega_idx[rank+1:end])]))\n") - if target_space == :τ - interp_err = (norm(G - G_analy)) - #interp_err = sqrt(Interp.integrate1D((G - G_analy) .^ 2, target_fine_grid )) - else - interp_err = (norm(G - G_analy)) - end - print("IR err: $(maximum(abs.(GIR - Gsample)))\n") - print("full err: $(maximum(abs.(G - G_analy)))\n") - - print("Exact Green err: $(interp_err/dlr.rtol)\n") - end - - end - - -if abspath(PROGRAM_FILE) == @__FILE__ - # dlr = DLRGrid(Euv=lambda, β=beta, isFermi=true, rtol=1e-12, symmetry=:sym) - - datatype = Float64 - #setprecision(128) - #atatype = BigFloat - isFermi = true - symmetry = :none - beta = datatype(1.0) - Lambda = datatype(100000) - eps = datatype(1e-6) - n_trunc = datatype(10) #omega_n is truncated at n_trunc * Lambda - expan_trunc = 100 - omega_grid, tau_grid, n_grid = generate_grid(eps, Lambda, n_trunc, :τ, false, datatype(Lambda), true) - #omega_grid, tau_grid, n_grid = generate_grid_expan(eps, Lambda, expan_trunc, :τ, false, datatype(Lambda)) - #generate_grid_resum(eps, Lambda, n_trunc, false, datatype(Lambda)) -end diff --git a/test/svd_sym.jl b/test/svd_sym.jl deleted file mode 100644 index 7ca06eb..0000000 --- a/test/svd_sym.jl +++ /dev/null @@ -1,499 +0,0 @@ -include("kernel_svd_sym.jl") -include("../src/vertex3/QR.jl") -include("../src/vertex3/omega.jl") -include("../src/vertex3/matsubara.jl") -using DoubleFloats - - -function IR_new(space, U, finegrid, sym::T, Lambda::T, eps::T) where {T} - _, idx = size(U) - if space ==:τ - finemesh = TauFineMesh(U, finegrid; sym=sym) - basis = FQR.Basis{TauGrid, T, T}(Lambda, eps, finemesh) - elseif space == :n - finemesh = MatsuFineMesh{T}(U, finegrid, true; sym=sym) - basis = FQR.Basis{MatsuGrid, T, Complex{T}}(Lambda, eps, finemesh) - elseif space==:ω - finemesh = TauFineMesh(U, finegrid; sym=sym) - basis = FQR.Basis{TauGrid, T, T}(Lambda, eps, finemesh) - end - - idxlist = [] - while length(idxlist) < idx - selected = findall(basis.mesh.selected .== false) - maxResidual, idx_inner = findmax(basis.mesh.residual[selected]) - idx_inner = selected[idx_inner] - FQR.addBasisBlock!(basis, idx_inner, false) - append!(idxlist, Int(idx_inner)) - if sym>0 - idx_mirror = length(finegrid) - idx_inner +1 - #print("idxmirror $(idx_mirror)\n") - append!(idxlist, Int(idx_mirror)) - end - end - return sort(idxlist), finegrid[sort(idxlist)] -end -function generate_grid(eps::T, Lambda::T, n_trunc::T, space::Symbol=:τ, regular = false, - omega0::T = Lambda, hasweight =false) where {T} - sym = T(1.0) - # generate frequency finegrid - w_grid = fine_ωGrid_test(T(Lambda), 24, T(1.5)) - weight_w = zeros(T,length(w_grid)) - #calculate grid weights - for i in 1:length(w_grid) - data = zeros(T,length(w_grid)) - data[i] = 1.0 - weight_w[i] = Interp.integrate1D(data, w_grid) - end - - #symmetrize the grid - wgrid = vcat(-w_grid.grid[end:-1:1], w_grid.grid) - weight_w = vcat(weight_w[end:-1:1], weight_w) - - #generate tau fine grid - t_grid = fine_τGrid_test(T(Lambda),128, T(1.5)) - weight_t = zeros(T,length(t_grid)) - for i in 1:length(t_grid) - data = zeros(T,length(t_grid)) - data[i] = 1.0 - weight_t[i] = Interp.integrate1D(data, t_grid) - end - tgrid = t_grid.grid - - # generate fine n grid - - ngrid = nGrid_test(true, T(Lambda), 12, T(1.5)) - #ngrid = Int.(uni_ngrid(true, T(n_trunc*Lambda))) - fine_ngrid = Int.(uni_ngrid(true, T(n_trunc*Lambda))) - omega = (2*ngrid.+1)* π - #dlr = DLRGrid(Lambda, beta, eps, true, :none, dtype=T) - dlr = DLRGrid(Lambda, beta, eps, true, :none, dtype=T) - Kn = Kfunc_freq(wgrid, Int.(ngrid), weight_w, regular, omega0) - if hasweight - Ktau = Kfunc(wgrid, tgrid, weight_w, weight_t, regular, omega0) - else - Ktau = Kfunc(wgrid, tgrid, regular, omega0) - end - - left = searchsortedfirst(omega, -n_trunc*Lambda/10) - right = searchsortedfirst(omega, n_trunc*Lambda/10) - - if space == :n - eig = svd(Kn, full = true) - elseif space == :τ - eig = svd(Ktau, full = true) - end - - #idx = searchsortedfirst(eig.S./eig.S[1], eps, rev=true) - # if idx0 - idx_mirror = length(finegrid) - idx_inner +1 - #print("idxmirror $(idx_mirror)\n") - append!(idxlist, Int(idx_mirror)) - end - end - selected = findall(basis.mesh.selected .== false) - return vcat(sort(idxlist), selected), finegrid[sort(idxlist)] -end - -function generate_grid(eps::T, Lambda::T, n_trunc::T, space::Symbol=:τ, regular = false, - omega0::T = Lambda, hasweight =false) where {T} - sym = T(1.0) - # generate frequency finegrid - w_grid = fine_ωGrid_test(T(Lambda), 24, T(1.5)) - weight_w = zeros(T,length(w_grid)) - #calculate grid weights - for i in 1:length(w_grid) - data = zeros(T,length(w_grid)) - data[i] = 1.0 - weight_w[i] = Interp.integrate1D(data, w_grid) - end - - #symmetrize the grid - wgrid = vcat(-w_grid.grid[end:-1:1], w_grid.grid) - weight_w = vcat(weight_w[end:-1:1], weight_w) - - #generate tau fine grid - t_grid = fine_τGrid_test(T(Lambda),128, T(1.5)) - weight_t = zeros(T,length(t_grid)) - for i in 1:length(t_grid) - data = zeros(T,length(t_grid)) - data[i] = 1.0 - weight_t[i] = Interp.integrate1D(data, t_grid) - end - tgrid = t_grid.grid - - # generate fine n grid - - ngrid = nGrid_test(true, T(Lambda), 12, T(1.5)) - #ngrid = Int.(uni_ngrid(true, T(n_trunc*Lambda))) - fine_ngrid = Int.(uni_ngrid(true, T(n_trunc*Lambda))) - omega = (2*ngrid.+1)* π - #dlr = DLRGrid(Lambda, beta, eps, true, :none, dtype=T) - dlr = DLRGrid(Lambda, beta, eps, true, :none, dtype=T) - Kn = Kfunc_freq(wgrid, Int.(ngrid), sqrt.(weight_w)) - Ktau = Kfunc(wgrid, tgrid, sqrt.(weight_w), sqrt.(weight_t)) - - left = searchsortedfirst(omega, -n_trunc*Lambda/10) - right = searchsortedfirst(omega, n_trunc*Lambda/10) - - if space == :n - eig = svd(Kn, full = true) - elseif space == :τ - eig = svd(Ktau, full = true) - end - - idx = searchsortedfirst(eig.S./eig.S[1], eps, rev=true) - # if idx