text:
"""
            Functions for sph mapping to grid.

            sphMapping_2D without property conservation works like SPLASH by
            Daniel Price: http://users.monash.edu.au/~dprice/splash/.

            The property conservation is based on Smac by Dolag et. al. 2005:
            https://ui.adsabs.harvard.edu/abs/2005MNRAS.363...29D/abstract

            sphAdaptiveMapping is adapted from SPHMapper by Dr. Alexander Arth
            (private conversations). A similar implementation is used in Pygad.
            Publication: https://arxiv.org/abs/1803.03652

    Author: Ludwig Böss
    Contact: lboess@usm.lmu.de
    Created: 2018-12-12

"""

using Distributed
using Statistics
using ProgressMeter
using Base.Threads

@inline function find_position_periodic(pos::Vector{Float64}, k::Int64, bosize::Float64)

    x = !(k & 0x1) ? pos[1] : ( pos[1] > 0.0 ? pos[1] - boxsize : pos[1] + boxsize)
    y = !(k & 0x2) ? pos[2] : ( pos[2] > 0.0 ? pos[2] - boxsize : pos[2] + boxsize)
    z = !(k & 0x4) ? pos[3] : ( pos[3] > 0.0 ? pos[3] - boxsize : pos[3] + boxsize)

    return [x, y, z]
end

@inline function get_d_hsml_2D(dx::Float64, dy::Float64, hsml_inv::Float64)
    sqrt( dx*dx + dy*dy ) * hsml_inv
end

@inline function get_d_hsml_3D(dx::Float64, dy::Float64, dz::Float64,
                               hsml_inv::Float64)
    sqrt( dx*dx + dy*dy + dz*dz ) * hsml_inv
end

@inline function check_in_image(pos::Float64, hsml::Float64,
                                minCoords::Float64, maxCoords::Float64)

    if ( (minCoords - hsml) <= pos <= (maxCoords + hsml) )
        return true
    else
        return false
    end
end

@inline function find_min_pixel(pos::Float64, hsml::Float64,
                                minCoords::Float64,
                                pixsize_inv::Float64)

    pix = floor(Int64, (pos - hsml - minCoords) * pixsize_inv )

    return max(pix, 1)
end

@inline function find_max_pixel(pos::Float64, hsml::Float64,
                                minCoords::AbstractFloat, pixsize_inv::Float64,
                                max_pixel::Integer)

    pix = floor(Int64, (pos + hsml - minCoords) * pixsize_inv )

    return min(pix, max_pixel)
end

