diff --git a/src/vertex3/QR.jl b/src/vertex3/QR.jl index eda846c..7b9c1ce 100644 --- a/src/vertex3/QR.jl +++ b/src/vertex3/QR.jl @@ -101,10 +101,10 @@ function addBasisBlock!(basis::Basis, idx, verbose) basis.mesh.residual[idx] = 0 # the selected mesh grid has zero residual #print("$(mirror(basis.mesh, idx))\n") gridmirror , idxmirror = mirror(basis.mesh, idx) - print(idxmirror) + #print("mirror $(idx) $(idxmirror)\n") for (gi,grid) in enumerate(gridmirror) addBasis!(basis, grid, verbose) - L2Res = L2Residual(basis.mesh) + #L2Res = L2Residual(basis.mesh) # if typeof(basis.mesh)<:FreqFineMesh # if output && num % 5 == 0 && typeof(basis.mesh)<: TauFineMesh#MatsuFineMesh # pic = plot(ylabel = "residual") @@ -168,24 +168,24 @@ function updateResidual!(basis::Basis{Grid, Mesh, F, D}) where {Grid,Mesh,F,D} end end end - if basis.mesh.simplegrid - Threads.@threads for idx in 1:length(mesh.L2grids) - candidate = mesh.L2grids[idx] - pp = sum(q[j] * dot(mesh, candidate, basis.grid[j]) for j in 1:basis.N) - #pp = q[basis.N] * dot(mesh, candidate, basis.grid[basis.N]) - _residual = mesh.residual_L2[idx] - abs(pp) * abs(pp) - # @assert isnan(_residual) == false "$pp and $([q[j] for j in 1:basis.N]) => $([dot(mesh, basis.grid[j], candidate) for j in 1:basis.N])" - # println("working on $candidate : $_residual") - if _residual < 0 - if _residual < -basis.rtol - @warn("warning: residual smaller than 0 at $candidate got $(mesh.residual_L2[idx]) - $(abs(pp)^2) = $_residual") - end - mesh.residual_L2[idx] = 0 - else - mesh.residual_L2[idx] = _residual - end - end - end + # if basis.mesh.simplegrid + # Threads.@threads for idx in 1:length(mesh.L2grids) + # candidate = mesh.L2grids[idx] + # pp = sum(q[j] * dot(mesh, candidate, basis.grid[j]) for j in 1:basis.N) + # #pp = q[basis.N] * dot(mesh, candidate, basis.grid[basis.N]) + # _residual = mesh.residual_L2[idx] - abs(pp) * abs(pp) + # # @assert isnan(_residual) == false "$pp and $([q[j] for j in 1:basis.N]) => $([dot(mesh, basis.grid[j], candidate) for j in 1:basis.N])" + # # println("working on $candidate : $_residual") + # if _residual < 0 + # if _residual < -basis.rtol + # @warn("warning: residual smaller than 0 at $candidate got $(mesh.residual_L2[idx]) - $(abs(pp)^2) = $_residual") + # end + # mesh.residual_L2[idx] = 0 + # else + # mesh.residual_L2[idx] = _residual + # end + # end + # end end # function updateResidual!(basis::Basis{Grid, Mesh, F, D}, candidate) where {Grid,Mesh,F,D} @@ -227,9 +227,9 @@ function GramSchmidt(basis::Basis{G,M,F,D}) where {G,M,F,D} newgrid = basis.grid[end] overlap = [dot(basis.mesh, basis.grid[j], newgrid) for j in 1:basis.N-1] - if !isempty(overlap) - println( "$(maximum(imag(overlap)))\n" ) - end + # if !isempty(overlap) + # println( "$(maximum(imag(overlap)))\n" ) + # end for qi in 1:basis.N-1 _R[qi, end] = basis.Q[:, qi]' * overlap _Q[:, end] -= _R[qi, end] * _Q[:, qi] # q @@ -238,7 +238,7 @@ function GramSchmidt(basis::Basis{G,M,F,D}) where {G,M,F,D} _norm = dot(basis.mesh, newgrid, newgrid) - _R[:, end]' * _R[:, end] _norm = sqrt(abs(_norm)) - @assert _norm > eps(F(1)) * 100 "$_norm is too small as a denominator!\nnewgrid = $newgrid\nexisting grid = $(basis.grid)\noverlap=$overlap\nR=$_R\nQ=$_Q" + @assert _norm > eps(F(1)) * 100 "$_norm is too small as a denominator!"#\nnewgrid = $newgrid\nexisting grid = $(basis.grid)\noverlap=$overlap\nR=$_R\nQ=$_Q" _R[end, end] = _norm _Q[:, end] /= _norm diff --git a/src/vertex3/frequency.jl b/src/vertex3/frequency.jl index 83b3c49..54f3e00 100644 --- a/src/vertex3/frequency.jl +++ b/src/vertex3/frequency.jl @@ -67,8 +67,8 @@ struct FreqFineMesh{D,Float,Double} <: FQR.FineMesh # elseif D == 3 elseif D == 1 if simplegrid - candidate_grid = log_ωGrid(Float(init), Float(factor*Λ), Float(ratio))# Float(1.35^(log(1e-6) / log(rtol)))) - #candidate_grid = matsu_ωGrid(50, Float(1.0)) + #candidate_grid = log_ωGrid(Float(init), Float(factor*Λ), Float(ratio))# Float(1.35^(log(1e-6) / log(rtol)))) + candidate_grid = matsu_ωGrid(80, Float(1.0)) #fine_ωGrid(Float(10Λ), 1, Float(1.3)) else candidate_grid = _finegrid diff --git a/src/vertex3/generate_basis.jl b/src/vertex3/generate_basis.jl index a96f1e6..54323f8 100644 --- a/src/vertex3/generate_basis.jl +++ b/src/vertex3/generate_basis.jl @@ -260,12 +260,12 @@ end if abspath(PROGRAM_FILE) == @__FILE__ D = 1 - Err = [-6, ] + Err = [-10, ] #Λ = [1e4, ] #Err = [-6, -8, -10, -12, -14] #Err = [-5, -7, -9, -11, -13] - Λ = [1e2] # [1e3, 1e4, 1e5, 1e6,1e7] - setprecision(128) + Λ = [1e4] # [1e3, 1e4, 1e5, 1e6,1e7] + setprecision(512) Float = BigFloat Double = BigFloat #Float = Float64 @@ -299,7 +299,7 @@ if abspath(PROGRAM_FILE) == @__FILE__ # for rat in ratio_scan # for init in init_scan #mesh = FreqFineMesh{D,Float, Double}(lambda, rtol, sym=sym, degree = degree1, ratio = ratio1, simplegrid=true) - mesh = FreqFineMesh{D,Float, Double}(lambda, rtol, sym=sym, degree = degree1, ratio = ratio1, factor = 1000, init=init, simplegrid=false) + mesh = FreqFineMesh{D,Float, Double}(lambda, rtol, sym=sym, degree = degree1, ratio = ratio1, factor = 1000, init=init, simplegrid=true) basis = FQR.Basis{FreqGrid, Float ,Double}(lambda, rtol, mesh) L2Res = L2Residual(basis.mesh) diff --git a/src/vertex3/matsubara.jl b/src/vertex3/matsubara.jl index f286e78..15b9204 100644 --- a/src/vertex3/matsubara.jl +++ b/src/vertex3/matsubara.jl @@ -29,7 +29,7 @@ struct MatsuFineMesh{Float} <: FQR.FineMesh _finegrid = Int.(nGrid(isFermi, Float(Λ), degree, Float(ratio) )) #print("ngrid $(_finegrid)\n") # separationTest(_finegrid) - mesh = new(isFermi, sym, [], [], [], _finegrid) + mesh = new{Float}(isFermi, sym, [], [], [], _finegrid) for (xi, x) in enumerate(_finegrid) coord = xi @@ -48,6 +48,22 @@ struct MatsuFineMesh{Float} <: FQR.FineMesh return mesh end + + function MatsuFineMesh{Float}(U::AbstractMatrix, finegrid, isFermi; sym=1) where {Float} + mesh = new{Float}(isFermi, sym, [], [], [], finegrid) + + for (xi, x) in enumerate(finegrid) + coord = xi + #if irreducible(coord, sym, Nfine) # if grid point is in the reducible zone, then skip residual initalization + vec = U[xi, :] + g = MatsuGrid(x, coord, vec) + push!(mesh.candidates, g) + push!(mesh.residual, real.(FQR.dot(mesh, g, g)) ) + push!(mesh.selected, false) + #end + end + return mesh + end end function Freq2Index(isFermi, ωnList) diff --git a/src/vertex3/omega.jl b/src/vertex3/omega.jl new file mode 100644 index 0000000..d715d39 --- /dev/null +++ b/src/vertex3/omega.jl @@ -0,0 +1,244 @@ +using Lehmann +using StaticArrays, Printf +using CompositeGrids +using LinearAlgebra +using DelimitedFiles +#const Float = BigFloat #FQR.Float +#const Double =BigFloat #FQR.Double +#const DotF = BigFloat +#const Tiny = DotF(1e-5) + +struct TauGrid{Float} <: FQR.Grid + tau::Float # actual location of the grid point + coord::Int # integer coordinate of the grid point on the fine meshes + vec::Vector{Float} +end + +Base.show(io::IO, grid::TauGrid) = print(io, "τ = ($(@sprintf("%12.4f", grid.tau[1])))") + +struct TauFineMesh{Float} <: FQR.FineMesh + symmetry::Int # symmetrize (omega1, omega2) <-> (omega2, omega1) + candidates::Vector{TauGrid{Float}} # vector of grid points + selected::Vector{Bool} + residual::Vector{Float} + + ## for frequency mesh only ### + fineGrid::Vector{Float} # fine grid for each dimension + + function TauFineMesh(U::AbstractMatrix{Float}, finegrid; sym=1) where {Float} + # initialize the residual on fineGrid with + #_finegrid = Float.(τChebyGrid(Λ)) + # separationTest(_finegrid) + mesh = new{Float}(sym, [], [], [], finegrid) + + for (xi, x) in enumerate(finegrid) + coord = xi + #if irreducible(coord, sym, Nfine) # if grid point is in the reducible zone, then skip residual initalization + vec = U[xi,:] + g = TauGrid(x, coord, vec) + push!(mesh.candidates, g) + push!(mesh.residual, FQR.dot(mesh, g, g)) + push!(mesh.selected, false) + #end + end + + println("fine mesh initialized.") + return mesh + end +end + + + +# function τChebyGrid(Λ, degree=24, print = true) +# npt = Int(ceil(log(Λ) / log(2.0))) - 2 # subintervals on [0,1/2] in tau space (# subintervals on [0,1] is 2*npt) + +# pbpt = zeros(Float, 2npt + 1) +# pbpt[1] = 0.0 +# for i = 1:npt +# pbpt[i+1] = 1.0 / 2^(npt - i + 1) +# end +# pbpt[npt+2:2npt+1] = 1 .- pbpt[npt:-1:1] +# #println("fine grid size: $(length(grid)) within [$(grid[1]), $(grid[end])]") +# finegrid = Lehmann.Discrete.CompositeChebyshevGrid(degree, pbpt).grid +# #println("$(finegrid)\n")#[1:(length(finegrid)÷2+1)]) +# #println("$(finegrid+reverse(finegrid))\n") +# return finegrid +# end + +""" +composite expoential grid +""" +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( + :cheb,# 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 + ) + #print(grid[1:length(grid)÷2+1]) + #print(grid+reverse(grid)) + # println("Composite expoential grid size: $(length(grid))") + println("fine grid size: $(length(grid)) within [$(grid[1]), $(grid[end])]") + return grid + + ############# DLR based fine grid ########################################## + # dlr = DLRGrid(Euv=Float64(Λ), beta=1.0, rtol=Float64(rtol) / 100, isFermi=true, symmetry=:ph, rebuild=true) + # # println("fine basis number: $(dlr.size)\n", dlr.ω) + # degree = 4 + # grid = Vector{Double}(undef, 0) + # panel = Double.(dlr.τ) + # for i in 1:length(panel)-1 + # uniform = [panel[i] + (panel[i+1] - panel[i]) / degree * j for j in 0:degree-1] + # append!(grid, uniform) + # end + + # println("fine grid size: $(length(grid)) within [$(grid[1]), $(grid[2])]") + # return grid +end + +# """ +# Test the finegrids do not overlap +# """ +# function separationTest(finegrid) +# epsilon = eps(DotF(1)) * 10 +# for (i, f) in enumerate(finegrid) +# # either zero, or sufficiently large +# @assert abs(f) < epsilon || abs(f) > Tiny "$i: $f should either smaller than $epsilon or larger than $Tiny" +# for (j, g) in enumerate(finegrid) +# # two frequencies are either the same, or well separated +# @assert abs(f - g) < epsilon || abs(f - g) > Tiny "$i: $f and $j: $g should either closer than $epsilon or further than $Tiny" +# fg = f + g +# for (k, l) in enumerate(finegrid) +# @assert abs(l - fg) < epsilon || abs(l - fg) > Tiny "$i: $f + $j: $g = $fg and $k: $l should either closer than $epsilon or further than $Tiny" +# end +# end +# end +# return +# end + +# function irreducible(coord, symmetry,length) +# @assert iseven(length) "The fineGrid should have even number of points" +# if symmetry == 0 +# return true +# else +# return coord ", length(mirror(mesh, i))) + + # for grid in FQR.mirror(mesh, i) + # if grid.sector == 1 + # xp, yp = grid.coord + # residual[xp, yp] = residual[x, y] + # # println(xp, ", ", yp) + # end + # end + # end + # end + + # for i in 1:Nfine + # for j in 1:Nfine + # println(io, residual[i, j]) + # end + # end + # end +end diff --git a/src/vertex3/tau.jl b/src/vertex3/tau.jl index cdc4cfc..329fbf5 100644 --- a/src/vertex3/tau.jl +++ b/src/vertex3/tau.jl @@ -46,11 +46,33 @@ struct TauFineMesh{Float} <: FQR.FineMesh #end end + println("fine mesh initialized.") + return mesh + end + function TauFineMesh(U::AbstractMatrix{Float}, finegrid; sym=1) where {Float} + # initialize the residual on fineGrid with + #_finegrid = Float.(τChebyGrid(Λ)) + # separationTest(_finegrid) + mesh = new{Float}(sym, [], [], [], finegrid) + + for (xi, x) in enumerate(finegrid) + coord = xi + #if irreducible(coord, sym, Nfine) # if grid point is in the reducible zone, then skip residual initalization + vec = U[xi,:] + g = TauGrid(x, coord, vec) + push!(mesh.candidates, g) + push!(mesh.residual, FQR.dot(mesh, g, g)) + push!(mesh.selected, false) + #end + end + println("fine mesh initialized.") return mesh end end + + # function τChebyGrid(Λ, degree=24, print = true) # npt = Int(ceil(log(Λ) / log(2.0))) - 2 # subintervals on [0,1/2] in tau space (# subintervals on [0,1] is 2*npt) diff --git a/test/data_generate.jl b/test/data_generate.jl index 95bcb4c..8a84927 100644 --- a/test/data_generate.jl +++ b/test/data_generate.jl @@ -135,7 +135,7 @@ function L2normτ(value_dlr, dlr, case, grid, poles=nothing, weight=nothing, spa # println("fine grid size: $(length(grid)) within [$(grid[1]), $(grid[2])]") # return grid - fineGrid = fine_τGrid(dlr.Euv, 12, typeof(dlr.Euv)(1.5), :gauss) + fineGrid = fine_τGrid(dlr.Euv, 128, typeof(dlr.Euv)(1.5), :gauss) if space == :n value = real(matfreq2tau(dlr, value_dlr, fineGrid, grid)) else @@ -320,7 +320,7 @@ function test_err_cheby(case, isFermi, symmetry, Euv, β, eps, poles, weights; d eta2 = 0.0 block = zeros(dtype, 10) print(dtype) - omega_grid, tau_grid, n_grid= generate_grid(dlr10.rtol, dlr10.Euv, dtype(10.0), :t, false, dlr10.Euv) + omega_grid, tau_grid, n_grid= generate_grid(dlr10.rtol, dlr10.Euv, dtype(10.0), :τ, false, dlr10.Euv) #e3 e-8 3 1.9 / e5 e-8 5 2.3 110 86/ e5 -12 7 2.3 140 118/ degree = 3 ratio = dtype(2.5^(log(1e-4) / log(eps))) @@ -352,6 +352,8 @@ function test_err_cheby(case, isFermi, symmetry, Euv, β, eps, poles, weights; d # println("ngrid: $(length(sample_ngrid)) $(length(dlr.n))") # println("taugrid: $(length(sample_taugrid)) $(length(dlr.τ))") if case == MultiPole + maxerr= [] + maxerr2 = [ ] for i in 1:N eff_poles = poles[i, :] weight = weights[i, :] @@ -371,8 +373,9 @@ function test_err_cheby(case, isFermi, symmetry, Euv, β, eps, poles, weights; d value2, err2, max_err2 = L2normτ(Gndlr, dlr10, case, dlr10.n, eff_poles, weight, :n) #modulus = abs(sum(dlreff)) - - println("eta: $(err/eps) $(max_err) $(err2/eps) $(max_err2)\n") + append!(maxerr, err) + append!(maxerr2, err2) + #println("eta: $(err/eps) $(max_err) $(err2/eps) $(max_err2)\n") if output @@ -396,6 +399,7 @@ function test_err_cheby(case, isFermi, symmetry, Euv, β, eps, poles, weights; d max_err_sum += max_err block[(i-1)÷(N÷10)+1] += err / N * 10 end + println("max: $(maximum(maxerr)/eps) $(maximum(maxerr2)/eps)\n") else Gtaudlr = case(dlr10, dlr10.τ, :τ) @@ -440,11 +444,11 @@ function test_err_cheby(case, isFermi, symmetry, Euv, β, eps, poles, weights; d end -#cases = [MultiPole] -cases = [SemiCircle] -Λ = [1.6e3] #, 1e4, 1e5, 1e6, 1e7] +cases = [MultiPole] +#cases = [SemiCircle] +Λ = [1e4] #, 1e4, 1e5, 1e6, 1e7] -rtol = [1e-6,]# 1e-12] +rtol = [1e-8,]# 1e-12] #rtol = [1e-4, 1e-6, 1e-8, 1e-10, 1e-12]# 1e-6 , 1e-8, 1e-10, 1e-12] #Λ = [1e3, 1e5, 1e6, 1e7] @@ -452,8 +456,8 @@ rtol = [1e-6,]# 1e-12] #rtol = [1e-10] -N_poles = 100 -N = 100 +N_poles = 2 +N = 10000 dtype = Float64 #dtype = BigFloat #setprecision(128) @@ -472,7 +476,8 @@ for i in 1:N weights[i, :] = rand(dtype, N_poles)#2.0 * rand(dtype, N_poles) .- 1.0 weights[i, :] = weights[i, :] / sum(abs.(weights[i, :])) end - +print("poles $(poles[1:10,1])\n") +print("weights $(weights[1:10,1])\n") for case in cases for l in Λ for r in rtol diff --git a/test/grid.jl b/test/grid.jl index 02b5a46..837abe6 100644 --- a/test/grid.jl +++ b/test/grid.jl @@ -3,7 +3,7 @@ function fine_ωGrid_test(Λ::Float, degree, ratio::Float) where {Float} N = Int(floor(log(Λ) / log(ratio) + 1)) grid = CompositeGrid.LogDensedGrid( - :cheb,# The top layer grid is :gauss, optimized for integration. For interpolation use :cheb + :gauss,# 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 diff --git a/test/kernel_svd_sym.jl b/test/kernel_svd_sym.jl new file mode 100644 index 0000000..a63073d --- /dev/null +++ b/test/kernel_svd_sym.jl @@ -0,0 +1,442 @@ +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 new file mode 100644 index 0000000..e161075 --- /dev/null +++ b/test/kernel_svd_wh.jl @@ -0,0 +1,528 @@ +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_sym.jl b/test/svd_sym.jl new file mode 100644 index 0000000..7ca06eb --- /dev/null +++ b/test/svd_sym.jl @@ -0,0 +1,499 @@ +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