using LinearAlgebra, DelimitedFiles, SparseArrays

@enum coordinate_system jacobi src

"Sum over array while minimizing catastrophic cancellation as much as possible"
function better_sum(arr::Array{T}) where T<:Real
    pos_arr = arr[arr .> 0]
    neg_arr = arr[arr .< 0]

    sort!(pos_arr)
    sort!(neg_arr, rev=true)

    return sum(pos_arr) + sum(neg_arr)
end

better_sum(arr::Array{ComplexF64}) = better_sum(real.(arr)) + 1im * better_sum(imag.(arr))

"The triangle inequality for angular momenta"
triangle_ineq(l1, l2, L) = abs(l1 - l2) ≤ L && L ≤ (l1 + l2)

"Index of the nearest value in a list to a given reference point"
nearestIndex(list::Array, ref) = argmin(norm.(list .- ref))

"Nearest value in a list to a given reference point"
nearest(list::Array, ref) = list[nearestIndex(list, ref)]

"Simple implementation of the Kronecker sum"
function kron_sum(A::AbstractMatrix, B::AbstractMatrix)
    @assert size(A, 1) == size(A, 2) && size(B, 1) == size(B, 2) "Matrices should be square"
    return kron(A, I(size(B, 1))) + kron(I(size(A, 1)), B)
end

"Flattened vector version of Iterators.product(...) with index hierachy reversed -- leftmost index has the highest hierachy"
iter_prod(args...) = reverse.(collect(Iterators.product(reverse(args)...))[:])

"Export CSV data for a scatter plot taking in data as a list of complex vectors (x=Re, y=Im), and a list of corresponding labels, typically used for EC results"
function exportCSV(file::String, data, labels=nothing)
    columns = ["x" "y" "label"]

    if isnothing(labels)
        columns[end] = ""
        labels = fill("", length(data))
    end
    
    open(file, "w") do f
        writedlm(f, columns)
        for (d, label) in zip(data, labels)
            l = fill(label, length(d))
            writedlm(f, hcat(reim(d)..., l))
        end
    end
end

"In-place c-orthogonalization via (modified) Gram-Schmidt. Only significant vectors are returned (c-normalized).
The number of significant vectors to return are determined by the original singular values (compared to the threshold)."
function gram_schmidt!(vecs::Vector{Vector{T}}, ws=ones(length(vecs[1])), threshold=1e-5; verbose=false) where T<: Number
    c_product(i, j) = sum(vecs[i] .* ws .* vecs[j])
    norm(i) = c_product(i, i)
    proj(i, j) = (c_product(i, j) / norm(i)) .* vecs[i] # component of vec[i] in vec[j]

    # initial normalization
    for i in eachindex(vecs)
        vecs[i] ./= sqrt(norm(i))
    end

    verbose && println("Absolute singular values = $(round.(c_singular_values(vecs, ws); sigdigits=1))")
    target_dim = c_rank(vecs, ws, threshold)
    verbose && println("Target dimensionality = $target_dim")

    selected_vecs_i = Integer[]

    while length(selected_vecs_i) < target_dim
        i = argmax(i -> abs(norm(i)), setdiff(eachindex(vecs), selected_vecs_i)) # find the largest vector from the remaining
        push!(selected_vecs_i, i)
        for j in setdiff(eachindex(vecs), selected_vecs_i)
            vecs[j] .-= proj(i, j)
        end
    end
    
    verbose && println("Absolute norms of selected vectors = $(round.(selected_vecs_i .|> norm .|> abs; sigdigits=1))")
    verbose && println("Absolute norms of dropped vectors = $(round.(setdiff(eachindex(vecs), selected_vecs_i) .|> norm .|> abs; sigdigits=1))")

    # final normalization
    for i in selected_vecs_i
        vecs[i] ./= sqrt(norm(i))
    end

    return vecs[selected_vecs_i]
end

"Same as above but the basis is provided as a matrix instead of an array of vectors"
gram_schmidt!(vecs::Matrix{T}, ws=ones(size(vecs, 1)), threshold=1e-5; verbose=false) where T<: Number = hcat(gram_schmidt!([vecs[:, i] for i in axes(vecs, 2)], ws, threshold; verbose=verbose)...)

"Rank of a basis set (pre-normalized) under c-product"
c_rank(vecs, ws, threshold=1e-8) = count(c_singular_values(vecs, ws) .> threshold)

"Singular values (magnitudes) of a basis set (pre-normalized) under c-product"
function c_singular_values(vecs, ws)
    basis = hcat(vecs...)
    N = transpose(basis) * spdiagm(ws) * basis
    singular_values = eigvals(N) .|> abs .|> sqrt
    return singular_values
end