function sphMapping_2D(Pos, HSML, M, ρ, Bin_Quant;
                       param::mappingParameters, kernel,
                       conserve_quantities::Bool=false,
                       show_progress::Bool=false)

    # if this is not a float it has units, which need to be stripped
    if !(typeof(Bin_Quant[1,1]) <: AbstractFloat)

        if show_progress
            @info "Stripping units..."
        end

        Pos       = ustrip(Pos)
        HSML      = ustrip(HSML)
        M         = ustrip(M)
        ρ         = ustrip(ρ)
        Bin_Quant = ustrip(Bin_Quant)

    end

    N = length(M)  # number of particles

    image = zeros(length(param.x), length(param.y))

    minCoords = [param.x[1], param.y[1], param.z[1]]
    maxCoords = [param.x[end], param.y[end], param.z[end]]

    max_pixel = [length(param.x), length(param.y)]

    particles_in_image = 0

    if show_progress
        P = Progress(N)
        idx = 0
        #P_lock = SpinLock()  # uncomment to make thread-safe if needed in the future
    end

    @inbounds for p = 1:N

        # save stuff from array to single variables
        pos  = Float64.(Pos[p,:])
        hsml = Float64(HSML[p])

        in_image = false

        @inbounds for dim = 1:3

            in_image = check_in_image(pos[dim], hsml,
                                      minCoords[dim], maxCoords[dim])

            # exit the loop if the particle is not in the image frame
            if !in_image
                break
            end
        end

        # only calculate the properties if the particle is in the image
        if in_image

            # store this here for performance increase
            pixsize_inv = 1.0/param.pixelSideLength

            particles_in_image += 1

            # save rest of variables
            bin_q       = Float64(Bin_Quant[p])
            hsml_inv    = Float64(1.0/hsml)
            m           = Float64(M[p])
            rho_inv     = Float64(1.0/ρ[p])

            pixmin = Vector{Int64}(undef,2)
            pixmax = Vector{Int64}(undef,2)

            @inbounds for dim = 1:2

                pixmin[dim] = find_min_pixel(pos[dim], hsml, minCoords[dim],
                                             pixsize_inv)

                pixmax[dim] = find_max_pixel(pos[dim], hsml, minCoords[dim],
                                             pixsize_inv, max_pixel[dim])

            end

            if conserve_quantities

                # calculate pixel area
                xp1 = (pos[1] - hsml) * pixsize_inv
                xp2 = (pos[1] + hsml) * pixsize_inv
                yp1 = (pos[2] - hsml) * pixsize_inv
                yp2 = (pos[2] + hsml) * pixsize_inv

                pix_area = (xp2 - xp1) * (yp2 - yp1) * param.pixelArea

                d3 = (m * rho_inv) / ( pix_area * param.pixelSideLength)

                # number of pixels over which the particle is distributed
                N_distr = Int64((pixmax[1] - pixmin[1] + 1 ) * (pixmax[2] - pixmin[2] + 1))

                # allocate arrays for weights
                kernel_tab = zeros(N_distr)
                d1_tab     = zeros(N_distr)
                d2_tab     = zeros(N_distr)
                dx_tab     = zeros(N_distr)

                # weights table
                wit1    = 0.
                witd    = 0.
                wit2    = 0
                wit2tot = 0

                # pixel count for attributing weights to pixels
                N_count = 1

                # first loop to calculate weights
                @inbounds for i = pixmin[2]:pixmax[2]
                    dy = param.y[i] - pos[2]
                    djmin = max(yp1, i - 1.0)
                    djmax = min(yp2, i)

                    @inbounds for j = pixmin[1]:pixmax[1]

                        dx = param.x[j] - pos[1]
                        dimin = max(xp1, j - 1.0)
                        dimax = min(xp2, j)

                        d1 = dimax - dimin
                        d2 = djmax - djmin

                        # compute distance to pixel center in units of hsml
                        dx_tab[N_count] = get_d_hsml_2D(dx, dy, hsml_inv)
                        # update pixel value
                        wi = kernel_value_2D(kernel, dx_tab[N_count], hsml_inv)

                        wit1 += wi * d1 * d2
                        witd += d1 * d2

                        if wi > 0
                            wit2 += 1
                        end

                        # store data in arrays
                        d1_tab[N_count]     = d1
                        d2_tab[N_count]     = d2
                        kernel_tab[N_count] = wi

                        wit2tot += 1

                        N_count += 1

                    end # end x-loop
                end # end y-loop

                #fak = 1.0/( (pixmax[2] - pixmin[2]) * (pixmax[1] - pixmin[1]) )
                if wit1 > 0
                    fak = wit2/wit1
                else
                    wit2 = wit2tot
                    fak = wit2tot / witd
                end

                fak_hsml = pix_area/(wit2 * param.pixelArea)

                # reset the counter
                N_count = 1

            end # conserve_quantities

            # to reduce number of calculations. Only needed if quantity conservation
            # is switched off. Otherwise it is incapsulated in d3
            if !conserve_quantities
                bin_prefac = bin_q * m * rho_inv
            end

            # second loop to calculate value
            @inbounds for i = pixmin[2]:pixmax[2]

                @inbounds for j = pixmin[1]:pixmax[1]

                    if !conserve_quantities

                        # calculate simple distance to pixel center
                        dx = param.x[j] - pos[1]
                        dy = param.y[i] - pos[2]
                        distance_hsml = get_d_hsml_2D(dx, dy, hsml_inv)

                        # update pixel value
                        image[j, i] += bin_prefac * kernel_value_2D(kernel, distance_hsml, hsml_inv)
                    else

                        if wit1 <= 1
                            wi = fak * fak_hsml
                        else
                            wi = kernel_tab[N_count] * fak * fak_hsml
                        end
                        # update pixel value with weights
                        image[j, i] += bin_q * wi *
                                       d1_tab[N_count] * d2_tab[N_count] * d3
                        N_count += 1
                    end

                end # end x-loop
            end # end y-loop

        end # end check if in image

        # update for ProgressMeter
        if show_progress
            #lock(P_lock)  # uncomment to make thread-safe if needed in the future
            idx += 1
            ProgressMeter.update!(P, idx)
            #unlock(P_lock)
        end
    end

    if show_progress
        @info "Mapped $particles_in_image / $N particles."
    end

    # tilt image
   return copy(transpose(image))
