diff --git a/EC.jl b/EC.jl index 372cc2f..30647a1 100644 --- a/EC.jl +++ b/EC.jl @@ -3,8 +3,8 @@ include("helper.jl") "EC model for a Hamiltonian family H(c) = H0 + c * H1" mutable struct affine_EC - H0::SparseMatrixCSC{ComplexF64} - H1::SparseMatrixCSC{ComplexF64} + H0::AbstractMatrix{ComplexF64} + H1::AbstractMatrix{ComplexF64} weights::Vector{ComplexF64} trained::Bool @@ -16,7 +16,7 @@ mutable struct affine_EC exact_E::Vector{ComplexF64} extrapolated_E::Vector{ComplexF64} - affine_EC(H0::SparseMatrixCSC{ComplexF64}, H1::SparseMatrixCSC{ComplexF64}, weights::Vector{ComplexF64}=ones(ComplexF64, size(H0, 1))) = new(H0, H1, weights, false, nothing, nothing, nothing, ComplexF64[], ComplexF64[], ComplexF64[]) + affine_EC(H0::AbstractMatrix{ComplexF64}, H1::AbstractMatrix{ComplexF64}, weights::Vector{ComplexF64}=ones(ComplexF64, size(H0, 1))) = new(H0, H1, weights, false, nothing, nothing, nothing, ComplexF64[], ComplexF64[], ComplexF64[]) end "Train an EC model for a given range of c values. @@ -67,7 +67,7 @@ end "Extrapolate using a trained EC model for a given range of c values If a list is provided for ref_eval, they are used as reference values for picking the closest eigenvalues at each point. If a single number is provided for ref_eval, it is used as a reference for the first point, and the previous eigenvalue is used as the reference for each successive point." -function extrapolate!(EC::affine_EC, c_vals; ref_eval=EC.training_E[end], verbose=true, tol=1e-5) +function extrapolate!(EC::affine_EC, c_vals; ref_eval=EC.training_E[end], verbose=true, tol=1e-5, precalculated_exact_E=nothing) @assert EC.trained "EC model must be trained using train() before extrapolation" for c in c_vals @@ -81,10 +81,14 @@ function extrapolate!(EC::affine_EC, c_vals; ref_eval=EC.training_E[end], verbos current_E = popfirst!(ref_eval) end - H = EC.H0 + c .* EC.H1 - - evals, _ = eigs(H, sigma=current_E, maxiter=5000, tol=tol, ritzvec=false, check=1) - current_E = nearest(evals, current_E) + if isnothing(precalculated_exact_E) + verbose && println("Extact solution for c = $c") + H = EC.H0 + c .* EC.H1 + evals, _ = eigs(H, sigma=current_E, maxiter=5000, tol=tol, ritzvec=false, check=1) + current_E = nearest(evals, current_E) + else + current_E = popfirst!(precalculated_exact_E) + end push!(EC.exact_E, current_E) @@ -101,9 +105,16 @@ end exportCSV(EC::affine_EC, filename) = exportCSV(filename, (EC.training_E, EC.exact_E, EC.extrapolated_E), ("training", "exact", "extrapolated")) "Plot EC data and optionally save figure to a file" -function plot(EC::affine_EC, save_fig_filename=nothing) +function plot(EC::affine_EC, save_fig_filename=nothing; basis_points=nothing, basis_contour=nothing, xlims=nothing, ylims=nothing) scatter(real.(EC.training_E), imag.(EC.training_E), label="training") scatter!(real.(EC.exact_E), imag.(EC.exact_E), label="exact") scatter!(real.(EC.extrapolated_E), imag.(EC.extrapolated_E), label="extrapolated") + + isnothing(basis_points) || scatter!(real.(basis_points), imag.(basis_points), m=:x, label="basis") + isnothing(basis_contour) || plot!(real.(basis_contour), imag.(basis_contour), label="contour") + + isnothing(xlims) || xlims!(xlims...) + isnothing(ylims) || ylims!(ylims...) + isnothing(save_fig_filename) || savefig(save_fig_filename) end diff --git a/calculations/B2R_Berggren_poles.jl b/calculations/B2R_Berggren_poles.jl index 76c1a43..ba6967a 100644 --- a/calculations/B2R_Berggren_poles.jl +++ b/calculations/B2R_Berggren_poles.jl @@ -1,12 +1,12 @@ -using Plots +include("../EC.jl") include("../helper.jl") include("../p_space.jl") # contour p, w = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128]) -# ResonanceEC: Eq. (20) -V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) +μ = 0.5 +V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20) # generating a Berggren basis with a pole using the same system basis_c = 0.6 @@ -16,48 +16,18 @@ N_berg = sqrt.(diag(transpose(berg_basis .* w) * berg_basis)) berg_basis = berg_basis ./ transpose(N_berg) berg_basis_w = berg_basis .* w -training_points = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5) +H0 = transpose(berg_basis_w) * get_T_matrix(p, μ) * berg_basis +V = transpose(berg_basis_w) * get_V_matrix(V_system(1), p, w) * berg_basis -training_E = Vector{ComplexF64}(undef, length(training_points)) -EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points)) +training_c = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5) +extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8) -# training -for (j, c) in enumerate(training_points) - H = get_H_matrix(V_system(c), p, w) - H_berg = transpose(berg_basis_w) * H * berg_basis - evals, evecs = eigen(H_berg) - i = argmin(real.(evals)) -# i = identify_pole_i(basis_p, evals) - training_E[j] = evals[i] - EC_basis[:, j] = evecs[:, i] -end +training_ref = -0.26 +extrapolating_ref = [quick_pole_E(V_system(c)) for c in extrapolating_c] -EC_basis = hcat(EC_basis, conj.(EC_basis)) # CA-EC -EC_basis = gram_schmidt!(EC_basis) +EC = affine_EC(H0, V) +train!(EC, training_c; ref_eval=training_ref, CAEC=true) +extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref) -extrapolate_points = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8) - -exact_E = Vector{ComplexF64}(undef, length(extrapolate_points)) -extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points)) - -# extrapolating -for (j, c) in enumerate(extrapolate_points) - exact_E[j] = quick_pole_E(V_system(c)) - - H = get_H_matrix(V_system(c), p, w) - H_berg = transpose(berg_basis_w) * H * berg_basis - H_EC = transpose(EC_basis) * H_berg * EC_basis - evals = eigvals(H_EC) - i = argmin(abs.(evals .- exact_E[j])) - extrapolate_E[j] = evals[i] -end - -exportCSV("temp/2b_GSM_B2R.csv", (training_E, exact_E, extrapolate_E), ("training", "exact", "extrapolated")) - -scatter(real.(training_E), imag.(training_E), label="training") -scatter!(real.(exact_E), imag.(exact_E), label="exact") -scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated") -scatter!(real.(basis_E), imag.(basis_E), m=:x, label="Berggren basis") -xlims!(0,0.3) -ylims!(-0.120,0.020) -savefig("temp/2b_GSM_B2R.pdf") +exportCSV(EC, "temp/2b_GSM_B2R.csv") +plot(EC, "temp/2b_GSM_B2R.pdf"; basis_points=basis_E, xlims=(0, 0.3), ylims=(-0.120, 0.020)) diff --git a/calculations/B2R_comparison.jl b/calculations/B2R_comparison.jl index 2459e07..f063fab 100644 --- a/calculations/B2R_comparison.jl +++ b/calculations/B2R_comparison.jl @@ -1,52 +1,30 @@ -using Plots +include("../EC.jl") include("../helper.jl") include("../p_space.jl") berggren_mesh = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128]) csm_mesh = get_mesh([0, 8 - 3im], [512]) -for mesh in (berggren_mesh, csm_mesh) +for (mesh, name) in zip((berggren_mesh, csm_mesh), ("beggren", "csm")) p, w = mesh mesh_E = p.*p ./ (2*0.5) - # ResonanceEC: Eq. (20) - V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) + μ = 0.5 + V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20) + + H0 = get_T_matrix(p, μ) + V = get_V_matrix(V_system(1), p, w) + + training_c = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5) + extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8) - training_points = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5) - training_E = Vector{ComplexF64}(undef, length(training_points)) - EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points)) + training_ref = [quick_pole_E(V_system(c)) for c in training_c] + exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c] - for (j, c) in enumerate(training_points) - evals, evecs = eigen(get_H_matrix(V_system(c), p, w)) - i = identify_pole_i(p, evals) - training_E[j] = evals[i] - EC_basis[:, j] = evecs[:, i] - end - - EC_basis = hcat(EC_basis, conj.(EC_basis)) # CA-EC - EC_basis = gram_schmidt!(EC_basis, w) - EC_basis_w = EC_basis .* w - - extrapolate_points = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8) - - exact_E = Vector{ComplexF64}(undef, length(extrapolate_points)) - extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points)) - - for (j, c) in enumerate(extrapolate_points) - exact_E[j] = quick_pole_E(V_system(c)) - - H = get_H_matrix(V_system(c), p, w) - H_EC = transpose(EC_basis_w) * H * EC_basis - evals = eigvals(H_EC) - i = argmin(abs.(evals .- exact_E[j])) - extrapolate_E[j] = evals[i] - end - - scatter(real.(training_E), imag.(training_E), label="training") - scatter!(real.(exact_E), imag.(exact_E), label="exact") - scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated") - plot!(real.(mesh_E), imag.(mesh_E), label="contour") - xlims!(-0.3,0.3) - ylims!(-0.120,0.020) - savefig("temp/" * string(rand(UInt16)) * ".pdf") + EC = affine_EC(H0, V, w) + train!(EC, training_c; ref_eval=training_ref, CAEC=true) + extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E) + + #exportCSV(EC, "temp/2b_comparison_$name.csv") + plot(EC, "temp/2b_comparison_$name.pdf"; basis_contour=mesh_E, xlims=(-0.3,0.3), ylims=(-0.120,0.020)) end \ No newline at end of file diff --git a/calculations/R2R_Berggren.jl b/calculations/R2R_Berggren.jl index f5656f4..dc21130 100644 --- a/calculations/R2R_Berggren.jl +++ b/calculations/R2R_Berggren.jl @@ -1,4 +1,4 @@ -using Plots +include("../EC.jl") include("../helper.jl") include("../p_space.jl") @@ -8,42 +8,21 @@ subdivisions = [128, 128, 128] p, w = get_mesh(vertices, subdivisions) mesh_E = p.*p ./ (2*0.5) -# ResonanceEC: Eq. (20) -V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) +μ = 0.5 +V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20) -training_points = range(0.75, 0.45, 5) -training_E = Vector{ComplexF64}(undef, length(training_points)) -EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points)) +H0 = get_T_matrix(p, μ) +V = get_V_matrix(V_system(1), p, w) -for (j, c) in enumerate(training_points) - evals, evecs = eigen(get_H_matrix(V_system(c), p, w)) - i = identify_pole_i(p, evals) - training_E[j] = evals[i] - EC_basis[:, j] = evecs[:, i] -end +training_c = range(0.75, 0.45, 5) +extrapolating_c = range(0.40, 0.25, 5) -extrapolate_points = range(0.40, 0.25, 5) -ref_E = 0.2 - 0.1im +training_ref = [quick_pole_E(V_system(c)) for c in training_c] +exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c] -exact_E = Vector{ComplexF64}(undef, length(extrapolate_points)) -extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points)) +EC = affine_EC(H0, V, w) +train!(EC, training_c; ref_eval=training_ref, CAEC=false) +extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E) -EC_basis = gram_schmidt!(EC_basis, w) -EC_basis_w = EC_basis .* w - -for (j, c) in enumerate(extrapolate_points) - exact_E[j] = quick_pole_E(V_system(c)) - - H = get_H_matrix(V_system(c), p, w) - H_EC = transpose(EC_basis_w) * H * EC_basis - evals = eigvals(H_EC) - i = argmin(abs.(evals .- ref_E)) - global ref_E = evals[i] - extrapolate_E[j] = evals[i] -end - -scatter(real.(training_E), imag.(training_E), label="training") -scatter!(real.(exact_E), imag.(exact_E), label="exact") -scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated") -plot!(real.(mesh_E), imag.(mesh_E), label="contour") -xlims!(0,1) +#exportCSV(EC, "temp/2b_R2R.csv") +plot(EC, "temp/2b_R2R.pdf"; basis_contour=mesh_E, xlims=(0, 1)) diff --git a/calculations/R2R_Berggren_poles.jl b/calculations/R2R_Berggren_poles.jl index 1d68787..9e3afdd 100644 --- a/calculations/R2R_Berggren_poles.jl +++ b/calculations/R2R_Berggren_poles.jl @@ -1,66 +1,33 @@ -using Plots +include("../EC.jl") include("../helper.jl") include("../p_space.jl") # contour -p, w = get_mesh([0, 0.4 - 0.08im, 0.8, 6], [128, 128, 128]) -contour_E = p.^2 +p, w = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128]) -# ResonanceEC: Eq. (20) -V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) +μ = 0.5 +V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20) # generating a Berggren basis with a pole using the same system basis_c = 0.6 basis_E, berg_basis = eigen(get_H_matrix(V_system(basis_c), p, w); permute=false, scale=false) -pole_E = quick_pole_E(V_system(basis_c)) # basis only has 1 pole basis_p = sqrt.(basis_E) N_berg = sqrt.(diag(transpose(berg_basis .* w) * berg_basis)) berg_basis = berg_basis ./ transpose(N_berg) berg_basis_w = berg_basis .* w -training_points = range(0.79, 0.66, 4) # original: range(1.35, 0.9, 5) +H0 = transpose(berg_basis_w) * get_T_matrix(p, μ) * berg_basis +V = transpose(berg_basis_w) * get_V_matrix(V_system(1), p, w) * berg_basis -training_E = Vector{ComplexF64}(undef, length(training_points)) -EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points)) +training_c = range(0.79, 0.66, 4) # original: range(1.35, 0.9, 5) +extrapolating_c = range(0.62, 0.40, 6) # original: range(0.75, 0.40, 8) -# training -for (j, c) in enumerate(training_points) - H = get_H_matrix(V_system(c), p, w) - H_berg = transpose(berg_basis_w) * H * berg_basis - evals, evecs = eigen(H_berg) - i = identify_pole_i(basis_p, evals) - training_E[j] = evals[i] - EC_basis[:, j] = evecs[:, i] -end +training_ref = [quick_pole_E(V_system(c)) for c in training_c] +extrapolating_ref = [quick_pole_E(V_system(c)) for c in extrapolating_c] -EC_basis = gram_schmidt!(EC_basis) +EC = affine_EC(H0, V) +train!(EC, training_c; ref_eval=training_ref, CAEC=false) +extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref) -extrapolate_points = range(0.62, 0.40, 6) # original: range(0.75, 0.40, 8) - -exact_E = Vector{ComplexF64}(undef, length(extrapolate_points)) -extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points)) - -# extrapolating -for (j, c) in enumerate(extrapolate_points) - exact_E[j] = quick_pole_E(V_system(c)) - - H = get_H_matrix(V_system(c), p, w) - H_berg = transpose(berg_basis_w) * H * berg_basis - H_EC = transpose(EC_basis) * H_berg * EC_basis - evals = eigvals(H_EC) - i = argmin(abs.(evals .- exact_E[j])) - extrapolate_E[j] = evals[i] -end - -exportCSV("temp/2b_GSM_R2R.csv", (training_E, exact_E, extrapolate_E, [pole_E]), ("training", "exact", "extrapolated", "basis")) - -contour_E_export = contour_E[real.(contour_E) .< 1] # to trim unnecessary data outside axis limits -exportCSV("temp/2b_GSM_R2R_contour.csv", contour_E_export) - -scatter(real.(training_E), imag.(training_E), label="training") -scatter!(real.(exact_E), imag.(exact_E), label="exact") -scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated") -scatter!(real.(basis_E), imag.(basis_E), m=:x, label="Berggren basis") -xlims!(0,0.3) -ylims!(-0.120,0.020) -savefig("temp/2b_GSM_R2R.pdf") +exportCSV(EC, "temp/2b_GSM_R2R.csv") +plot(EC, "temp/2b_GSM_R2R.pdf"; basis_points=basis_E, xlims=(0, 0.3), ylims=(-0.120, 0.020)) diff --git a/calculations/XZ_technique.jl b/calculations/XZ_technique.jl index 98389e3..e74c3c0 100644 --- a/calculations/XZ_technique.jl +++ b/calculations/XZ_technique.jl @@ -1,5 +1,4 @@ -using LinearAlgebra, Plots -include("../helper.jl") +include("../EC.jl") include("../ho_basis.jl") include("../p_space.jl") @@ -22,35 +21,10 @@ n_EC = 8 train_cs = (0.7 .+ 0.05 * randn(n_EC)) - 1im * (0.2 .+ 0.05 * randn(n_EC)) target_cs = range(0.77, 0.22, 6) -train_E = zeros(ComplexF64, n_EC) -EC_basis = zeros(ComplexF64, (n_max + 1, length(train_cs))) -exact_E = zeros(ComplexF64, length(target_cs)) -extrapolate_E = similar(exact_E) - near_E = 0.2 + 0.2im -for (j, c) in enumerate(train_cs) - H = T + c .* V - evals, evecs = eigen(H) - i = argmin(abs.(evals .- near_E)) - train_E[j] = evals[i] - EC_basis[:, j] = evecs[:, i] -end +EC = affine_EC(T, V) +train!(EC, train_cs; ref_eval=near_E, CAEC=false) +extrapolate!(EC, target_cs) -EC_basis = gram_schmidt!(EC_basis) - -for (j, c) in enumerate(target_cs) - exact_E[j] = quick_pole_E((p, q) -> c*(V1*g0(R1, p, q) + V2*g0(R2, p, q)), μ; cs_angle=0.5) - - H = T + c .* V - H_EC = transpose(EC_basis) * H * EC_basis - evals = eigvals(H_EC) - i = argmin(abs.(evals .- exact_E[j])) - extrapolate_E[j] = evals[i] -end - -scatter(real.(train_E), imag.(train_E), label="training") -scatter!(real.(exact_E), imag.(exact_E), label="exact") -scatter!(real.(extrapolate_E), imag.(extrapolate_E), label="extrapolated") -xlims!(-0.2,0.3) -ylims!(-0.3,0.3) +plot(EC, "temp/XZ.pdf"; xlims=(-0.2,0.3), ylims=(-0.3,0.3))