using LinearAlgebra, Plots
include("ho_basis.jl")
include("p_space.jl")

μ = 0.5
l = 0
V1 = -8.6
R1 = 1.75
V2 = 3.4
R2 = 3.2
n_max = 20

ns = collect(0:n_max)
ls = fill(l, n_max + 1)

T = 1/μ * sp_T_matrix(ns, ls)
V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns)) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns))

n_EC = 16
train_cs = (2.2 - 0.7im) .+ (0.1 .* randn(n_EC)) .+ (0.1im * randn(n_EC))
target_cs = range(2.0, 0.8, 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_bound = 1.5 + 1.2im

for (j, c) in enumerate(train_cs)
    H = T + c .* V
    evals, evecs = eigen(H)
    i = argmin(abs.(evals .- near_E_bound))
    train_E[j] = evals[i]
    EC_basis[:, j] = evecs[:, i]
end

N_EC = transpose(EC_basis) * EC_basis

near_E_res = 1.5 - 0.4im

for (j, c) in enumerate(target_cs)
    mesh_p, mesh_w = get_mesh([0, 8 * exp(-0.5im)], 256)
    p_space_H = get_H_matrix((p, q) -> c*(V1*g0(R1, p, q) + V2*g0(R2, p, q)), mesh_p, mesh_w, μ)
    evals = eigvals(p_space_H)
    i = argmin(abs.(evals .- near_E_res))
    exact_E[j] = evals[i]

    H = T + c .* V
    H_EC = transpose(EC_basis) * H * EC_basis
    evals = eigvals(H_EC, N_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.5,2.0)
ylims!(-1.5,1.5)