end

function sphMapping_3D(Pos, HSML, M, ρ, Bin_Quant;
                       param::mappingParameters, kernel::SPHKernel,
                       conserve_quantities::Bool=false,
                       show_progress::Bool=false)

    # if this is not a float it has units, which need to be stripped
    if !(typeof(Bin_Quant[1,1]) <: AbstractFloat)

        if show_progress
            @info "Stripping units..."
        end

        Pos       = ustrip(Pos)
        HSML      = ustrip(HSML)
        M         = ustrip(M)
        ρ         = ustrip(ρ)
        Bin_Quant = ustrip(Bin_Quant)

    end

    N = length(M)  # number of particles

    image = zeros(length(param.x), length(param.y), length(param.z))

    minCoords = [param.x[1], param.y[1], param.z[1]]
    maxCoords = [param.x[end], param.y[end], param.z[end]]

    max_pixel = [length(param.x), length(param.y), length(param.z)]

    particles_in_image = 0

    if show_progress
        P = Progress(N)
        idx = 0
        #P_lock = SpinLock()  # uncomment to make thread-safe if needed in the future
    end

    @inbounds for p = 1:N

        # save stuff from array to single variables
        bin_q = Float64(Bin_Quant[p])
        pos   = Float64.(Pos[p,:])
        hsml  = Float64(HSML[p])

        in_image = false

        @inbounds for dim = 1:3

            in_image = check_in_image(pos[dim], hsml,
                                      minCoords[dim], maxCoords[dim])

            # exit the loop if the particle is not in the image frame
            if !in_image
                break
            end
        end

        # only calculate the properties if the particle is in the image
        if in_image

            # store this here for performance increase
            pixsize_inv = 1.0/param.pixelSideLength

            particles_in_image += 1

            # save rest of variables
            hsml_inv    = Float64(1.0/hsml)
            m           = Float64(M[p])
            rho_inv     = Float64(1.0/ρ[p])

            pixmin = Vector{Int64}(undef,3)
            pixmax = Vector{Int64}(undef,3)

            @inbounds for dim = 1:3

                pixmin[dim] = find_min_pixel(pos[dim], hsml, minCoords[dim],
                                             pixsize_inv)

                pixmax[dim] = find_max_pixel(pos[dim], hsml, minCoords[dim],
                                             pixsize_inv, max_pixel[dim])

            end

            bin_prefac = bin_q * m * rho_inv

            # second loop to calculate value
            @inbounds for i = pixmin[3]:pixmax[3]
                dz = param.z[i] - pos[3]

                @inbounds for j = pixmin[2]:pixmax[2]
                    dy = param.y[j] - pos[2]

                    @inbounds for k = pixmin[1]:pixmax[1]
                        dx = param.x[k] - pos[1]

                        # calculate simple distance to pixel center
                        distance_hsml = get_d_hsml_3D(dx, dy, dz, hsml_inv)

                        # update pixel value
                        image[k, j, i] += bin_prefac * kernel_value_3D(kernel, distance_hsml, hsml_inv)

                    end # end x-loop
                end # end y-loop
            end # end z-loop

        end # end check if in image

        # update for ProgressMeter
        if show_progress
            #lock(P_lock)  # uncomment to make thread-safe if needed in the future
            idx += 1
            ProgressMeter.update!(P, idx)
            #unlock(P_lock)
        end
    end

    if show_progress
        @info "Mapped $particles_in_image / $N particles."
    end

   return copy(transpose(image))
end

function get_bounds(pos::Array{Float32,2}, min::Float64, max::Float64, axis::Int64,
                    range_arr::Vector{Int64})

    k = findall(pos[range_arr,axis] .>= min  )
    m = findall(pos[range_arr[k],axis] .<= max  )

    return range_arr[k[m]]
end

