BergEC-jl/EC.jl

using Statistics, SparseArrays, LinearAlgebra, Arpack, Plots
include("common.jl")

"EC model for a Hamiltonian family H(c) = H0 + c * H1"
mutable struct affine_EC
    H0::AbstractMatrix{ComplexF64}
    H1::AbstractMatrix{ComplexF64}
    weights::Vector{ComplexF64}
    trained::Bool

    H0_EC
    H1_EC
    N_EC

    ensemble_size::Int
    H0_EC_ensemble
    H1_EC_ensemble
    N_EC_ensemble

    training_E::Vector{ComplexF64}
    exact_E::Vector{ComplexF64}
    extrapolated_E::Vector{ComplexF64}
    extrapolated_CI::Vector{ComplexF64}

    affine_EC(H0::AbstractMatrix{ComplexF64}, H1::AbstractMatrix{ComplexF64}, weights::Vector{ComplexF64}=ones(ComplexF64, size(H0, 1)); ensemble_size=0) =
        new(H0, H1, weights, false, nothing, nothing, nothing, ensemble_size, Matrix[], Matrix[], Matrix[], ComplexF64[], ComplexF64[], ComplexF64[], ComplexF64[])
end

"Train an 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 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.
If orthonormalize_threshold > 0, Gram-Schmidt orthonormalization is performed, using this value as the threshold for dropping redundant vectors."
function train!(EC::affine_EC, c_vals; ref_eval=-10.0, CAEC=false, gram_schmidt_threshold=0, verbose=true, tol=1e-5)
    training_vecs = Vector{ComplexF64}[]

    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 = popfirst!(ref_eval)
        end

        H = EC.H0 + c .* EC.H1

        evals, evecs = eigs(H, sigma=current_E, maxiter=5000, tol=tol, ritzvec=true, check=1)
        current_E = nearest(evals, current_E)

        push!(EC.training_E, current_E)
        push!(training_vecs, evecs[:, nearestIndex(evals, current_E)])
    end

    CAEC && append!(training_vecs, conj.(training_vecs))

    (EC.H0_EC, EC.H1_EC, EC.N_EC) = get_reduced_matrices(EC, training_vecs, gram_schmidt_threshold; verbose=true)

    for _ in 1:EC.ensemble_size
        (H0_EC, H1_EC, N_EC) = get_reduced_matrices(EC, training_vecs, gram_schmidt_threshold, resample(length(training_vecs)); verbose=false)
        push!(EC.H0_EC_ensemble, H0_EC)
        push!(EC.H1_EC_ensemble, H1_EC)
        push!(EC.N_EC_ensemble, N_EC)
    end

    EC.trained = true
end

function get_reduced_matrices(EC::affine_EC, training_vecs, gram_schmidt_threshold, subsample=1:length(training_vecs); verbose=false)
    vecs = deepcopy(training_vecs[subsample])

    if gram_schmidt_threshold > 0; vecs = gram_schmidt!(sampled_vecs, EC.weights, gram_schmidt_threshold; verbose=verbose); end

    EC_basis = hcat(orth_vecs...)
    weights_mat = spdiagm(EC.weights)
    H0_EC = transpose(EC_basis) * weights_mat * EC.H0 * EC_basis
    H1_EC = transpose(EC_basis) * weights_mat * EC.H1 * EC_basis
    N_EC = transpose(EC_basis) * weights_mat * EC_basis

    return (H0_EC, H1_EC, N_EC)
end

resample(n::Int) = rand(1:n, n) |> unique |> sort

"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.
If precalculated_exact_E is provided, ref_eval is ignored.
If pseudo_inv_rtol > 0, the GEVP is avoided using Moore-Penrose psuedoinverse, using this value as the relative tolerance for dropping redundant vectors."
function extrapolate!(EC::affine_EC, c_vals; ref_eval=EC.training_E[end], pseudo_inv_rtol=1e-6, 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
        global current_E

        if isnothing(precalculated_exact_E)
            if ref_eval isa Number
                current_E = ref_eval
                ref_eval = nothing
            elseif !isnothing(ref_eval)
                current_E = popfirst!(ref_eval)
            end

            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)

        verbose && println("Extrapolating for c = $c")

        evals = get_extrapolated_evals(EC.H0_EC, EC.H1_EC, EC.N_EC, c, pseudo_inv_rtol)
        push!(EC.extrapolated_E, nearest(evals, current_E))

        if EC.ensemble_size > 0
            E_ensemble = zeros(ComplexF64, EC.ensemble_size)
            for i in 1:EC.ensemble_size
                evals = get_extrapolated_evals(EC.H0_EC_ensemble[i], EC.H1_EC_ensemble[i], EC.N_EC_ensemble[i], c, pseudo_inv_rtol)
                E_ensemble[i] = nearest(evals, current_E)
            end
            re_CI = std(real.(E_ensemble))
            im_CI = std(imag.(E_ensemble))
            push!(EC.extrapolated_CI, complex(re_CI, im_CI))
        end
    end
end

"Solve the GEVP with or without Moore-Penrose psuedoinverse"
function get_extrapolated_evals(H0_EC, H1_EC, N_EC, c, pseudo_inv_rtol)
    H_EC = H0_EC + c .* H1_EC
    if pseudo_inv_rtol > 0
        inv_N_EC = pinv(N_EC; rtol=pseudo_inv_rtol)
        H_EC = inv_N_EC * H_EC
        return eigvals(H_EC)
    else
        return eigvals(H_EC, N_EC)
    end
end

"Export EC data as CSV file"
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; 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")
    if EC.ensemble_size > 0
        scatter!(real.(EC.extrapolated_E), imag.(EC.extrapolated_E), xerror=real.(EC.extrapolated_CI), yerror=imag.(EC.extrapolated_CI), label="extrapolated")
    else
        scatter!(real.(EC.extrapolated_E), imag.(EC.extrapolated_E), label="extrapolated")
    end

    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