From 29bbceac030c77a857ed04a1ff371f823662e2b8 Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Fri, 10 Jan 2025 15:35:16 -0500 Subject: [PATCH] EC.jl implemented for all 3-body systems --- EC.jl | 15 +++---- calculations/3body_Berggren_R2R_EC.jl | 63 +++++---------------------- calculations/3body_HO_B2R_EC.jl | 4 +- calculations/3body_HO_R2R_EC.jl | 58 ++++-------------------- calculations/3body_dis_HO_EC.jl | 61 +++----------------------- 5 files changed, 35 insertions(+), 166 deletions(-) diff --git a/EC.jl b/EC.jl index cf9e8c7..372cc2f 100644 --- a/EC.jl +++ b/EC.jl @@ -23,18 +23,18 @@ end If a list is provided for ref_eval, they are used as reference values for picking the closest eigenvalues at each sampling 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. Set orthonormalize_threshold=0 to skip Gram-Schmidt orthonormalization and use GEVP. Otherwise this value is used as the threshold for dropping redundant vectors." -function train!(EC::affine_EC, c_vals; ref_eval=nothing, orthonormalize_threshold=1e-5, CAEC=false, verbose=true, tol=1e-5) +function train!(EC::affine_EC, c_vals; ref_eval=-10.0, orthonormalize_threshold=1e-5, CAEC=false, verbose=true, tol=1e-5) training_vecs = Vector{ComplexF64}[] - current_E = -10.0 # randomly chosen reference value for c in c_vals verbose && println("Training for c = $c") + global current_E if ref_eval isa Number current_E = ref_eval ref_eval = nothing elseif !isnothing(ref_eval) - current_E = pop!(ref_eval) + current_E = popfirst!(ref_eval) end H = EC.H0 + c .* EC.H1 @@ -64,22 +64,21 @@ function train!(EC::affine_EC, c_vals; ref_eval=nothing, orthonormalize_threshol EC.trained = true end -"Extrapolate using a train EC model for a given range of c values +"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=nothing, verbose=true, tol=1e-5) +function extrapolate!(EC::affine_EC, c_vals; ref_eval=EC.training_E[end], verbose=true, tol=1e-5) @assert EC.trained "EC model must be trained using train() before extrapolation" - current_E = -10.0 # randomly chosen reference value - for c in c_vals verbose && println("Extact solution for c = $c") + global current_E if ref_eval isa Number current_E = ref_eval ref_eval = nothing elseif !isnothing(ref_eval) - current_E = pop!(ref_eval) + current_E = popfirst!(ref_eval) end H = EC.H0 + c .* EC.H1 diff --git a/calculations/3body_Berggren_R2R_EC.jl b/calculations/3body_Berggren_R2R_EC.jl index 431f20c..4c59beb 100644 --- a/calculations/3body_Berggren_R2R_EC.jl +++ b/calculations/3body_Berggren_R2R_EC.jl @@ -1,20 +1,20 @@ -using Plots +include("../EC.jl") training_c = [1.1, 0.9, 0.7, 0.5] extrapolating_c = 0.0 : 0.2 : 1.2 -training_ref = reverse([1.4750633616275919 - 0.0003021770706749637im +training_ref = [1.4750633616275919 - 0.0003021770706749637im 1.9567078295375822 - 0.0007646829108872369im 2.4351117758403076 - 0.001281037843108658im - 2.9096543462392357 - 0.002962488527470604im]) + 2.9096543462392357 - 0.002962488527470604im] -exact_ref = reverse([4.076662025307587-0.012709842443350328im, +extrapolating_ref = [4.076662025307587-0.012709842443350328im, 3.613318119833891-0.007335804709990623im, 3.1453431847006783-0.004030580410326795im, 2.672967129943755-0.00211498327461944im, 2.196542557810288-0.0010719835443437104im, 1.7164583929199813-0.0005455212208182736im, - 1.233088227541505-0.0003070320106485624im]) + 1.233088227541505-0.0003070320106485624im] include("../p_space_3body_resonance.jl") H0 = H @@ -32,54 +32,11 @@ end @time "Vp2" Vp2 = W_left * Vp2_HO * transpose(W_right) @time "Vp" Vp = Vpb + Vp2 -# free memory -basis = Hb_blocks = Hb = basis_ho = V2_HO = W_right = W_left = V2 = nothing -GC.gc() - -exact = ComplexF64[] -training = ComplexF64[] -extrapolated = ComplexF64[] -training_vecs = Vector{ComplexF64}[] - -for c in training_c - println("Training for c = $c") - global current_E = pop!(training_ref) - - local H = H0 + c .* Vp - local evals, evecs = eigs(H, sigma=current_E, maxiter=5000, tol=1e-5, ritzvec=true, check=1) - - global current_E = nearest(evals, current_E) - push!(training, current_E) - push!(training_vecs, evecs[:, nearestIndex(evals, current_E)]) -end - weights = repeat(kron(ws, ws), jmax + 1) -weights_mat = spdiagm(weights) -println("Original EC dimensionality = $(length(training_vecs))") -@time "Gram-Schmidt" training_vecs = gram_schmidt!(training_vecs, weights; verbose=true) # orthonormalization +EC = affine_EC(H0, Vp, weights) +train!(EC, training_c; ref_eval=training_ref, CAEC=false) +extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref) -EC_basis = hcat(training_vecs...) -H0_EC = transpose(EC_basis) * weights_mat * H0 * EC_basis -Vp_EC = transpose(EC_basis) * weights_mat * Vp * EC_basis - -for c in extrapolating_c - println("Extrapolating for c = $c") - global current_E = pop!(exact_ref) - - local H = H0 + c .* Vp - local evals, _ = eigs(H, sigma=current_E, maxiter=5000, tol=1e-5, ritzvec=false, check=1) - - global current_E = nearest(evals, current_E) - push!(exact, current_E) - - # extrapolation - H_EC = H0_EC + c .* Vp_EC - evals = eigvals(H_EC) - push!(extrapolated, nearest(evals, current_E)) -end - -scatter(real.(training),imag.(training), label="training") -scatter!(real.(exact),imag.(exact), label="exact") -scatter!(real.(extrapolated),imag.(extrapolated), label="extrapolated") -savefig("temp/Berggren_R2R.pdf") \ No newline at end of file +exportCSV(EC, "temp/Berggren_R2R.csv") +plot(EC, "temp/Berggren_R2R.pdf") \ No newline at end of file diff --git a/calculations/3body_HO_B2R_EC.jl b/calculations/3body_HO_B2R_EC.jl index 779c89c..cd4acee 100644 --- a/calculations/3body_HO_B2R_EC.jl +++ b/calculations/3body_HO_B2R_EC.jl @@ -1,7 +1,7 @@ include("../EC.jl") training_ref = -0.72763 -exact_ref = 4.0766890719636635 - 0.01275892774109674im +extrapolating_ref = 4.0766890719636635 - 0.01275892774109674im training_c = [2.0, 1.9, 1.8] extrapolating_c = 0.0 : 0.2 : 1.2 @@ -15,7 +15,7 @@ Vp_of_r(r) = -exp(-(r/3)^2) EC = affine_EC(H0, Vp) train!(EC, training_c; ref_eval=training_ref, CAEC=true) -extrapolate!(EC, extrapolating_c; ref_eval=exact_ref) +extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref) exportCSV(EC, "temp/HO_B2R.csv") plot(EC, "temp/HO_B2R.pdf") diff --git a/calculations/3body_HO_R2R_EC.jl b/calculations/3body_HO_R2R_EC.jl index 3a9aba9..6c85598 100644 --- a/calculations/3body_HO_R2R_EC.jl +++ b/calculations/3body_HO_R2R_EC.jl @@ -1,56 +1,16 @@ -using DelimitedFiles, Plots -include("../ho_basis_3body_resonance.jl") +include("../EC.jl") -current_E = 5.9673 - 0.0006im +ref_E = 5.9673 - 0.0006im training_c = 2.0 : -0.2 : 1.2 extrapolating_c = 1.05 .- [0.0 : 0.1 : 0.4; 0.45 : 0.05 : 0.60] -@time "H0" H0 = T1 + T2 +include("../ho_basis_3body_resonance.jl") +@time "H0" H0 = T1 + T2 + T_cross -# free memory -basis = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing -GC.gc() +EC = affine_EC(H0, V) +train!(EC, training_c; ref_eval=ref_E, CAEC=false) +extrapolate!(EC, extrapolating_c) -exact = ComplexF64[] -training = ComplexF64[] -extrapolated = ComplexF64[] -training_vecs = Vector{ComplexF64}[] - -for c in training_c - println("Training for c = $c") - local H = H0 + c .* V - local evals, evecs = eigs(H, nev=3, ncv=24, which=:LI, maxiter=5000, tol=1e-5, ritzvec=true, check=1) - - global current_E = nearest(evals, current_E) - push!(training, current_E) - push!(training_vecs, evecs[:, nearestIndex(evals, current_E)]) -end - -println("Original EC dimensionality = $(length(training_vecs))") -@time "Gram-Schmidt" training_vecs = gram_schmidt!(training_vecs; verbose=true) # orthonormalization - -EC_basis = hcat(training_vecs...) -H0_EC = transpose(EC_basis) * H0 * EC_basis -V_EC = transpose(EC_basis) * V * EC_basis - -for c in extrapolating_c - println("Extrapolating for c = $c") - local H = H0 + c .* V - local evals, evecs = eigs(H, nev=3, ncv=24, which=:LI, maxiter=5000, tol=1e-5, ritzvec=true, check=1) - - global current_E = nearest(evals, current_E) - push!(exact, current_E) - - # extrapolation - H_EC = H0_EC + c .* V_EC - evals = eigvals(H_EC) - push!(extrapolated, nearest(evals, current_E)) -end - -exportCSV("temp/NCSM.csv", (training, exact, extrapolated), ("training", "exact", "extrapolated")) - -scatter(real.(training),imag.(training), label="training") -scatter!(real.(exact),imag.(exact), label="exact") -scatter!(real.(extrapolated),imag.(extrapolated), label="extrapolated") -savefig("temp/NCSM.pdf") \ No newline at end of file +exportCSV(EC, "temp/HO_R2R.csv") +plot(EC, "temp/HO_R2R.pdf") diff --git a/calculations/3body_dis_HO_EC.jl b/calculations/3body_dis_HO_EC.jl index 523d1b6..925fa58 100644 --- a/calculations/3body_dis_HO_EC.jl +++ b/calculations/3body_dis_HO_EC.jl @@ -1,6 +1,5 @@ -using Arpack, SparseArrays, LRUCache -using DelimitedFiles, Plots include("../ho_basis.jl") +include("../EC.jl") Λ = 0 m = 1.0 @@ -31,62 +30,16 @@ println("Basis size = ", basis.dim) @time "Ha" Ha = T1 + T2 + Va @time "Eigenvalues" target_evals, _ = eigs(Ha, nev=5, ncv=50, which=:SR, maxiter=5000, tol=1e-5, ritzvec=false, check=1) - display(target_evals) -# free memory -basis = T1 = T2 = V1_cache = V_relative_cache = V1 = V_relative = U = V2 = nothing -GC.gc() - training_c = [-0.5, -0.65, -0.8, -1, -1.2] extrapolating_c = [0.8, 0.6, 0.4, 0.2, 0.1, 0.0, -0.1, -0.2, -0.3] -current_E = -0.5173809356244544 +ref_E = -0.5173809356244544 -exact = ComplexF64[] -training = ComplexF64[] -extrapolated = ComplexF64[] -training_vecs = Vector{ComplexF64}[] +EC = affine_EC(Ha, Vb) +train!(EC, training_c; ref_eval=ref_E, CAEC=true) # try CAEC=false !!! +extrapolate!(EC, extrapolating_c) -for c in training_c - print("Training for c = $c: ") - H = Ha + c .* Vb - evals, evecs = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=true, check=1) - - global current_E = nearest(evals, current_E) - println(current_E) - push!(training, current_E) - push!(training_vecs, evecs[:, nearestIndex(evals, current_E)]) -end - -training_vecs = vcat(training_vecs, conj(training_vecs)) # CA-EC -println("Original EC dimensionality = $(length(training_vecs))") -@time "Gram-Schmidt" training_vecs = gram_schmidt!(training_vecs; verbose=true) # orthonormalization - -EC_basis = hcat(training_vecs...) -Ha_EC = transpose(EC_basis) * Ha * EC_basis -Vb_EC = transpose(EC_basis) * Vb * EC_basis - -current_E = -0.3005521915662689 - 0.13612069020686351im - -for c in extrapolating_c - print("Extrapolating for c = $c: ") - H = Ha + c .* Vb - evals, evecs = eigs(H, nev=3, ncv=24, which=:SR, maxiter=5000, tol=1e-5, ritzvec=true, check=1) - - global current_E = nearest(evals, current_E) - println(current_E) - push!(exact, current_E) - - # extrapolation - H_EC = Ha_EC + c .* Vb_EC - evals = eigvals(H_EC) - push!(extrapolated, nearest(evals, current_E)) -end - -# exportCSV("temp/dis_HO_B2R.csv", (training, exact, extrapolated), ("training", "exact", "extrapolated")) - -scatter(real.(training),imag.(training), label="training") -scatter!(real.(exact),imag.(exact), label="exact") -scatter!(real.(extrapolated),imag.(extrapolated), label="extrapolated") -savefig("temp/dis_HO_B2R.pdf") \ No newline at end of file +exportCSV(EC, "temp/dis_HO_B2R.csv") +plot(EC, "temp/dis_HO_B2R.pdf") \ No newline at end of file