function sphAdaptiveMapping(Pos, HSML, M, ρ, Bin_Quant; param::mappingParameters, kernel)

        """ Integral conserving SPH mapping. https://arxiv.org/abs/1803.03652
            !!! Diverges at the moment for large particle numbers due to lack
            of kernel normalisation !!!
        """

    N = length(M)
    pix_Vol = param.pixelSideLength^3
    boundary = false

    # println("mapping ", N, " particles")
    # println("Running on ", Threads.nthreads(), " Threads.")

    val = zeros(length(param.x), length(param.y),length(param.z))

    minCoords = [param.x[1], param.y[1], param.z[1]]

    #@showprogress 1 "Mapping..." for p = 1:N
    @threads for p = 1:N

        @inbounds pos = Pos[p,:]
        @inbounds hsml = HSML[p]
        @inbounds mass = M[p]
        @inbounds rho = ρ[p]
        @inbounds bin_quant = Bin_Quant[p]

        pixmin = Vector{Int64}(undef,3)
        pixmax = Vector{Int64}(undef,3)

        for dim = 1:3

            @inbounds pixmin[dim] = Int64( round((pos[dim] - hsml - minCoords[dim]) / param.pixelSideLength ))

            if pixmin[dim] < 1

                boundary = true
                pixmin[dim] = 1

            end

            pixmax[dim] = Int64( round((pos[dim] + hsml - minCoords[dim]) / param.pixelSideLength ))

            if dim == 1
                max = length(param.x)
            elseif dim == 2
                max = length(param.y)
            else
                max = length(param.z)
            end

            if pixmax[dim] > max
                boundary = true
                pixmax[dim] = max
            end

        end

        norm = 0.
        loop = true

        if boundary && hsml > 0.5 * param.pixelSideLength
            norm = 1.0
        else
            # closestPixel = Vector{Int64}(undef,3)
            # closestDistance = Inf

            for i = pixmin[1]:pixmax[1]
                for j = pixmin[2]:pixmax[2]
                    for k = pixmin[3]:pixmax[3]

                        @inbounds distance = sqrt.( (param.x[i] - pos[1])^2 +
                                                    (param.y[j] - pos[2])^2 +
                                                    (param.z[k] - pos[3])^2 )

                        # if distance < closestDistance
                        #     closestPixel = [i, j, k]            # ?
                        #     closestDistance = distance          # ?
                        # end

                        norm += kernel_value_3D(kernel, distance/hsml, Float64.(hsml))

                    end
                end
            end

            norm *= pix_Vol

            if norm < 0.001
                # Line 468 ???
                #val = mass / ( rho * pix_Vol)
            else
                norm = 1.0/norm
            end

        end

        if loop == true        # ??

            closestPixel = Vector{Int64}(undef,3)
            closestDistance = Inf

            for i = pixmin[1]:pixmax[1]
                for j = pixmin[2]:pixmax[2]
                    for k = pixmin[3]:pixmax[3]

                        @inbounds distance = sqrt.( (param.x[i] - pos[1])^2 +
                                          (param.y[j] - pos[2])^2 +
                                          (param.z[k] - pos[3])^2 )

                        if distance < closestDistance
                            closestPixel = [i, j, k]
                            closestDistance = distance
                        end

                        val[i,j,k] += bin_quant * mass / rho * kernel_value_3D(kernel, distance/hsml, Float64.(hsml))  * norm

                    end
                end
            end

        end

    end
    #
    # x = zeros(length(param.x)*length(param.y))
    # y = zeros(length(param.x)*length(param.y))
    # c = zeros(length(param.x)*length(param.y))
    # x = zeros(length(param.x_max), length(param.y_max))
    # y = zeros(length(param.x_max), length(param.y_max))
    # c = zeros(length(param.x), length(param.y))
     c = zeros(length(param.y), length(param.x))

    #cell_count = 1
    for i = 1:length(param.x)
        for j = 1:length(param.y)

            @inbounds c[i,j] = sum(val[i,j,:])
            #c[j,i] = sum(val[i,j,:]) # transpose for imshow()

            # x[cell_count] = param.x[i]
            # y[cell_count] = param.y[j]
            # c[cell_count] = sum(val[i,j,:])
            #
            # cell_count += 1

        end
    end

    return c

end

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
# code for interoperability with exciting

import DelimitedFiles: readdlm
import HDF5: h5open, readmmap
import LightXML: parse_file, root, attribute
import Random: MersenneTwister, shuffle
import Mmap

mutable struct BSEProblemExciting <: AbstractBSEProblem
    input_xml

    N_rs
    r_lattice
    r_cartesian
    a_mat
    Ω0_vol
    atoms

    N_ks
    k_bz
    b_mat
    BZ_vol

    E_v
    E_c
    u_v
    u_c

    N_k_diffs
    q_bz
    q_2bz
    q_2bz_ind
    q_2bz_shift
    w_hat
    gqmax

    pmat
    #TODO: add scissor shift here?
