EC.jl implemented for all

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}
+    H0::AbstractMatrix{ComplexF64}
-    H1::SparseMatrixCSC{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
 
+        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

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)
+μ = 0.5
-V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
+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))
+training_c = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5)
-EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))
+extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
 
-# training
+training_ref = -0.26
-for (j, c) in enumerate(training_points)
+extrapolating_ref = [quick_pole_E(V_system(c)) for c in extrapolating_c]
-    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
 
-EC_basis = hcat(EC_basis, conj.(EC_basis)) # CA-EC
+EC = affine_EC(H0, V)
-EC_basis = gram_schmidt!(EC_basis)
+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)
+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))
-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")

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)
+    μ = 0.5
-    V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
+    V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q)) # ResonanceEC: Eq. (20)
     
-    training_points = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5)
+    H0 = get_T_matrix(p, μ)
-    training_E = Vector{ComplexF64}(undef, length(training_points))
+    V = get_V_matrix(V_system(1), p, w)
-    EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))
     
-    for (j, c) in enumerate(training_points)
+    training_c = range(1.1, 0.9, 5) # original: range(1.35, 0.9, 5)
-        evals, evecs = eigen(get_H_matrix(V_system(c), p, w))
+    extrapolating_c = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
-        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
+    training_ref = [quick_pole_E(V_system(c)) for c in training_c]
-    EC_basis = gram_schmidt!(EC_basis, w)
+    exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c]
-    EC_basis_w = EC_basis .* w
 
-    extrapolate_points = range(0.78, 0.45, 7) # original: range(0.75, 0.40, 8)
+    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)
     
-    exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))
+    #exportCSV(EC, "temp/2b_comparison_$name.csv")
-    extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))
+    plot(EC, "temp/2b_comparison_$name.pdf"; basis_contour=mesh_E, xlims=(-0.3,0.3), ylims=(-0.120,0.020))
-
-    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")
 end

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)
+μ = 0.5
-V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
+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)
+H0 = get_T_matrix(p, μ)
-training_E = Vector{ComplexF64}(undef, length(training_points))
+V = get_V_matrix(V_system(1), p, w)
-EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))
 
-for (j, c) in enumerate(training_points)
+training_c = range(0.75, 0.45, 5)
-    evals, evecs = eigen(get_H_matrix(V_system(c), p, w))
+extrapolating_c = range(0.40, 0.25, 5)
-    i = identify_pole_i(p, evals)
-    training_E[j] = evals[i]
-    EC_basis[:, j] = evecs[:, i]
-end
 
-extrapolate_points = range(0.40, 0.25, 5)
+training_ref = [quick_pole_E(V_system(c)) for c in training_c]
-ref_E = 0.2 - 0.1im
+exact_E = [quick_pole_E(V_system(c)) for c in extrapolating_c]
 
-exact_E = Vector{ComplexF64}(undef, length(extrapolate_points))
+EC = affine_EC(H0, V, w)
-extrapolate_E = Vector{ComplexF64}(undef, length(extrapolate_points))
+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)
+#exportCSV(EC, "temp/2b_R2R.csv")
-EC_basis_w = EC_basis .* w
+plot(EC, "temp/2b_R2R.pdf"; basis_contour=mesh_E, xlims=(0, 1))
-
-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)

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])
+p, w = get_mesh([0, 0.4 - 0.15im, 0.8, 6], [128, 128, 128])
-contour_E = p.^2
 
-# ResonanceEC: Eq. (20)
+μ = 0.5
-V_system(c) = (p, q) -> c*(-5*g0(sqrt(3), p, q) + 2*g0(sqrt(10), p, q))
+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))
+training_c = range(0.79, 0.66, 4) # original: range(1.35, 0.9, 5)
-EC_basis = Matrix{ComplexF64}(undef, length(p), length(training_points))
+extrapolating_c = range(0.62, 0.40, 6) # original: range(0.75, 0.40, 8)
 
-# training
+training_ref = [quick_pole_E(V_system(c)) for c in training_c]
-for (j, c) in enumerate(training_points)
+extrapolating_ref = [quick_pole_E(V_system(c)) for c in extrapolating_c]
-    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
 
-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)
+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))
-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")

calculations/XZ_technique.jl

@@ -1,5 +1,4 @@
-using LinearAlgebra, Plots
+include("../EC.jl")
-include("../helper.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)
+EC = affine_EC(T, V)
-    H = T + c .* V
+train!(EC, train_cs; ref_eval=near_E, CAEC=false)
-    evals, evecs = eigen(H)
+extrapolate!(EC, target_cs)
-    i = argmin(abs.(evals .- near_E))
-    train_E[j] = evals[i]
-    EC_basis[:, j] = evecs[:, i]
-end
 
-EC_basis = gram_schmidt!(EC_basis)
+plot(EC, "temp/XZ.pdf"; xlims=(-0.2,0.3), ylims=(-0.3,0.3))
-
-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)