end

function BSEProblemExciting(N_core, N_v, N_c, N_ks, N_rs, path)
    input_xml = parse_file(path * "/input.xml")

    r_lattice, r_cartesian, a_mat, Ω0_vol, k_bz, b_mat, BZ_vol = read_r_points_k_points(N_rs, path)

    atoms = []
    for species in root(input_xml)["structure"][1]["species"]
        species_name = attribute(species, "speciesfile")[1:(end - 4)]
        mt_radius = parse(Float64, attribute(species, "rmt"))
        for atom in species["atom"]
            position = [parse(Float64, x) for x in split(attribute(atom, "coord"))]
            push!(atoms, (species = species_name, position = position, mt_radius = mt_radius))
        end
    end

    N_k = size(k_bz, 2)
    E_v, E_c, u_v, u_c = read_eigenvalues_eigenfunctions(N_core, N_v, N_c, N_k, N_rs, k_bz, path)

    N_ks, N_k_diffs, q_bz, q_2bz, q_2bz_ind, q_2bz_shift, w_hat = read_q_points_screenedcoulomb(N_ks, Ω0_vol, path)
    gqmax = parse(Float64, attribute(root(input_xml)["xs"][1], "gqmax"))

    pmat = read_pmat(N_core, N_v, N_c, N_k, path)

    return BSEProblemExciting(input_xml, N_rs, r_lattice, r_cartesian, a_mat, Ω0_vol, atoms, N_ks, k_bz, b_mat, BZ_vol, E_v, E_c, u_v, u_c, N_k_diffs, q_bz, q_2bz, q_2bz_ind, q_2bz_shift, w_hat, gqmax, pmat)
end

function read_r_points_k_points(N_rs, path) # TODO: make into type constructor?
    lattice_out = readdlm(path * "/LATTICE.OUT")
    a_mat = zeros(3, 3)
    a_mat .= lattice_out[8:10, 1:3]
    b_mat = zeros(3, 3)
    b_mat .= lattice_out[23:25, 1:3]

    Ω0_vol = abs(det(a_mat))
    BZ_vol = abs(det(b_mat))

    rl_1d_1 = range(0.0, stop = 1.0, length = N_rs[1] + 1)[1:N_rs[1]]
    rl_1d_2 = range(0.0, stop = 1.0, length = N_rs[2] + 1)[1:N_rs[2]]
    rl_1d_3 = range(0.0, stop = 1.0, length = N_rs[3] + 1)[1:N_rs[3]]
    rl_1 = kron(ones(N_rs[3]), ones(N_rs[2]), rl_1d_1)
    rl_2 = kron(ones(N_rs[3]), rl_1d_2, ones(N_rs[1]))
    rl_3 = kron(rl_1d_3, ones(N_rs[2]), ones(N_rs[1]))
    r_lattice = transpose(hcat(rl_1, rl_2, rl_3))
    r_cartesian = zero(r_lattice)
    for ir in 1:size(r_lattice, 2)
        r_cartesian[:, ir] = a_mat * r_lattice[:, ir]
    end

    k_points_out = readdlm(path * "/KPOINTS_QMT001.OUT")

    k_bz = float.(transpose(k_points_out[2:end, 2:4]))

    r_lattice, r_cartesian, a_mat, Ω0_vol, k_bz, b_mat, BZ_vol
end

function read_eigenvalues_eigenfunctions(N_core, N_v, N_c, N_k, N_rs, k_bz, path)
    N_r = prod(N_rs)

    eigval_out = readdlm(path * "/EIGVAL_QMT001.OUT")

    N_states = eigval_out[2, 1]
    eigval = zeros(N_states, N_k)
    for ik in 1:N_k
        eigval[:, ik] = eigval_out[(4 + 1 + (N_states + 2) * (ik - 1)):(4 + N_states + (N_states + 2) * (ik - 1)), 2]
    end

    E_v = eigval[(N_core + 1):(N_core + N_v), :]
    E_c = eigval[(N_core + N_v + 1):(N_core + N_v + N_c), :]

    mask_boundary = trues(N_rs[1] + 1, N_rs[2] + 1, N_rs[3] + 1)
    mask_boundary[N_rs[1] + 1, :, :] .= false
    mask_boundary[:, N_rs[2] + 1, :] .= false
    mask_boundary[:, :, N_rs[3] + 1] .= false
    mask_boundary = vec(mask_boundary)

    file = h5open(path * "/bse_output.h5", "r")

    s = open(path * "/u_v.bin", "w+")
    for ik in 1:N_k
        phase = kron(
            exp.(-2 * pi * im * range(0.0, stop = 1.0, length = N_rs[3] + 1)[1:N_rs[3]] .* k_bz[3, ik]),
            exp.(-2 * pi * im * range(0.0, stop = 1.0, length = N_rs[2] + 1)[1:N_rs[2]] .* k_bz[2, ik]),
            exp.(-2 * pi * im * range(0.0, stop = 1.0, length = N_rs[1] + 1)[1:N_rs[1]] .* k_bz[1, ik]))
        for iv in 1:N_v
            c = read(file, "wfplot/" * lpad(string(N_core + iv), 4, string(0)) * "/" * lpad(string(ik), 4, string(0)) * "/data")

            write(s, (c[1, mask_boundary] .+ im * c[2, mask_boundary]) .* phase)
        end
    end
    close(s)
    s_v = open(path * "/u_v.bin")
    u_v = Mmap.mmap(s_v, Array{Complex{Float64}, 3}, (N_r, N_v, N_k))

    s = open(path * "/u_c.bin", "w+")
    for ik in 1:N_k
        phase = kron(
            exp.(-2 * pi * im * range(0.0, stop = 1.0, length = N_rs[3] + 1)[1:N_rs[3]] .* k_bz[3, ik]),
            exp.(-2 * pi * im * range(0.0, stop = 1.0, length = N_rs[2] + 1)[1:N_rs[2]] .* k_bz[2, ik]),
            exp.(-2 * pi * im * range(0.0, stop = 1.0, length = N_rs[1] + 1)[1:N_rs[1]] .* k_bz[1, ik]))
        for ic in 1:N_c
            c = read(file, "wfplot/" * lpad(string(N_core + N_v + ic), 4, string(0)) * "/" * lpad(string(ik), 4, string(0)) * "/data")

            write(s, (c[1, mask_boundary] .+ im * c[2, mask_boundary]) .* phase)
        end
    end
    close(s)
    s_c = open(path * "/u_c.bin")
    u_c = Mmap.mmap(s_c, Array{Complex{Float64}, 3}, (N_r, N_c, N_k))

    close(file)

    return E_v, E_c, u_v, u_c
end

function read_q_points_screenedcoulomb(N_ks, Ω0_vol, path)
    q_points_out = readdlm(path * "/QPOINTS_SCR.OUT")
    q_bz = float.(transpose(q_points_out[2:end, 2:4]))

    N_k_diffs = size_q(N_ks)

    # q grid [0, 1[ ∪ [-1, 0[
    q_2bz = q_lattice(N_ks)

    N_k_diff = size(q_2bz, 2)

    q_2bz_ind = zeros(Int, N_k_diff)
    q_2bz_shift = zeros(Int, 3, N_k_diff)

    for iq in 1:N_k_diff
        q_2bz_ind[iq] = findfirst(vec(prod(mod.(q_bz .- q_2bz[:, iq] .+ 1e-5, 1) .< 1e-3; dims=1)))
        q_2bz_shift[:, iq] = round.(Int, q_bz[:, q_2bz_ind[iq]] - q_2bz[:, iq])
    end

    G_vec = Array{Int64, 2}[]
    for iq in 1:size(q_bz, 2)
        G = round.(Int, readdlm(path * "/GQPOINTS/GQPOINTS_SCR_Q" * lpad(string(iq), 5, string(0)) * ".OUT")[2:end, 2:4]' .- q_bz[:, iq])
        push!(G_vec, G)
    end

    w_hat = []
    file = h5open(path * "/bse_output.h5", "r")
    for iq in 1:size(q_bz, 2)
        wqq = read(file, "screenedpotential/" * lpad(string(iq), 4, string(0)) * "/wqq")
        G_len = size(G_vec[iq], 2)
        wqq_complex = wqq[1, 1:G_len, 1:G_len] + im * wqq[2, 1:G_len, 1:G_len]
        push!(w_hat, (wqq_complex, G_vec[iq])) # TODO: permutedim necessary?
    end
    close(file)

    return N_ks, N_k_diffs, q_bz, q_2bz, q_2bz_ind, q_2bz_shift, w_hat
end

function read_pmat(N_core, N_v, N_c, N_k, path)
    pmat = zeros(Complex{Float64}, N_v * N_c * N_k, 3)
    pmat_reshaped = reshape(pmat, N_v, N_c, N_k, 3)

    file = h5open(path * "/bse_output.h5", "r")
    for ik in 1:N_k
        pmat_block = permutedims(read(file, "pmat/" * lpad(string(ik), 4, string(0)) * "/pmat/"), [3, 4, 2, 1])
        pmat_reshaped[:, :, ik, :] = pmat_block[(N_core + 1):(N_core + N_v), (N_core + N_v + 1):(N_core + N_v + N_c), :, 1] -
                                im * pmat_block[(N_core + 1):(N_core + N_v), (N_core + N_v + 1):(N_core + N_v + N_c), :, 2]
    end
    close(file)

    return pmat
end

function size(prob::BSEProblemExciting)
    return (size(prob.E_v, 1), size(prob.E_c, 1), prod(prob.N_ks))
end

function size_k(prob::BSEProblemExciting)
    return prob.N_ks
end

function size_r(prob::BSEProblemExciting)
    return prob.N_rs
end

function energies(prob::BSEProblemExciting)
    return prob.E_v, prob.E_c
end

function orbitals(prob::BSEProblemExciting)
    return prob.u_v, prob.u_c
end

function lattice_matrix(prob::BSEProblemExciting)
    return prob.a_mat
end

function compute_v_hat(prob::BSEProblemExciting)
    return compute_v_hat(prob, prob.gqmax)
end

function compute_w_hat(prob::BSEProblemExciting)
    return prob.w_hat, prob.q_2bz_ind, prob.q_2bz_shift
end

function optical_absorption_vector(prob::BSEProblemExciting, direction)
    N_v, N_c, N_k = size(prob)
    ev = vec([prob.E_c[ic, ik] - prob.E_v[iv, ik] for iv in 1:N_v, ic in 1:N_c, ik in 1:N_k])
    return prob.pmat[:, direction] ./ ev
end

#TODO: maybe include this in the tests? currently this code is not used
function read_reference(prob, example_path)
    N_v, N_c, N_k = size(prob)
    f = h5open(example_path * "bse_matrix.h5")

    reordered_indices = vec(permutedims(reshape(1:(N_c * N_v * N_k), N_c, N_v, N_k), [2, 1, 3]))

    H_reference = read(f["H_BSE"])[1, reordered_indices, reordered_indices] - im * read(f["H_BSE"])[2, reordered_indices, reordered_indices]

    V_reference_reordered = zeros(Complex{Float64}, N_c * N_v * N_k, N_c * N_v * N_k)
    V_reshaped = reshape(V_reference_reordered, N_c * N_v, N_k, N_c * N_v, N_k)
    n = 1
    for ik in 1:N_k
        for jk in ik:N_k
            V_reshaped[:, ik, :, jk] = read(f, "0001/EXCLI_BSE-BAR_QMT001.OUT/" * lpad(string(n), 6, string(0)))[1, :, :] - im * read(f, "0001/EXCLI_BSE-BAR_QMT001.OUT/" * lpad(string(n), 6, string(0)))[2, :, :]
            n += 1
        end
    end
    V_reference = V_reference_reordered[reordered_indices, reordered_indices]
    for i in 1:size(V_reference, 1)
        V_reference[i, i] = real(V_reference[i, i])
    end

    W_reference_reordered = zeros(Complex{Float64}, N_c * N_v * N_k, N_c * N_v * N_k)
    W_reshaped = reshape(W_reference_reordered, N_c * N_v, N_k, N_c * N_v, N_k)
    n = 1
    for ik in 1:N_k
        for jk in ik:N_k
            W_reshaped[:, ik, :, jk] = read(f, "0001/SCCLI_QMT001.OUT/" * lpad(string(n), 6, string(0)))[1, :, :] - im * read(f, "0001/SCCLI_QMT001.OUT/" * lpad(string(n), 6, string(0)))[2, :, :]
            n += 1
        end
    end
    W_reference = W_reference_reordered[reordered_indices, reordered_indices]
    for i in 1:size(W_reference, 1)
        W_reference[i, i] = real(W_reference[i, i])
    end

    close(f)

    return H_reference, Hermitian(V_reference), Hermitian(W_reference)
end

ext: jl