text:
# BSE Hamiltonian
using FFTW, LinearMaps, ProgressMeter

# methods for setting up the matrix free Hamiltonian
"""
    setup_D(prob)

Assemble a sparse matrix of size `(N_v * N_c * N_K, N_v * N_c * N_k)`
for the diagonal part of the BSE Hamiltonian. The entries on the
diagonal are given by the single particle energy differences.
"""
function setup_D(prob)
    E_v, E_c = energies(prob)
    N_v, N_c, N_k = size(prob)

    D = spdiagm(0 => vec([E_c[ic, ik] - E_v[iv, ik] for iv in 1:N_v, ic in 1:N_c, ik in 1:N_k]))
    return D
end

"""
    setup_V(prob, isdf)

Create a linear operator for the application of `V` to a vector.
The result is a hermitian `LinearMap` of size
`(N_v * N_c * N_K, N_v * N_c * N_k)`.
"""
function setup_V(prob::AbstractBSEProblem, isdf::ISDF)
    N_rs = size_r(prob)
    N_v, N_c, N_k = size(prob)
    Ω0_vol = Ω0_volume(prob)
    N_μ = size(isdf)[3]

    V_tilde = assemble_V_tilde(prob, isdf)
    V_workspace = create_V_workspace(prob, isdf)

    u_v_vc_conj = conj.(interpolation_coefficients(isdf)[3])
    u_c_vc = interpolation_coefficients(isdf)[4]

    V = LinearMap{Complex{Float64}}(
        x -> V_times_vector(x, V_tilde, u_v_vc_conj, u_c_vc, V_workspace),
        N_v * N_c * N_k; ishermitian=true)

    return V
end

"""
    assemble_V_tilde(prob, isdf)

Assembles a Matrix representation of ``V`` in the basis given by interpolation vectors of the ISDF.
This corresponds to equation (2.33) in the paper.
Compared to (2.33), the scaling factor 1 / N_k is included here.
"""
function assemble_V_tilde(prob, isdf)
    N_rs = size_r(prob)
    N_r = prod(N_rs)
    N_k = size(prob)[3]
    Ω0_vol = Ω0_volume(prob)
    ζ_vc = interpolation_vectors(isdf)[3]

    v_hat = compute_v_hat(prob)

    ζ_vc_hat = zeros(Complex{Float64}, size(ζ_vc))
    for iμ in 1:size(ζ_vc, 2)
        ζ_vc_hat[:, iμ] = 1 / N_r * vec(fft(reshape(ζ_vc[:, iμ], N_rs)))
    end

    V_tilde = Ω0_vol / N_k * ζ_vc_hat' * (v_hat .* ζ_vc_hat)

    return V_tilde
end

"""
    create_V_workspace(N_v, N_c, N_k, N_μ)

Allocates momory for efficiently computing the application of ``V`` to
a vector. The output is a 5-tuple of arrays of different dimension and
sizes.
"""
function create_V_workspace(prob, isdf)
    N_v, N_c, N_k = size(prob)
    N_μ = size(isdf)[3]

    X = complex(zeros(N_v, N_c, N_k))
    B = complex(zeros(N_μ, N_c, N_k))
    C = complex(zeros(N_μ))
    D = complex(zeros(N_μ))
    E = complex(zeros(N_v, N_c, N_k))
    V_workspace = (X, B, C, D, E)

    return V_workspace
end

"""
    setup_W(prob, isdf)

Create a linear operator for the application of `W` to a vector.
The result is a hermitian `LinearMap` of size
`(N_v * N_c * N_K, N_v * N_c * N_k)`.
This corresponds to equation (2.35) in the paper. Instead of the
final expression in the equation, the expression in the third row is used.
"""
function setup_W(prob::AbstractBSEProblem, isdf)
    N_v, N_c, N_k = size(prob)
    u_v_vv_conj = conj.(interpolation_coefficients(isdf)[1]) #TODO: test whether the conjugation is really necessary for performance
    u_c_cc = interpolation_coefficients(isdf)[2]

    W_tilde_hat = assemble_W_tilde(prob, isdf)
    W_workspace = create_W_workspace(prob, isdf)

    W = LinearMap{Complex{Float64}}(
        x -> W_times_vector(x, W_tilde_hat, u_v_vv_conj, u_c_cc, W_workspace),
        N_v * N_c * N_k; ishermitian=true)

    return W
end

"""
    assemble_W_tilde(prob, isdf)

Assembles a Matrix representation of ``W`` in the basis given by interpolation vectors of the ISDF.
""" #TODO: better dispatch on problem type
function assemble_W_tilde(prob::AbstractBSEProblem, isdf::ISDF)
    w_hat, q_2bz_ind, q_2bz_shift = compute_w_hat(prob)

    N_rs = size_r(prob)
    N_k = size(prob)[3]
    N_qs = size_q(prob)
    Ω0_vol = Ω0_volume(prob)

    ζ_vv = interpolation_vectors(isdf)[1]
    ζ_cc = interpolation_vectors(isdf)[2]

    return assemble_W_tilde(w_hat, ζ_vv, ζ_cc, Ω0_vol, N_rs, N_k, N_qs, q_2bz_ind, q_2bz_shift)
end
function assemble_W_tilde(w_hat, ζ_vv, ζ_cc, Ω0_vol, N_rs, N_k, N_qs, q_2bz_ind, q_2bz_shift)
    N_r = size(ζ_vv, 1)
    N_ν = size(ζ_vv, 2)
    N_μ = size(ζ_cc, 2)
    N_q = prod(N_qs)

    ζ_vv_hat = zeros(Complex{Float64}, size(ζ_vv))
    for iμ in 1:size(ζ_vv, 2)
        ζ_vv_hat[:, iμ] = 1 / N_r * vec(fft(reshape(ζ_vv[:, iμ], N_rs)))
    end
    ζ_cc_hat = zeros(Complex{Float64}, size(ζ_cc))
    for iμ in 1:size(ζ_cc, 2)
        ζ_cc_hat[:, iμ] = 1 / N_r * vec(fft(reshape(ζ_cc[:, iμ], N_rs)))
    end

    W_tilde = complex(zeros(N_qs[1], N_qs[2], N_qs[3], N_μ, N_ν))
    W_tilde_reshaped = reshape(W_tilde, N_q, N_μ, N_ν)

    # TODO: can this be written in terms of more general data structures as W_tilde[iq, :, :] = 1 / (Ω0_vol^2 * N_k) * ζ_cc_hat' * w_hat[iq] * ζ_vv_hat?
    for iq in 1:N_q
        G_shift_indices = vec(mapslices(G -> G_vector_to_index(G, N_rs), w_hat[q_2bz_ind[iq]][2] .+ q_2bz_shift[:, iq]; dims=1))
        W_tilde_reshaped[iq, :, :] = Ω0_vol / N_k *
            ζ_cc_hat[G_shift_indices, :]' * w_hat[q_2bz_ind[iq]][1] * ζ_vv_hat[G_shift_indices, :]
    end

    return fft!(W_tilde, (1, 2, 3))
end

"""
    create_W_workspace(prob, isdf)

Allocates momory for efficiently computing the application of ``W`` to
a vector. The output is a 17-tuple of FFTW plans, arrays of different
dimension and sizes, and tuples.
""" #TODO: better dispatch on problem type
function create_W_workspace(prob::AbstractBSEProblem, isdf::ISDF)
    N_v, N_c, N_k = size(prob)
    N_ks = size_k(prob)
    N_k = prod(N_ks)
    Ω0_vol = Ω0_volume(prob)
    N_qs = size_q(prob)

    N_ν = size(isdf)[1]
    N_μ = size(isdf)[2]

    p = plan_fft(zeros(Complex{Float64}, N_qs); flags = FFTW.PATIENT, timelimit=60)
    p_back = plan_bfft(zeros(Complex{Float64}, N_qs); flags = FFTW.PATIENT, timelimit=60)
    X = zeros(Complex{Float64}, N_v, N_c, N_k)
    B = zeros(Complex{Float64}, N_ν, N_c, N_k)
    C = zeros(Complex{Float64}, N_k, N_μ, N_ν)
    D = zeros(Complex{Float64}, N_k, N_μ, N_ν)
    E = zeros(Complex{Float64}, N_ν, N_c, N_k)
    F = zeros(Complex{Float64}, N_v, N_c, N_k)
    c_reshaped = zeros(Complex{Float64}, N_ks..., Threads.nthreads())
    c_padded = zeros(Complex{Float64}, N_qs..., Threads.nthreads())
    c_transformed = zeros(Complex{Float64}, N_qs..., Threads.nthreads())
    d_reshaped = zeros(Complex{Float64}, N_ks..., Threads.nthreads())
    d_padded = zeros(Complex{Float64}, N_qs..., Threads.nthreads())
    w_tilde_reshaped = zeros(Complex{Float64}, N_qs..., Threads.nthreads())
    w_tilde_hat = zeros(Complex{Float64}, N_qs..., Threads.nthreads())
    W_workspace = (p, p_back, X, B, C, D, E, F, c_reshaped, c_padded, d_reshaped, d_padded, c_transformed, w_tilde_reshaped, w_tilde_hat, N_ks, N_qs)

    return W_workspace
end

# MIT licensed via julia
"""
    w_conv!(b, w, a, ap, bp, cp, w_hat, p, p_back)

Compute the convolution of `w` and `a` and store it in `b`. This means

``b_i = \\sum_j a_j w_{i - j + 1}``

with periodic `w`. The function uses Fourier transforms in the form
of FFTW forward and backward plans `p` and `p_back`. The array `a` is
zero padded to the size of `w`.
"""
function w_conv!(b, w, a, ap, bp, cp, w_hat, p, p_back)
    dim = ndims(a)
    sa = size(a)
    sw = size(w)

    indices = CartesianIndices(sa)

    copyto!(ap, indices, a, indices)

    mul!(cp, p, ap)
    # mul!(w_hat, p, w)
    cp .*= (1 / length(w)) .* w_hat

    mul!(bp, p_back, cp)

    copyto!(b, indices, bp, indices)

    return b
end

"""
    V_times_vector(x, V_tilde, u_v_vc_conj, u_c_vc, V_workspace)

Compute the product of ``V`` and `x`.
Corresponds to equation (4.5) in the paper.
"""
function V_times_vector(x, V_tilde, u_v_vc_conj, u_c_vc, V_workspace)
    N_μ = size(V_tilde, 1)
    N_k = size(u_v_vc_conj, 3)
    N_v = size(u_v_vc_conj, 2)
    N_c = size(u_c_vc, 2)
    X, B, C, D, E = V_workspace

    X[:] .= x
    @views for jk in 1:N_k
        mul!(B[:, :, jk], u_v_vc_conj[:, :, jk], X[:, :, jk])
    end
    for jμ in 1:N_μ
        C[jμ] = 0.0
        for jk in 1:N_k
            for jc in 1:N_c
                C[jμ] += u_c_vc[jμ, jc, jk] * B[jμ, jc, jk]
            end
        end
    end
    mul!(D, V_tilde, C)
    for ik in 1:N_k
        for ic in 1:N_c
            for iv in 1:N_v
                E[iv, ic, ik] = 0.
                for iμ in 1:N_μ
                    E[iv, ic, ik] += conj(u_c_vc[iμ, ic, ik] * u_v_vc_conj[iμ, iv, ik]) * D[iμ]
                end
            end
        end
    end

    return vec(E)
end

"""
    W_times_vector(x, W_tilde_hat, u_v_vv_conj, u_c_cc, W_workspace)

Compute the product of ``W`` and `x`.
Corresponds to equation (4.6) in the paper.
"""
function W_times_vector(x, W_tilde_hat, u_v_vv_conj, u_c_cc, W_workspace)
    N_v = size(u_v_vv_conj, 2)
    N_c = size(u_c_cc, 2)
    N_k = size(u_v_vv_conj, 3)
    N_ν = size(u_v_vv_conj, 1)
    N_μ = size(u_c_cc, 1)
    p, p_back, X, B, C, D, E, F, c_reshaped_threads, c_padded_threads, d_reshaped_threads, d_padded_threads, c_transformed_threads, w_tilde_reshaped_threads, w_tilde_hat_threads, N_ks, N_qs = W_workspace

    X[:] .= x
    Threads.@threads for jk in 1:N_k
        mul!(@view(B[:, :, jk]), @view(u_v_vv_conj[:, :, jk]), @view(X[:, :, jk]))
        mul!(@view(C[jk, :, :]), @view(u_c_cc[:, :, jk]), transpose(@view(B[:, :, jk])))
    end
    Threads.@threads for iν in 1:N_ν
        @views for iμ in 1:N_μ
            c_reshaped_threads[:, :, :, Threads.threadid()] .= reshape(C[:, iμ, iν], N_ks)
            w_tilde_hat_threads[:, :, :, Threads.threadid()] .= W_tilde_hat[:, :, :, iμ, iν]
            w_conv!(d_reshaped_threads[:, :, :, Threads.threadid()], w_tilde_reshaped_threads[:, :, :, Threads.threadid()], c_reshaped_threads[:, :, :, Threads.threadid()], c_padded_threads[:, :, :, Threads.threadid()], d_padded_threads[:, :, :, Threads.threadid()], c_transformed_threads[:, :, :, Threads.threadid()], w_tilde_hat_threads[:, :, :, Threads.threadid()], p, p_back)
            D[:, iμ, iν] .= vec(d_reshaped_threads[:, :, :, Threads.threadid()])
        end
    end
    Threads.@threads for ik in 1:N_k
        mul!(@view(E[:, :, ik]), transpose(@view(D[ik, :, :])), conj.(@view(u_c_cc[:, :, ik])))
        mul!(@view(F[:, :, ik]), adjoint(@view(u_v_vv_conj[:, :, ik])), @view(E[:, :, ik]))
    end

    return copy(vec(F))
end

"""
    setup_H(prob, isdf)

Create a linear operator for the application of ``H = D + 2 V - W``
to a vector. The result is a hermitian `LinearMap` of size
`(N_v * N_c * N_K, N_v * N_c * N_k)`.
"""
function setup_H(prob, isdf)
    D = setup_D(prob)
    V = setup_V(prob, isdf)
    W = setup_W(prob, isdf)
    return D + 2 * V - W # TODO: add scissor shift here?
end

###############################################################################

# Some methods which are only used for consitency checks in the tests.

"""
    V_entry(v_hat, iv, ic, ik, jv, jc, jk, u_v, u_c, Ω0_vol, N_rs, N_ks)

Compute the matrix entry in ``V`` corresponding to the multiindex
`(iv ic ik, jv jc jk)`.
"""
function V_entry(v_hat, iv, ic, ik, jv, jc, jk, u_v, u_c, Ω0_vol, N_rs, N_ks)
    N_r = prod(N_rs)
    N_k = prod(N_ks)

    U_vc_1_hat = 1 / N_r * vec(fft(reshape(u_c[:, ic, ik] .* conj.(u_v[:, iv, ik]), N_rs)))
    U_vc_2_hat = 1 / N_r * vec(fft(reshape(u_c[:, jc, jk] .* conj.(u_v[:, jv, jk]), N_rs)))

    v = Ω0_vol / N_k * U_vc_1_hat' * (v_hat .* U_vc_2_hat)

    return v
end

"""
    assemble_exact_V(prob)

Assemble the full matrix for ``V``. The entries are computed using the function `V_entry`.
This should only be used for very small problems.
"""
function assemble_exact_V(prob)
    N_v, N_c, N_k = size(prob)
    N_rs = size_r(prob)
    N_ks = size_k(prob)
    Ω0_vol = Ω0_volume(prob)
    u_v, u_c = orbitals(prob)

    v_hat = compute_v_hat(prob)

    V_reshaped = zeros(Complex{Float64}, N_v, N_c, N_k, N_v, N_c, N_k)
    @showprogress 1 "Assemble exact V ..." for jk in 1:N_k, jc in 1:N_c, jv in 1:N_v, ik in 1:N_k, ic in 1:N_c, iv in 1:N_v
        V_reshaped[iv, ic, ik, jv, jc, jk] =
            V_entry(v_hat, iv, ic, ik, jv, jc, jk, u_v, u_c, Ω0_vol, N_rs, N_ks)
    end
    V = reshape(V_reshaped, N_v * N_c * N_k, N_v * N_c * N_k)
    for i in 1:(N_v * N_c * N_k)
        V[i, i] = real(V[i, i])
    end

    return V
end

"""
    W_entry(w_hat, iv, ic, ik, jv, jc, jk, u_v, u_c, Ω0_vol, N_rs, ikkp2iq, q_2bz_ind, q_2bz_shift)

Compute the matrix entry in ``V`` corresponding to the multiindex
`(iv ic ik, jv jc jk)`. This method is not effient for large problems
since `ikkp2iq` is quadratic in `N_k`.
"""
function W_entry(w_hat, iv, ic, ik, jv, jc, jk, u_v, u_c, Ω0_vol, N_rs, ikkp2iq, q_2bz_ind, q_2bz_shift)
    N_r = prod(N_rs)
    N_k = size(ikkp2iq, 1)

    iq = ikkp2iq[ik, jk]

    U_c_hat = 1 / N_r * vec(fft(reshape(u_c[:, ic, ik] .* conj.(u_c[:, jc, jk]), N_rs)))
    U_v_hat = 1 / N_r * vec(fft(reshape(u_v[:, iv, ik] .* conj.(u_v[:, jv, jk]), N_rs)))

    G_shift_indices = vec(mapslices(G -> G_vector_to_index(G, N_rs), w_hat[q_2bz_ind[iq]][2] .+ q_2bz_shift[:, iq]; dims=1))
    w = Ω0_vol / N_k *
        U_c_hat[G_shift_indices]' * (w_hat[q_2bz_ind[iq]][1] * U_v_hat[G_shift_indices])

    return w
end

"""
    assemble_exact_W(prob)

Assemble the full matrix for ``W``. The entries are computed using the function `W_entry`.
This should only be used for very small problems.
"""
function assemble_exact_W(prob)
    N_v, N_c, N_k = size(prob)
    N_rs = size_r(prob)
    N_ks = size_k(prob)
    Ω0_vol = Ω0_volume(prob)
    u_v, u_c = orbitals(prob)

    w_hat, q_2bz_ind, q_2bz_shift = compute_w_hat(prob)
    ikkp2iq = ikkp2iq_matrix(k_lattice(prob), q_lattice(prob))

    W_reshaped = zeros(Complex{Float64}, N_v, N_c, N_k, N_v, N_c, N_k)
    @showprogress 1 "Assemble exact W ..." for jk in 1:N_k, jc in 1:N_c, jv in 1:N_v, ik in 1:N_k, ic in 1:N_c, iv in 1:N_v
        W_reshaped[iv, ic, ik, jv, jc, jk] =
            W_entry(w_hat, iv, ic, ik, jv, jc, jk, u_v, u_c, Ω0_vol, N_rs, ikkp2iq, q_2bz_ind, q_2bz_shift)
    end
    W = reshape(W_reshaped, N_v * N_c * N_k, N_v * N_c * N_k)
    for i in 1:(N_v * N_c * N_k)
        W[i, i] = real(W[i, i])
    end

    return W
end

"""
    ikkp2iq_matrix(k_bz, q_2bz)

Compute the all indices in `q_bz` corresponding to differences of points in `k_bz`. The resulting matrix is such that the entry `l = ikkp2iq[i, j]` satisfies `k_bz[i] - k_bz[j] == q_2bz[l]`.

This method is not effient for large problems since it is quadratic in `N_k`.
"""
function ikkp2iq_matrix(k_bz, q_2bz)
    N_k = size(k_bz, 2)
    ikkp2iq = zeros(Int, N_k, N_k)
    for ik in 1:N_k
        for jk in 1:N_k
            k_diff = k_bz[:, ik] - k_bz[:, jk]
            ikkp2iq[ik, jk] = findfirst(vec(prod(abs.(k_diff .- q_2bz) .< 1e-3; dims=1)))
        end
    end
    return ikkp2iq
end


text:
# Application Product Test
module APT
    include("../src/sdistance.jl")
    include("../src/environments.jl")
    include("../test/draw.jl")
    include("../src/SignedDistanceFunction.jl")
    # import .Sdistance: signedDistance2D
    using .SignedDistanceFunction
    import .Draw: draw, draw2x2
    import Plots
    using Plots

    """
        Ns: 分割数の入った配列
        _csv_datafiles: CSV Data
        circle_n: 閉曲線データの種類 {"multi", others}

        Ns[start]からNs[end]の分割数でsubdomainを分割し、
        signedDistance2D()で_csv_datafilesで与えられた閉曲線の符合付距離函数を生成、
        それらをcontoursでplot.
    """
    function plots_contours(Ns::Array, _csv_datafiles, circle_n::Union{String,Nothing} = nothing)
        plots = []
        for (i, N) in enumerate(Ns)
            _x = [i for i = -L:2*L/N:L] # len:N+1
            _y = [i for i = -L:2*L/N:L] # len:N+1
            _phi = signedDistance2D(_csv_datafiles, N, circle_n)
            push!(plots, plot(_x, _y, _phi, title = "N=$N", st = :contour,grid=false))
        end
        println(plots)
        # plot(plots..., size = (2300, 2500), layout = (3, 2))
        plot(plots..., size = (3500, 1500), layout = (2, 3),titlefont=font(18),grid=false)
        if circle_n == "multi"
            savefig("test/image/NsPlots_floodfill.png")
        else
            savefig("test/image/NsPlots_isinside.png")
        end
    end

    """
        Ns: 分割数の入った配列
        _csv_datafiles: CSV Data
        circle_n: 閉曲線データの種類 {"multi", others}

        Ns[start]からNs[end]の分割数でsubdomainを分割し、
        signedDistance2D()で_csv_datafilesで与えられた閉曲線の符合付距離函数を生成、
        それらをwireframeでplot.
    """
    function plots_wireframe(Ns::Array, _csv_datafiles, circle_n::Union{String,Nothing} = nothing)
        plots = []
        for (i, N) in enumerate(Ns)
            _x = [i for i = -L:2*L/N:L] # len:N+1
            _y = [i for i = -L:2*L/N:L] # len:N+1
            _phi = signedDistance2D(_csv_datafiles, N, circle_n)
            push!(plots, plot(_x, _y, _phi, title = "N=$N", st = :wireframe))
        end
        println(plots)
        plot(plots..., size = (2300, 2500), layout = (3, 2))
        if circle_n == "multi"
            savefig("test/image/NsPlots_floodfill_wireframe.png")
        else
            savefig("test/image/NsPlots_isinside_wireframe.png")
        end
    end

    """
    N: field splits
    _csv_datafile: CSV Data files
    circle_n: Kinds of circle data{"multi" or "single"}
    描画
    """
    function plot_for_debug(N::Int = 1000, _csv_datafile::String = "../test/mock_csv_data/interface.csv", circle_n::Union{String,Nothing} = nothing)
        csvfile_name = match(r"\./test/mock_csv_data/(.*)", _csv_datafile[1:end-4]).captures
        #===  case: double circle ===#
        if circle_n == "multi"
            _x = [i for i = -L:2*L/N:L] # len:N+1
            _y = [i for i = -L:2*L/N:L] # len:N+1
            _phi = signedDistance2D(_csv_datafile, N, "multi")
            tmpname = csvfile_name[1]
            if JULIA_MULTI_PROCESS
                filename = "$tmpname" * "_multicurves_multiprocess_" * "$(N)"
                draw(_x, _y, _phi, filename)
                # draw2x2(_x, _y, _phi, filename)
            else
                _filename = "$tmpname" * "_multicurves_normalprocess_" * "$(N)"
                draw(_x, _y, _phi, _filename)
            end

            return _phi

            #=== case: simple circle ===#
        else
            # scatter(_gamma[:,1], _gamma[:,2], markersize = 2)
            # savefig("test/image/the_data.png")
            _x = [i for i = -L:2*L/N:L] # len:N+1
            _y = [i for i = -L:2*L/N:L] # len:N+1
            _phi = signedDistance2D(_csv_datafile, N, "single")
            # DataFrame(_phi, :auto) |> CSV.write("./test/result/interface_result_tmp.csv", header = false)
            tmpname = csvfile_name[1]
            if JULIA_MULTI_PROCESS
                filename = "$tmpname" * "_jordancurve_multiprocess_" * "$(N)"
                draw(_x, _y, _phi, filename)
                # draw2x2(_x, _y, _phi, filename)
            else
                _filename = "$tmpname" * "_jordancurve_normalprocess_" * "$(N)"
                draw(_x, _y, _phi, _filename)
            end
            return _phi
        end
    end
    export plots_contours, plots_wireframe, plot_for_debug
end

ext: jl

max_stars_repo_name: jabelic/Signed-distance-function

max_issues_repo_name: jabelic/Signed-distance-function

max_forks_repo_name: jabelic/Signed-distance-function


text:
using FourierTools, BenchmarkTools

function main()
    x = randn((133, 513, 33))
    y = copy(x)
    @btime $y .= real.(ifft(fft($x)));
    @btime $y .= real.(ifft(ifftshift(fftshift(fft($x)))));
    @btime $y .= real.(iffts(ffts($x)));
    @btime $y .= real.(ift(ft($x)));

    return
end

main()


text:
#Heatmap plots

function PlotHeatmap(label_id, x::AbstractArray{T}, rows::Integer, cols::Integer, scale_min::Real = 0., scale_max::Real = 0., label_fmt::String = "&.1f", bounds_min = ImPlot_ImPlotPoint(0.0,0.0),
                 bounds_max = ImPlot_ImPlotPoint(1.0,1.0)) where {T<:ImPlotData}
     ImPlot_PlotHeatmap(label_id, x, rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max)
end

function PlotHeatmap(label_id, x::AbstractArray{T}, rows::Integer, cols::Integer, scale_min::Real = 0., scale_max::Real = 0., label_fmt::String = "&.1f", bounds_min = ImPlot_ImPlotPoint(0.0,0.0),
                 bounds_max = ImPlot_ImPlotPoint(1.0,1.0)) where {T<:Float64}
    ImPlot_PlotHeatmap(label_id, Float64.(x), rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max)
end

function PlotHeatmap(x::AbstractArray{T}, rows, cols, scale_min = 0.0, scale_max = 1.0;
                 label_id::String = "", label_fmt=C_NULL, bounds_min = ImPlot_ImPlotPoint(0.0,0.0),
                 bounds_max = ImPlot_ImPlotPoint(1.0,1.0)) where {T<:ImPlotData}

    ImPlot_PlotHeatmap(label_id, x, rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max)
end

heatmap_functions = [
    Int16 => LibCImGui.ImPlot_PlotHeatmap_S16Ptr,
    Int32 => LibCImGui.ImPlot_PlotHeatmap_S32Ptr,
    Int64 => LibCImGui.ImPlot_PlotHeatmap_S64Ptr,
    UInt16 => LibCImGui.ImPlot_PlotHeatmap_U16Ptr,
    UInt32 => LibCImGui.ImPlot_PlotHeatmap_U32Ptr,
    UInt64 => LibCImGui.ImPlot_PlotHeatmap_U64Ptr,
    Cfloat => LibCImGui.ImPlot_PlotHeatmap_FloatPtr,
    Cdouble => LibCImGui.ImPlot_PlotHeatmap_doublePtr,
]
for (type, func) in heatmap_functions
    @eval ImPlot_PlotHeatmap(label_id, x::AbstractArray{<:$type}, rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max) = $func(label_id, x, rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max)
end


text:
# We wish to extend operations to identically named methods dispatched
# on `Machine`s. For example, we have from the model API
#
# `predict(model::M, fitresult, X) where M<:Supervised`
#
# but want also want to define
#
# 1. `predict(machine::Machine, X)` where `X` is concrete data
#
# and we would like the syntactic sugar (for `X` a node):
#
# 2. `predict(machine::Machine, X::Node) = node(predict, machine, X)`
#
# Finally, for a `model` that is `ProbabilisticComposite`,
# `DetermisiticComposite`, or `UnsupervisedComposite`, we want
#
# 3. `predict(model, fitresult, X) = fitresult.predict(X)`
#
# which makes sense because `fitresult` in those cases is a named
# tuple keyed on supported operations and with nodes as values.

## TODO: need to add checks on the arguments of
## predict(::Machine, ) and transform(::Machine, )

_err_rows_not_allowed() =
    throw(ArgumentError("Calling `transform(mach, rows=...)` when "*
                        "`mach.model isa Static` is not allowed, as no data "*
                        "is bound to `mach` in this case. Specify a explicit "*
                        "data or node, as in `transform(mach, X)`, or "*
                        "`transform(mach, X1, X2, ...)`. "))
_err_serialized(operation) =
    throw(ArgumentError("Calling $operation on a "*
                        "deserialized machine with no data "*
                        "bound to it. "))

warn_serializable_mach(operation) = "The operation $operation has been called on a "*
                        "deserialised machine mach whose learned parameters "*
                        "may be unusable. To be sure, first run restore!(mach)."

# 0. operations on machine, given rows=...:

for operation in OPERATIONS

    if operation != :inverse_transform

        ex = quote
            function $(operation)(mach::Machine{<:Model,false}; rows=:)
                # catch deserialized machine with no data:
                isempty(mach.args) && _err_serialized($operation)
                return ($operation)(mach, mach.args[1](rows=rows))
            end
            function $(operation)(mach::Machine{<:Model,true}; rows=:)
                # catch deserialized machine with no data:
                isempty(mach.args) && _err_serialized($operation)
                model = mach.model
                return ($operation)(model,
                                    mach.fitresult,
                                    selectrows(model, rows, mach.data[1])...)
            end
        end
        eval(ex)

    end
end

# special case of Static models (no training arguments):
transform(mach::Machine{<:Static}; rows=:) = _err_rows_not_allowed()

inverse_transform(mach::Machine; rows=:) =
            throw(ArgumentError("`inverse_transform()(mach)` and "*
                                "`inverse_transform(mach, rows=...)` are "*
                                "not supported. Data or nodes "*
                                "must be explictly specified, "*
                                "as in `inverse_transform(mach, X)`. "))

_symbol(f) = Base.Core.Typeof(f).name.mt.name

for operation in OPERATIONS

    ex = quote
        # 1. operations on machines, given *concrete* data:
        function $operation(mach::Machine, Xraw)
            if mach.state != 0
                mach.state == -1 && @warn warn_serializable_mach($operation)
                return $(operation)(mach.model,
                                    mach.fitresult,
                                    reformat(mach.model, Xraw)...)
            else
                error("$mach has not been trained.")
            end
        end

        function $operation(mach::Machine{<:Static}, Xraw, Xraw_more...)
            isdefined(mach, :fitresult) || (mach.fitresult = nothing)
            return $(operation)(mach.model, mach.fitresult,
                                    Xraw, Xraw_more...)
        end

        # 2. operations on machines, given *dynamic* data (nodes):
        $operation(mach::Machine, X::AbstractNode) =
            node($(operation), mach, X)

        $operation(mach::Machine{<:Static},
                   X::AbstractNode,
                   Xmore::AbstractNode...) =
                       node($(operation), mach, X, Xmore...)
    end
    eval(ex)
end

## SURROGATE AND COMPOSITE MODELS

const err_unsupported_operation(operation) = ErrorException(
    "The `$operation` operation has been applied to a composite model or learning "*
    "network machine that does not support it. "
)

for operation in [:predict,
                  :predict_joint,
                  :transform,
                  :inverse_transform]
    ex = quote
        function $operation(model::Union{Composite,Surrogate}, fitresult,X)
            if hasproperty(fitresult, $(QuoteNode(operation)))
                return fitresult.$operation(X)
            else
                throw(err_unsupported_operation($operation))
            end
        end
    end
    eval(ex)
end

for (operation, fallback) in [(:predict_mode, :mode),
                              (:predict_mean, :mean),
                              (:predict_median, :median)]
    ex = quote
        function $(operation)(m::Union{ProbabilisticComposite,ProbabilisticSurrogate},
                              fitresult,
                              Xnew)
            if hasproperty(fitresult, $(QuoteNode(operation)))
                return fitresult.$(operation)(Xnew)
            end
            return $(fallback).(predict(m, fitresult, Xnew))
        end
    end
    eval(ex)
end


text:
using DataFrames, Multisets

export get_vertex_counts

"""
Takes a GroupedDataFrame and a string indicating column name for vertex labels,
outputs a Vector{Multiset}, that is, a multisets of vertex counts for each group
in the data.
"""
function get_vertex_counts(
    X::GroupedDataFrame;
    vertex_col_name::String="page"
    )

    z = Vector{Multiset}()
    for group in X
        push!(z, Multiset(group[:, vertex_col_name]))
    end
    return z
end

#
# get_vertex_counts(df_g)

# ## Using

# function get_vertex_counts(
#     X::GroupedDataFrame;
#     vertex_col_name::String="page"
#     )
#
#     z = Vector{Multiset}()
#     for group in X
#         push!(z, countmap(group[:, vertex_col_name]))
#     end
#     return z
# end

# function hamming(
#     X::Dict{T,Int} where T,
#     Y::Dict{T,Int} where T
#     )
#     z = 0.0
#     for k in intersect(keys(X), keys(Y))
#         z += abs(X[k]-Y[k])
#     end
#     for k in setdiff(keys(X), keys(Y))
#         z += X[k]
#     end
#     for k in setdiff(keys(Y), keys(X))
#         z += Y[k]
#     end
#     return z
# end
#
#
# function sum_values(d::Dict{T, Int} where T)
#     z = 0
#     for val in values(d)
#         z += val
#     end
#     return z
# end
#
# sum_values(tmp1)
#
# # Jaccard distance, this was cheaper than using steinhaus on hamming
# function jaccard(
#     X::Dict{T, Int} where T,
#     Y::Dict{T, Int} where T
#     )
#
#     a = 0.0
#     for k in intersect(keys(X), keys(Y))
#         a += min(X[k], Y[k])
#     end
#     b = sum_values(mergewith(max, X, Y))
#     return 1-a/b
# end
#
#
# d1 = Dict("a"=>10, "b"=>2)
# d2 = Dict("a"=>1, "d"=>20)
#
# using BenchmarkTools
#
# @time hamming(tmp1, tmp2)
#
# @btime jaccard(tmp1, tmp2)
#
# ## Distance Matrices
#
#
# @time
#
# tmp = get_dist_mat(data)
#
# using UMAP, DataFrames, StatsPlots
#
# Random.seed!(1234);  # Seed for repoduceability
#
# min_dist = 0.001
# n_neighbors = 50
# embedding = umap(tmp;
#     metric=:precomputed,
#     min_dist=min_dist,
#     n_neighbors=n_neighbors)
#
# embedding = DataFrame(embedding')  # Make dataframe
#
# @df embedding scatter(
#     :x1, :x2,
#     ms=2.3, markerstrokewidth=0, markeralpha=0.4,
#     leg=false,
#     size = (800,600))
#
#
# multiset(

max_stars_repo_path: src/Data Processing/VertexCounts.jl

max_issues_repo_path: src/Data Processing/VertexCounts.jl

max_forks_repo_path: src/Data Processing/VertexCounts.jl


text:
# ------------------------------------------------------------------------------------------
# # Distributed Data
#
# JuliaDB can distribute datasets across multiple processes and can even work with data
# larger than available memory (RAM).
#
# First let's start by adding some worker Julia processes.  If you do not specify a number,
# `addprocs` will add as many workers as there are CPU cores available.
# ------------------------------------------------------------------------------------------

addprocs(2)

using JuliaDB

# print table rather than column summary
IndexedTables.set_show_compact!(false);

# ------------------------------------------------------------------------------------------
# ### When multiple processes are available, `loadtable` will create distributed tables.
#
# - Note the line above the column names printed below:
#
# `Distributed Table with 56023 rows in 2 chunks`
# ------------------------------------------------------------------------------------------

dt = loadtable("stocksample", filenamecol = :Ticker, indexcols = [:Ticker, :Date])

# ------------------------------------------------------------------------------------------
# ### Notable difference 1: No `getindex`
# ------------------------------------------------------------------------------------------

# We are no longer able to get the row based on row number

dt[1]

# ------------------------------------------------------------------------------------------
# ### Notable difference 2: Not iterable
# ------------------------------------------------------------------------------------------

# Similarly, we're no longer allowed to iterate over the rows

for row in dt
    println(row)
end

# ------------------------------------------------------------------------------------------
# # Bring Distributed Table Into Master Process
#
# - While not necessary for most operations, you may occasionally want to bring a dataset
# into the master process.
# - This is accomplished with the `collect` function.
# - Note that after `collect`ing, the header says `Table` instead of `Distributed Table`.
# ------------------------------------------------------------------------------------------

t = collect(dt)

# ------------------------------------------------------------------------------------------
# # Table Operations Still Work on a Distributed Table!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Functions that return a single value still return a single value.
# ------------------------------------------------------------------------------------------

reduce(+, dt; select = :Close)

# ------------------------------------------------------------------------------------------
# ### Functions that returned a Table now return a Distributed Table.
# ------------------------------------------------------------------------------------------

groupreduce(+, dt, :Ticker; select = :Close)

# ------------------------------------------------------------------------------------------
# ### Functions that returned an `Array` now return a `DArray`.
# ------------------------------------------------------------------------------------------

select(dt, :Close)

# ------------------------------------------------------------------------------------------
# - Similar to tables, `collect` will change a `DArray` to an `Array`:
# ------------------------------------------------------------------------------------------

collect(select(dt, :Close))

# ------------------------------------------------------------------------------------------
# # Out-of-Core Functionality
#
# - Doc Reference: http://juliadb.org/latest/manual/out-of-core.html
#
# ### JuliaDB can be used to load/query datasets that are too big to fit in memory (RAM).
#
# - The current out-of-core design is based around working with many files where each fits
# comfortably in memory.
# - In the future, JuliaDB will support individual files that are larger than memory.
#
# #### Out-of-core support is restricted to `reduce`, `groupreduce`, and (a subset of)
# `join`
#
# - Now let's pretend the stock data is much larger than it actually is and go through the
# steps of how to handle a larger-than-memory dataset.
#
# ### 1. Data is loaded into a distributed dataset in *chunks* that fit in memory
#
# - A *chunk* contains the data from at least one file (one file cannot be split into
# multiple chunks).
# - The `loadtable` and `loadndsparse` functions have the keyword arguments:
#   - `output`: Directory where the loaded data is written to an efficient binary format.
#   - `chunks`: Number of parts to split the data into
# - Data processing will occur on `number_of_files / chunks` files at a time.
#
#
# - Here we load the stock data in 8 chunks (one file per chunk) into the output directory
# `bin`.
# ------------------------------------------------------------------------------------------

loadtable("stocksample", output = "bin", chunks=8,
    filenamecol=:Ticker, indexcols=[:Ticker, :Date])

# ------------------------------------------------------------------------------------------
# ### 2. Data can then be `load`ed from the binary format
#
# - `bintable` is now a distributed table made up of chunks that are on disk.
# ------------------------------------------------------------------------------------------

bintable = load("bin")

# ------------------------------------------------------------------------------------------
# ### 3. We can use out-of-core supported functions like usual
#
# - When we call a function like `groupreduce`, each worker process will operate on one
# chunk at a time.
#
# - You want to be sure that both `nprocs() * (size of chunks)` and the output of
# `reduce`/`groupreduce` fits in memory!
# ------------------------------------------------------------------------------------------

groupreduce(+, bintable, :Ticker; select = :Volume)

# ------------------------------------------------------------------------------------------
# # Plotting Big Data
#
# - There are a few options for plotting JuliaDB datasets through the use of the Plots
# package.
#
# ## Plotting with StatPlots
#
# - The StatPlots package allows you to plot a variety of tabular data structures, including
# those in JuliaDB.
# - It gives you the entire power and flexibility of Plots, but **not for Distributed
# Tables**.
# ------------------------------------------------------------------------------------------

using StatPlots
gr()

@df collect(dt) plot(:Date, :Close; group=:Ticker, legend=:topleft)

# ------------------------------------------------------------------------------------------
# ## Plotting with `partitionplot`
#
# - For large datasets, it may not be feasible to plot every point.
# - JuliaDB offers an alternative with `partitionplot`, which uses types in OnlineStats to
# incrementally build summaries of the data that can be plotted for datasets of **any**
# size.
# - The motivating idea is that the human brain cannot comprehend a huge number of points,
# so we can portray nearly the same information with summaries of the data rather than the
# data itself.
#
#
# - The syntax for `partitionplot` is
#
#   `partitionplot(table[, x], y; by, nparts = 100, stat = Extrema(), dropmissing = false)`
#
#   where:
#
#   - `x`, `y`, and `by` are selections
#   - The number of summaries to be calculated are between `nparts` and `2 * nparts`
#   - The type of summary is `stat`
#   - Missing values are skipped if `dropmissing == true`
#   - Any additional keyword arguments will be passed to Plots
#
# - In the following example, we "recreate" the above plot by plotting the extrema (maximum
# and minimum) of the closing price over sections of time.
# ------------------------------------------------------------------------------------------

using OnlineStats

partitionplot(collect(bintable), :Date, :Close, by = :Ticker, nparts=100, stat = Extrema(),
    legend = :topleft)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
abstract type Piece end

abstract type AbstractPlayer end
abstract type AbstractGame end

mutable struct Pawn <: Piece
  player::AbstractPlayer
  row::Int
  col::Int
end

mutable struct AntiPawn <: Piece
  player::AbstractPlayer
  row::Int
  col::Int
end

mutable struct Bishop <: Piece
  player::AbstractPlayer
  row::Int
  col::Int
end

mutable struct Queen <: Piece
  player::AbstractPlayer
  row::Int
  col::Int
end

mutable struct Rook <: Piece
  player::AbstractPlayer
  row::Int
  col::Int
end

mutable struct King <: Piece
  player::AbstractPlayer
  row::Int
  col::Int
end

mutable struct Knight <: Piece
  player::AbstractPlayer
  row::Int
  col::Int
end

mutable struct Game <: AbstractGame
  window::Window
  players::Vector{AbstractPlayer}
  board::Matrix{Union{Nothing,Piece}}
  is_whites_turn::Bool
  history::Vector{AbstractString}
  halfmove_clock::Int
  fullmove_number::Int
end

mutable struct Player <: AbstractPlayer
  is_white::Bool
  pieces::Vector{Piece}
  game::AbstractGame
  can_castle_left::Bool
  can_castle_right::Bool
  has_castled::Bool
  anti_pawn::Union{Nothing,AntiPawn}
  is_human::Bool
end

function Player(is_white,game,is_human)
  Player(is_white,[],game,true,true,false,nothing,is_human)
end

function Game(window::Window; player_1_is_human=true, player_2_is_human=false)

  players = AbstractPlayer[]
  board = Matrix{Union{Nothing,Piece}}(nothing,8,8)

  game = Game(
    window,players,board,true,
    [], 0, 1
  )

  white = Player(true,game,player_1_is_human)
  black = Player(false,game,player_2_is_human)

  push!(game.players, white)
  push!(game.players, black)

  white_i = 2
  black_i = 7

  for j = 1:8
    push!(white.pieces, Pawn(white,white_i,j))
    push!(black.pieces, Pawn(black,black_i,j))
  end

  white_i = 1
  black_i = 8
  for j = [1,8]
    push!(white.pieces, Rook(white,white_i,j))
    push!(black.pieces, Rook(black,black_i,j))
  end

  white_i = 1
  black_i = 8
  for j = [2,7]
    push!(white.pieces, Knight(white,white_i,j))
    push!(black.pieces, Knight(black,black_i,j))
  end

  white_i = 1
  black_i = 8
  for j = [3,6]
    push!(white.pieces, Bishop(white,white_i,j))
    push!(black.pieces, Bishop(black,black_i,j))
  end

  j = 4
  push!(white.pieces, Queen(white,white_i,j))
  push!(black.pieces, Queen(black,black_i,j))

  j = 5
  push!(white.pieces, King(white,white_i,j))
  push!(black.pieces, King(black,black_i,j))

  for piece in [white.pieces..., black.pieces...]
    board[piece.row,piece.col] = piece
  end

  game

end


text:
import Base.CoreLogging:
    AbstractLogger,
    LogLevel,
    handle_message,
    min_enabled_level,
    shouldlog,
    global_logger

struct BaseLogger <: AbstractLogger
    min_level::LogLevel
end

min_enabled_level(logger::BaseLogger) = logger.min_level
shouldlog(logger::BaseLogger, args...) = true

function handle_message(::BaseLogger, cl_level, msg, mod, group, id, filepath, line; kwargs...)
    logger = getlogger(mod)
    level = lowercase(string(cl_level))
    log(logger, logger.record(logger.name, level, getlevels(logger)[level], msg))
end

function substitute!(level::LogLevel=min_enabled_level(global_logger()))
    global_logger(BaseLogger(level))
    notice(getlogger(@__MODULE__), "Substituting global logging with Memento")
end


text:
# Converts terms in the expression tree to .nl supported expressions
convert_formula(c) = c
convert_formula(c::LinearityExpr) = convert_formula(c.c)

function convert_formula(c::Expr)
    map!(convert_formula, c.args, c.args)

    if c.head == :comparison
        n = length(c.args)
        @assert isodd(n)

        @assert(n > 3)

        # If more than binary comparison, we need to chain them together
        # Get binary comparisons in sequence and chain with nested &&
        # It looks like we might be able to use an n-ary &&, but that's not how
        # Julia parses chained && expressions.
        if n > 3
            new_expr = extract_binary_comparison(c, 1)
            for i in 3:2:(n - 2)
                new_expr = Expr(:&&, new_expr, extract_binary_comparison(c, i))
            end
            c = new_expr
        end

    elseif c.head == :call
        op = c.args[1]
        # Distinguish unary, binary and n-ary plus
        if op == :+
            if length(c.args) > 3
                # n-ary plus to sum
                c.args[1] = :sum
            elseif length(c.args) == 2
                # Remove unary plus
                c = c.args[2]
            end
        # Distinguish unary and binary minus
        elseif op == :-
            n = length(c.args)
            if n == 2
                if c.args[2] == 0
                    c = 0
                else
                    c.args[1] = :neg
                end
            end
        # .nl has no n-ary multiplication so we need to convert to binary
        elseif op == :*
            # Loop for each arg after the 3rd (more than binary)
            for _ in 1:(length(c.args) - 3)
                # Combine last term with previous to form a binary * expression
                arg = pop!(c.args)
                c.args[end] = :($(c.args[end]) * $arg)
            end
        # Handle normal conversion cases
        else
            if length(c.args) == 2
                c = get(unary_special_cases, c.args[1],
                        (x) -> Expr(:call, c.args[1], x))(c.args[2])
            end
        end
    end
    c
end

# Extracts `expression relation expression` from a larger comparison expression
function extract_binary_comparison(c::Expr, start::Integer)
    @assert c.head == :comparison
    @assert start <= length(c.args) - 2

    Expr(:call, c.args[start + 1], c.args[start], c.args[start + 2])
end

const unary_special_cases = Dict(
:cbrt  => (x) -> :($x ^ (1 / 3)),
:abs2  => (x) -> :($x ^ 2),
:inv   => (x) -> :(1 / $x),
:log2  => (x) -> :(log($x) / log(2)),
:log1p => (x) -> :(log(1 + $x)),
:exp2  => (x) -> :(2 ^ $x),
:expm1 => (x) -> :(exp($x) - 1),

:sec   => (x) -> :(1 / cos($x)),
:csc   => (x) -> :(1 / sin($x)),
:cot   => (x) -> :(1 / tan($x)),

:asec  => (x) -> :(acos(1 / $x)),
:acsc  => (x) -> :(asin(1 / $x)),
:acot  => (x) -> :(pi / 2 - atan($x)),

:sind  => (x) -> :(sin(pi /180 * $x)),
:cosd  => (x) -> :(cos(pi /180 * $x)),
:tand  => (x) -> :(tan(pi /180 * $x)),
:secd  => (x) -> :(1 / cos(pi / 180 * $x)),
:cscd  => (x) -> :(1 / sin(pi / 180 * $x)),
:cotd  => (x) -> :(1 / tan(pi / 180 * $x)),

:asind => (x) -> :(asin($x) * 180 / pi),
:acosd => (x) -> :(acos($x) * 180 / pi),
:atand => (x) -> :(atan($x) * 180 / pi),
:asecd => (x) -> :(acos(1 / $x) * 180 / pi),
:acscd => (x) -> :(asin(1 / $x) * 180 / pi),
:acotd => (x) -> :((pi / 2 - atan($x)) * 180 / pi),

:sech  => (x) -> :(1 / cosh($x)),
:csch  => (x) -> :(1 / sinh($x)),
:coth  => (x) -> :(1 / tanh($x)),

:asech => (x) -> :(acosh(1 / $x)),
:acsch => (x) -> :(asinh(1 / $x)),
:acoth => (x) -> :(atanh(1 / $x)),
)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "XPA"
version = v"2.1.20"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/ericmandel/xpa.git", "923cc1bc7e761424b87049b1a20853eefe921388"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/xpa
if [[ "${target}" == *-freebsd* ]]; then
    export CFLAGS=-fPIC
fi
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --enable-shared=yes
make -j$(nproc)
if [[ ${target} == *mingw* ]]; then
    make mingw-dll
    mkdir -p ${libdir}
    cp libxpa.dll ${libdir}/libxpa.dll
fi
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("xpamb", :xpamb),
    ExecutableProduct("xpaget", :xpaget),
    ExecutableProduct("xpainfo", :xpainfo),
    ExecutableProduct("xpans", :xpans),
    ExecutableProduct("xpaset", :xpaset),
    ExecutableProduct("xpaaccess", :xpaaccess),
    LibraryProduct("libxpa", :libxpa)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
# NOTE objects don't have the overwrite property. so with or witout ! they don't
# overwrite the axes they get placed on.

"""
    text(text, anchor)

### Example

```julia
text(t"x^{\\star}", (1, 1.5); fontsize=15)
```
"""
function text(text::String, anchor::Tuple; axes=nothing, o...)
    axes = check_axes(axes)
    t = Text2D(text=text, anchor=fl(anchor))
    set_properties!(t; defer_preview=true, o...)
    push!(axes.objects, t)
    return PreviewFigure(gcf())
end

"""
    _line(...)

Internal function to draw a horizontal or vertical line from axis to axis. See [`vline`](@ref)
and [`hline`](@ref).
"""
function _line(anchor::Float64, horiz::Bool; axes=nothing, o...)
    axes = check_axes(axes)
    l = StraightLine2D(anchor=anchor, horiz=horiz)
    set_properties!(l; defer_preview=true, o...)
    push!(axes.objects, l)
    return PreviewFigure(gcf())
end

"""
    vline(anchor; opts...)

Add a vertical line at a given point from the xaxis to the x2axis. Note that this does not adapt
the xaxis limits so if the point specified is outside of the current xaxis limits you won't see
anything (until you adjust the limits accordingly).

### Example

```julia
vline(0; color="red")
```
"""
vline(anchor::Real; o...) = _line(fl(anchor), false; o...)

"""
    hline(anchor; opts...)

Add a horizontal line at a given point from the yaxis to the y2axis. Note that this does not adapt
the axis limits so if the point specified is outside of the current axis limits you won't see
anything (until you adjust the limits accordingly).

### Example

```julia
hline(0; color="red")
```
"""
hline(anchor::Real; o...) = _line(fl(anchor), true; o...)

"""
    line(from, to; options)

Add a line from `from` (in `[x, y]` format) to `to` (same format).

### Example

```julia
line([0, 0], [1, 1]; ls="--")
line((0, 0), (1, 1)) # also valid
```
"""
line(a::Union{AVR,T2R}, b::Union{AVR,T2R}; o...) = plot!([a[1],b[1]], [a[2],b[2]]; o...)

"""
    box(anchor; opts...)

Add a box at position `anchor`.
### Example

```julia
text(t"x^{\\star}", (1, 1.5); fontsize=15)
```
"""
function box(size::Tuple, anchor::Tuple; axes=nothing, o...)
    axes = check_axes(axes)
    b = Box2D(anchor=fl(anchor), size=fl(size))
    set_properties!(b; defer_preview=true, o...)
    push!(axes.objects, b)
    return PreviewFigure(gcf())
end

"""
    colorbar(zmin, zmax, cmap; o...)

Add a colorbar on a side of the axes corresponding to the span `[zmin,zmax]` with a given colormap
`cmap` (vector of colors).
"""
function colorbar(zmin::Real, zmax::Real, cmap::Vector{<:Color}; axes=nothing, o...)
    axes = check_axes(axes)
    # check if there is an existing colorbar, if so remove it
    mask = isa.(axes.objects, Colorbar)
    b = Colorbar(fl(zmin), fl(zmax), cmap)
    # axes adjustment
    if b.position ∈ ["left","right"] # horizontal
        axes.scale = "0.8 auto"
    else
        axes.scale = "auto 0.8"
    end
    set_properties!(b; defer_preview=true, o...)
    if any(mask)
        axes.objects[findfirst(mask)] = b
    else
        push!(axes.objects, b)
    end
    return PreviewFigure(gcf())
end

"""
    colorbar(h; o...)

Adds a colorbar to a drawing handle generated by a heatmap.

### Example

```julia
h = heatmap(randn(20, 20))
colorbar(h)
```
"""
colorbar(dh::DrawingHandle{Heatmap}; o...) =
    colorbar(dh.drawing.zmin, dh.drawing.zmax, dh.drawing.cmap; o...)


text:
module TestBenchmarks
    using Benchmark

    benchmark(() -> 10^2, "Math", "Squaring", 100)
    benchmark(() -> 10^2, "Squaring", 100)

    benchmarks(Vector[
        [() -> 10^2, "Math", "Squaring", 100],
        [() -> 10 - 10, "Math", "Subtraction", 100]])
end


text:
using ProgressMeter
using CUDA: allowscalar
using Flux: @epochs, onecold
using Statistics: mean
using BSON: @save, @load

# FIXME dir structure?
include("train_utils.jl")

function threshold(W, ε, one=false)
    res = copy(W)
    res[abs.(res) .< ε] .= 0
    # FIXME this .= is scalar operation, must be run on CPU
    # a = copy(gpu(randn(20)))
    # a[abs.(a) .< 0.2] .= 0.3
    #
    # threshold(gpu(randn(20)), 0.2)
    if one
        res[abs.(res) .>= ε] .= 1
    end
    res
end

# sup_graph_metrics([1 1 0 1], [1 0 0 1])

function sup_graph_metrics(out, y)
    # TODO try enforce some of this metrics in loss?
    # d = convert(Int, sqrt(size(y,1)))
    d = size(y,1)
    nbatch = size(y, 3)
    # FIXME threshold value
    #
    # FIXME UPDATE this should be 0.5 for binary sigmoid xent results. But I
    # need to monitor what are really the values out there
    mat(y) = threshold(reshape(y, d, d, :), 0.5, true)

    mout = mat(out)
    my = mat(y)
    nnz = sum(mout .!= 0) / nbatch
    nny = sum(my .!= 0) / nbatch

    # tp = sum(mout[mout .== my] .== 1)
    # fp = sum(mout[mout .== 1] .!= my[mout .== 1])
    tp = sum(mout[my .== 1] .== 1)
    fp = sum(mout[my .== 0] .== 1)
    tt = sum(my .== 1)
    ff = sum(my .== 0)

    mydiv(a,b) = if a == 0 0 else a / b end

    prec = mydiv(tp, tp+fp)
    recall = mydiv(tp, tt)

    tpr = mydiv(tp, tt)
    fpr = mydiv(fp, ff)
    fdr = mydiv(fp, sum(mout .== 1))

    shd = sum(my .!= mout) / nbatch

    GraphMetric(nnz, nny, tpr, fpr, fdr, shd, prec, recall)
end

function create_save_cb(model_dir, model)
    if !isdir(model_dir) mkpath(model_dir) end
    function save_cb(step)
        model_file = joinpath(model_dir, "step-"*string(step)*".bson")
        # @info "saving .."
        # sometimes user interrupt would crash the file, I'm saving it as
        # xxx-tmp, and perform move afterwards, hopefully fix this.
        @save model_file model=cpu(model)
        # CAUTION monitor saving time for potential performance issue
        # @time
    end
end

function create_print_cb(;logger=nothing)
    function f(step, ms)
        # println()
        # evaluate to extract from metrics. This will only call every seconds
        values = map(ms) do x
            x[1]=>get!(x[2])
        end # |> Dict
        # DEBUG removing stdout data logging
        # @info "data" values
        if typeof(logger) <: TBLogger
            for value in values
                log_value(logger, value[1], value[2], step=step)
            end
        end
    end
end

function myσxent(logŷ, y)
    # FIXME performance the gradient of this cost a tons of time, and seems to
    # be moving data around
    #
    # Related:
    # https://github.com/JuliaGPU/CuArrays.jl/issues/611
    # https://github.com/JuliaGPU/CuArrays.jl/pull/602
    # https://github.com/JuliaGPU/CuArrays.jl/issues/141
    xent = Flux.Losses.logitbinarycrossentropy.(logŷ, y)
    loss = sum(xent)
    return loss * 1 // size(y)[end]
end

"""This function is used purely for evaluation
"""
function sup_test(model, test_ds; nbatch=test_ds.nbatch)
    gm = MeanMetric{GraphMetric}()
    loss_metric = MeanMetric{Float64}()

    t1 = time()
    # FIXME testmode!
    @showprogress 0.1 "Inner testing..." for i in 1:nbatch
        x, y = next_batch!(test_ds) |> gpu
        out = model(x)
        loss = myσxent(out, y)
        metric = sup_graph_metrics(cpu(σ.(out)), cpu(y))
        add!(gm, metric)
        add!(loss_metric, loss)
    end
    t2 = time()
    t = (t2 - t1) / nbatch

    g_v = get!(gm)
    loss_v = get!(loss_metric)
    @info "data" g_v loss_v
    # return the metrics
    g_v, t
end

function sup_test_raw(model, test_ds; nbatch=test_ds.nbatch)
    error("raw test disabled!")
    gm = MeanMetric{GraphMetric}()
    loss_metric = MeanMetric{Float64}()

    t1 = time()
    # FIXME testmode!
    @showprogress 0.1 "Inner testing..." for i in 1:nbatch
        x, y = next_batch!(test_ds)

        # shit!!! this should be cor or cov or cov/norm according to the true data
        # I'm thus disabling this function
        x = cat([cor(x[:,:,i]) for i in 1:size(x, 3)]..., dims=3) |> gpu
        y = y |> gpu

        out = model(x)
        loss = myσxent(out, y)
        metric = sup_graph_metrics(cpu(σ.(out)), cpu(y))
        add!(gm, metric)
        add!(loss_metric, loss)
    end

    t2 = time()
    t = (t2 - t1) / nbatch
    # FIXME actually I need to mention the batch size in the paper
    # * test_ds.batch_size

    g_v = get!(gm)
    loss_v = get!(loss_metric)
    @info "data" g_v loss_v
    # return the metrics
    g_v, t
end

function create_test_cb(model, test_ds, msg; logger=nothing)
    function test_cb(step)
        test_run_steps = 20

        # println()
        # @info "testing for $test_run_steps steps .."
        gm = MeanMetric{GraphMetric}()
        loss_metric = MeanMetric{Float64}()

        # FIXME testmode!
        # @showprogress 0.1 "Inner testing..."
        for i in 1:test_run_steps
            x, y = next_batch!(test_ds) |> gpu
            out = model(x)

            loss = myσxent(out, y)

            # FIXME performance on CPU
            metric = sup_graph_metrics(cpu(σ.(out)), cpu(y))

            add!(gm, metric)
            add!(loss_metric, loss)
        end

        g_v = get!(gm)
        loss_v = get!(loss_metric)
        # @info msg loss_v to_named_tuple(g_v)...
        if typeof(logger) <: TBLogger
            # log_value(logger, "graphM/" * msg, value, step=step)
            # DEBUG using tuple
            log_value(logger, "loss", loss_v, step=step)
            log_value(logger, "graph", g_v, step=step)
            # log_value(logger, "graph", to_named_tuple(g_v), step=step)
        end
    end
    # execute every 10 seconds
    # Flux.throttle(test_cb, 10)
end

function sup_train!(model, opt, ds;
                    from_step=1,
                    train_steps=ds.nbatch,
                    print_cb=(i, m)->(),
                    save_cb=(i)->(),
                    test_cb=(i)->())
    train_steps = convert(Int, train_steps)

    ps=Flux.params(model)
    # weight decay or all params decay?
    weights = weight_params(model)

    loss_metric = MeanMetric{Float64}()
    gm = MeanMetric{GraphMetric}()

    @info "training for $(train_steps) steps .."
    @showprogress 0.1 "Training..." for step in from_step:train_steps
        # CAUTION this actually cost half of the computing time.
        x, y = next_batch!(ds) |> gpu

        # FIXME evaluating this too frequently might be inefficient
        #
        # FIXME and it turns out I cannot do this here, otherwise Zygote
        # complains. But now I have to run one more forward pass to get the
        # data.
        #
        # FIXME performance hell!
        # metric = sup_graph_metrics(cpu(σ.(model(x))), cpu(y))
        # add!(gm, metric)

        gs = gradient(ps) do
            out = model(x)
            loss = myσxent(out, y)

            # add a weight decay
            # l2 = sum((x)->sum(x.^2), weights)
            # show l2?
            # loss = loss + 1e-5 * l2
            add!(loss_metric, loss)

            loss
        end

        Flux.Optimise.update!(opt, ps, gs)

        # "graph"=>gm
        print_cb(step, ["loss"=>loss_metric])

        test_cb(step)
        save_cb(step)
    end
end

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
using Query
using DataFrames
using NamedTuples

df = DataFrame(a=[2,1,1,2,1,3],b=[2,2,1,1,3,2])

x = @from i in df begin
    @orderby descending(i.a), i.b
    @select i
    @collect DataFrame
end

println(x)


text:
#Function to show example fits for a given pixel in a SAIM image sequence
function plot_fit(image_file, JLD_file, row, col)

	println(image_file)
	println(JLD_file)

	#Load the constants from the JLD_file
	constants = load(JLD_file, "constants")
	angles = load(JLD_file, "angles")
	angles = RTD*angles		#convert from radians to degrees

	#Load the image sequence and extract the channel intensities at pixel (row, col)
	img_HWC = load(image_file);
    channels = channelview(float64.(img_HWC))*65535; 	#Each channel corresponds to a specific excitation angle
    ydata = channels[row,col,:]

    #Load the fit results
    A = load(JLD_file, "A")
    B = load(JLD_file, "B")
	H = load(JLD_file, "H")
    constants = load(JLD_file, "constants")
    p = [A[row,col], B[row,col], H[row,col]]

	#Plot the example fit
	yfit = model(constants, p)
    colors = palette(:tab10)
	plt = plot(angles, ydata, label="data", seriestype = :scatter, markercolor = colors[1])
	plot!(angles, yfit, label="fit", linecolor = colors[4], linewidth = 1.5)
    xlabel!("Incidence angle (degrees)")
	ylabel!("Intensity (A.U.)")

    savefig(plt, "example_fit.png")
	savefig(plt, "example_fit.pdf")
	savefig(plt, "example_fit.svg")
    display(plt)
end

#Plot a SAIM curve
function plot_curve()

	#-------------- OPTICAL PARAMETERS----------------------
	nB = 1.33;			#The refractive index of the ambient media / cytoplasm
	dOx = 1924.9;		#The thickness of the SiO2 layer in units of nm
	λ_Ex_1 = 642.0; 		#The wavelength of excitation laser #1 in units of nm
	nOx_1 = 1.463;		#The refractive index of SiO2 at excitation wavelength #1
	nSi_1 = 4.3638; 	#The complex refractive index of Si at excitation wavelength #1; #4.3 + 0.073im;

	opt1 = Optical(nB, nOx_1, nSi_1, dOx, λ_Ex_1)

	#-------------- ANGLES  ---------------------
	angles_deg = range(5.0, length=100, stop=43.75)	#Incendece angle in degree for each image frame
	angles_rad = angles_deg*DTR						#Incendece angle in radians for each image frame
	num_angles = length(angles_rad)

	#-------------- SIMULATION----------------------
	colors = palette(:tab10)

	constants = Calculate_Constants1c(opt1, angles_rad)

	Pex = Pex_from_Constants(94.0, constants)
	plt = plot(angles_deg, Pex, linecolor = colors[1],
		linewidth = 2, label = "H = 94 nm", legend=:outertopright)

	Pex = Pex_from_Constants(365.0, constants)
	plot!(angles_deg, Pex, linecolor = colors[4],
		linewidth = 2, label = "H = 365 nm")

	xlabel!("Incidence angle (degrees)")
	ylabel!("Intensity (A.U.)")
	display(plt)
	savefig(plt, "94_365_642nm.png")
	savefig(plt, "94_365_642nm.pdf")
	savefig(plt, "94_365_642nm.svg")

end


text:
#
# misc.jl --
#
# Implement XPA configuration methods and miscellaneous methods.
#
#------------------------------------------------------------------------------
#
# This file is part of XPA.jl released under the MIT "expat" license.
# Copyright (C) 2016-2019, Éric Thiébaut (https://github.com/emmt/XPA.jl).
#

#------------------------------------------------------------------------------
# CONFIGURATION METHODS

# The following default values are defined in "xpap.c" and can be changed by
# user environment variables.
const _DEFAULTS = Dict{String,Any}("XPA_MAXHOSTS" => 100,
                                   "XPA_SHORT_TIMEOUT" => 15,
                                   "XPA_LONG_TIMEOUT" => 180,
                                   "XPA_CONNECT_TIMEOUT" => 10,
                                   "XPA_TMPDIR" => "/tmp/.xpa",
                                   "XPA_VERBOSITY" => true,
                                   "XPA_IOCALLSXPA" => false)

"""
```julia
XPA.getconfig(key) -> val
```

yields the value associated with configuration parameter `key` (a string or a
symbol).  The following parameters are available (see XPA doc. for more
information):

| Key Name                | Default Value |
|:----------------------- |:------------- |
| `"XPA_MAXHOSTS"`        | `100`         |
| `"XPA_SHORT_TIMEOUT"`   | `15`          |
| `"XPA_LONG_TIMEOUT"`    | `180`         |
| `"XPA_CONNECT_TIMEOUT"` | `10`          |
| `"XPA_TMPDIR"`          | `"/tmp/.xpa"` |
| `"XPA_VERBOSITY"`       | `true`        |
| `"XPA_IOCALLSXPA"`      | `false`       |

Also see [`XPA.setconfig!`](@ref).

"""
function getconfig(key::AbstractString)
    haskey(_DEFAULTS, key) || error("unknown XPA parameter \"$key\"")
    def = _DEFAULTS[key]
    if haskey(ENV, key)
        val = haskey(ENV, key)
        return (isa(def, Bool) ? (parse(Int, val) != 0) :
                isa(def, Integer) ? parse(Int, val) :
                isa(def, AbstractString) ? val :
                error("unexpected type $(typeof(def)) for default value of \"$key\""))
    else
        return def
    end
end

"""
```julia
XPA.setconfig!(key, val) -> oldval
```

set the value associated with configuration parameter `key` to be `val`.  The
previous value is returned.

Also see [`XPA.getconfig`](@ref).

"""
function setconfig!(key::AbstractString,
                    val::T) where {T<:Union{Integer,Bool,AbstractString}}
    global _DEFAULTS, ENV
    old = getconfig(key) # also check validity of key
    def = _DEFAULTS[key]
    if isa(def, Integer) && isa(val, Integer)
        ENV[key] = string(val)
    elseif isa(def, Bool) && isa(val, Bool)
        ENV[key] = (val ? "1" : "0")
    elseif isa(def, AbstractString) && isa(val, AbstractString)
        ENV[key] = val
    else
        error("invalid type for XPA parameter \"$key\"")
    end
    return old
end

getconfig(key::Symbol) = getconfig(string(key))
setconfig!(key::Symbol, val) = setconfig!(string(key), val)

@deprecate config(key::AbstractString) getconfig(key)
@deprecate config!(key::AbstractString, val) setconfig!(key, val)

#------------------------------------------------------------------------------
# PRIVATE METHODS

"""
Private methods:

```julia
_get_field(T, ptr, off, def)
```

and

```julia
_get_field(T, ptr, off1, off2, def)
```

retrieve a field of type `T` at offset `off` (in bytes) with respect to address
`ptr`.  If two offsets are given, the first one refers to a pointer with
respect to which the second is applied.  If `ptr` is NULL, `def` is returned.

"""
_get_field(::Type{T}, ptr::Ptr{Cvoid}, off::Int, def::T) where {T} =
    (ptr == C_NULL ? def : unsafe_load(convert(Ptr{T}, ptr + off)))

_get_field(::Type{String}, ptr::Ptr{Cvoid}, off::Int, def::String) =
    (ptr == C_NULL ? def : unsafe_string(convert(Ptr{Ptr{Byte}}, ptr + off)))

function _get_field(::Type{T}, ptr::Ptr{Cvoid}, off1::Int, off2::Int,
                    def::T) where {T}
    _get_field(T, _get_field(Ptr{Cvoid}, ptr, off1, C_NULL), off2, def)
end

function _set_field(::Type{T}, ptr::Ptr{Cvoid}, off::Int, val) where {T}
    @assert ptr != C_NULL
    unsafe_store!(convert(Ptr{T}, ptr + off), val)
end

_get_comm(xpa::Handle) =
    _get_field(Ptr{Cvoid}, xpa.ptr, _offsetof_comm, C_NULL)

for (memb, T, def) in ((:name,      String, ""),
                       (:class,     String, ""),
                       (:send_mode, Cint,   Cint(0)),
                       (:recv_mode, Cint,   Cint(0)),
                       (:method,    String, ""),
                       (:sendian,   String, "?"))
    off = Symbol(:_offsetof_, memb)
    func = Symbol(:get_, memb)
    @eval begin
        $func(xpa::Handle) = _get_field($T, xpa.ptr, $off, $def)
    end
end

for (memb, T, def) in ((:comm_status,  Cint,       Cint(0)),
                       (:comm_cmdfd,   Cint,       Cint(-1)),
                       (:comm_datafd,  Cint,       Cint(-1)),
                       (:comm_ack,     Cint,       Cint(1)),
                       (:comm_cendian, String,     "?"),
                       (:comm_buf,     Ptr{Byte},  NULL),
                       (:comm_len,     Csize_t,    Csize_t(0)))
    off = Symbol(:_offsetof_, memb)
    func = Symbol(:get_, memb)
    @eval begin
        $func(xpa::Handle) = _get_field($T, _get_comm(xpa), $off, $def)
    end
    if memb == :comm_buf || memb == :comm_len
        func = Symbol(:_set_, memb)
        @eval begin
            $func(xpa::Handle, val) =
                unsafe_store!(convert(Ptr{$T}, _get_comm(xpa) + $off), val)
        end
    end
end

"""
```julia
_malloc(len)
```

dynamically allocates `len` bytes and returns the corresponding byte pointer
(type `Ptr{UInt8}`).

"""
function _malloc(len::Integer) :: Ptr{Byte}
    ptr = ccall(:malloc, Ptr{Byte}, (Csize_t,), len)
    ptr != NULL || throw(OutOfMemoryError())
    return ptr
end

"""
```julia
_free(ptr)
```

frees dynamically allocated memory at address givne by `ptr` unless it is NULL.

"""
_free(ptr::Ptr{T}) where T =
    (ptr == Ptr{T}(0) || ccall(:free, Cvoid, (Ptr{T},), ptr))

"""
```julia
_memcpy!(dst, src, len) -> dst
```

copies `len` bytes from address `src` to destination `dst`.

"""
function _memcpy!(dst::Ptr, src::Ptr, len::Integer)
    if len > 0
        ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
              dst, src, len)
    end
    return dst
end

function _memcpy!(dst::AbstractArray, src::Ptr, len::Integer)
    len == sizeof(dst) || error("bad number of bytes to copy")
    if len > 0
        ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t),
              dst, src, len)
    end
    return dst
end


text:
module Misc

using LinearAlgebra
export dropdim1, dropdim2, eye

dropdim1(x) = dropdims(x, dims=1)   # useful for pipes
dropdim2(x) = dropdims(x, dims=2)

eye(d) = Matrix(I, d, d)

end


text:
struct FontHeader
    flags::UInt16
    units_per_em::UInt16
    created::DateTime
    modified::DateTime
    xmin::Int16
    ymin::Int16
    xmax::Int16
    ymax::Int16
    mac_style::UInt16
    min_readable_size::Int16
    font_direction_hint::Int16
    index_to_loc_format::Int16
    glyph_data_format::Int16
end

function Base.parse(io::IO, ::Type{FontHeader})
    skip(io, 12)
    read(io, UInt32) == 0x5f0f3cf5 || error("Invalid magic number in font header")
    FontHeader(
        read(io, UInt16),
        read(io, UInt16),
        DateTime(1904, 01, 01) + Dates.Second(read(io, Int64)),
        DateTime(1904, 01, 01) + Dates.Second(read(io, Int64)),
        (read(io, T) for T in fieldtypes(FontHeader)[5:end])...
    )
end


text:
#
# Device functionality
#

## execution

struct ROCArrayBackend <: AbstractGPUBackend end

struct ROCKernelContext <: AbstractKernelContext end

function GPUArrays.gpu_call(::ROCArrayBackend, f, args, threads::Int, blocks::Int;
                            name::Union{String,Nothing})
    groupsize, gridsize = threads, blocks*threads
    wait(@roc groupsize=groupsize gridsize=gridsize f(ROCKernelContext(), args...))
end

## on-device

# indexing

for (f, froc) in (
        (:blockidx, :blockIdx),
        (:blockdim, :blockDim),
        (:threadidx, :threadIdx),
        (:griddim, :gridDimWG)
    )
    @eval GPUArrays.$f(::ROCKernelContext) = AMDGPU.$froc().x
end

# math

@inline GPUArrays.cos(ctx::ROCKernelContext, x) = cos(x)
@inline GPUArrays.sin(ctx::ROCKernelContext, x) = sin(x)
@inline GPUArrays.sqrt(ctx::ROCKernelContext, x) = sqrt(x)
@inline GPUArrays.log(ctx::ROCKernelContext, x) = log(x)

# memory

@inline function GPUArrays.LocalMemory(ctx::ROCKernelContext, ::Type{T}, ::Val{dims}, ::Val{id}) where {T,dims,id}
    ptr = AMDGPU.alloc_special(Val(id), T, AMDGPU.AS.Local, Val(prod(dims)))
    ROCDeviceArray(dims, ptr)
end

# synchronization

@inline function GPUArrays.synchronize_threads(::ROCKernelContext)
    AMDGPU.sync_workgroup()
    return
end

#
# Host abstractions
#

mutable struct ROCArray{T,N} <: AbstractGPUArray{T,N}
    buf::Mem.Buffer
    own::Bool

    dims::Dims{N}
    offset::Int

    function ROCArray{T,N}(buf::Mem.Buffer, dims::Dims{N}; offset::Integer=0, own::Bool=true) where {T,N}
        @assert isbitstype(T) "ROCArray only supports bits types"
        xs = new{T,N}(buf, own, dims, offset)
        if own
            hsaref!()
            Mem.retain(buf)
            finalizer(unsafe_free!, xs)
        end
        return xs
    end
end

function unsafe_free!(xs::ROCArray)
    Mem.release(xs.buf) && Mem.free(xs.buf)
    hsaunref!()
    return
end

## aliases

const ROCVector{T} = ROCArray{T,1}
const ROCMatrix{T} = ROCArray{T,2}
const ROCVecOrMat{T} = Union{ROCVector{T},ROCMatrix{T}}

## constructors

# type and dimensionality specified, accepting dims as tuples of Ints
function ROCArray{T,N}(::UndefInitializer, dims::Dims{N}) where {T,N}
    buf = Mem.alloc(prod(dims)*sizeof(T))
    ROCArray{T,N}(buf, dims)
end

# type and dimensionality specified, accepting dims as series of Ints
ROCArray{T,N}(::UndefInitializer, dims::Integer...) where {T,N} = ROCArray{T,N}(undef, dims)

# type but not dimensionality specified
ROCArray{T}(::UndefInitializer, dims::Dims{N}) where {T,N} = ROCArray{T,N}(undef, dims)
ROCArray{T}(::UndefInitializer, dims::Integer...) where {T} =
    ROCArray{T}(undef, convert(Tuple{Vararg{Int}}, dims))

# from Base arrays
function ROCArray{T,N}(x::Array{T,N}, dims::Dims{N}) where {T,N}
    r = ROCArray{T,N}(undef, size(x))
    Mem.upload!(r.buf, pointer(x), sizeof(x))
    return r
end

# type as first argument
# FIXME: Remove me!
#ROCArray(::Type{T}, dims::Dims{N}) where {T,N} = ROCArray{T,N}(undef, dims)

# empty vector constructor
ROCArray{T,1}() where {T} = ROCArray{T,1}(undef, 0)

Base.similar(a::ROCArray{T,N}) where {T,N} = ROCArray{T,N}(undef, size(a))
Base.similar(a::ROCArray{T}, dims::Base.Dims{N}) where {T,N} = ROCArray{T,N}(undef, dims)
Base.similar(a::ROCArray, ::Type{T}, dims::Base.Dims{N}) where {T,N} = ROCArray{T,N}(undef, dims)

## array interface

Base.elsize(::Type{<:ROCArray{T}}) where {T} = sizeof(T)

Base.size(x::ROCArray) = x.dims
Base.sizeof(x::ROCArray) = Base.elsize(x) * length(x)

## interop with Julia arrays

ROCArray{T,N}(x::AbstractArray{S,N}) where {T,N,S} =
    ROCArray{T,N}(convert(Array{T}, x), size(x))

# underspecified constructors
ROCArray{T}(xs::AbstractArray{S,N}) where {T,N,S} = ROCArray{T,N}(xs)
(::Type{ROCArray{T,N} where T})(x::AbstractArray{S,N}) where {S,N} = ROCArray{S,N}(x)
ROCArray(A::AbstractArray{T,N}) where {T,N} = ROCArray{T,N}(A)

# idempotency
ROCArray{T,N}(xs::ROCArray{T,N}) where {T,N} = xs

## conversions

Base.convert(::Type{T}, x::T) where T <: ROCArray = x

## memory operations

function Base.copyto!(dest::Array{T}, d_offset::Integer,
                      source::ROCArray{T}, s_offset::Integer,
                      amount::Integer) where T
    amount == 0 && return dest
    @boundscheck checkbounds(dest, d_offset+amount-1)
    @boundscheck checkbounds(source, s_offset+amount-1)
    Mem.download!(pointer(dest, d_offset),
                  Mem.view(source.buf, (s_offset-1)*sizeof(T)),
                  amount*sizeof(T))
    dest
end
function Base.copyto!(dest::ROCArray{T}, d_offset::Integer,
                      source::Array{T}, s_offset::Integer,
                      amount::Integer) where T
    amount == 0 && return dest
    @boundscheck checkbounds(dest, d_offset+amount-1)
    @boundscheck checkbounds(source, s_offset+amount-1)
    Mem.upload!(Mem.view(dest.buf, (d_offset-1)*sizeof(T)),
                pointer(source, s_offset),
                amount*sizeof(T))
    dest
end
function Base.copyto!(dest::ROCArray{T}, d_offset::Integer,
                      source::ROCArray{T}, s_offset::Integer,
                      amount::Integer) where T
    amount == 0 && return dest
    @boundscheck checkbounds(dest, d_offset+amount-1)
    @boundscheck checkbounds(source, s_offset+amount-1)
    Mem.transfer!(Mem.view(dest.buf, (d_offset-1)*sizeof(T)),
                  Mem.view(source.buf, (s_offset-1)*sizeof(T)),
                  amount*sizeof(T))
    dest
end

# TODO: Workaround for hanging copy() broadcast kernel
function Base.copy(X::ROCArray{T}) where T
    Xnew = ROCArray{T}(undef, size(X))
    copyto!(Xnew, 1, X, 1, length(X))
    Xnew
end

## views

# optimize view to return a ROCArray when contiguous

struct Contiguous end
struct NonContiguous end

# NOTE: this covers more cases than the I<:... in Base.FastContiguousSubArray
ROCIndexStyle() = Contiguous()
ROCIndexStyle(I...) = NonContiguous()
ROCIndexStyle(::Base.ScalarIndex...) = Contiguous()
ROCIndexStyle(i1::Colon, ::Base.ScalarIndex...) = Contiguous()
ROCIndexStyle(i1::AbstractUnitRange, ::Base.ScalarIndex...) = Contiguous()
ROCIndexStyle(i1::Colon, I...) = ROCIndexStyle(I...)

rocviewlength() = ()
@inline rocviewlength(::Real, I...) = rocviewlength(I...) # skip scalar
@inline rocviewlength(i1::AbstractUnitRange, I...) = (Base.unsafe_length(i1), rocviewlength(I...)...)
@inline rocviewlength(i1::AbstractUnitRange, ::Base.ScalarIndex...) = (Base.unsafe_length(i1),)

@inline function Base.view(A::ROCArray, I::Vararg{Any,N}) where {N}
    J = to_indices(A, I)
    @boundscheck begin
        # Base's boundscheck accesses the indices, so make sure they reside on the CPU.
        # this is expensive, but it's a bounds check after all.
        J_cpu = map(j->adapt(Array, j), J)
        checkbounds(A, J_cpu...)
    end
    J_gpu = map(j->adapt(ROCArray, j), J)
    unsafe_view(A, J_gpu, ROCIndexStyle(I...))
end

@inline function unsafe_view(A, I, ::Contiguous)
    unsafe_contiguous_view(Base._maybe_reshape_parent(A, Base.index_ndims(I...)), I, rocviewlength(I...))
end
@inline function unsafe_contiguous_view(a::ROCArray{T}, I::NTuple{N,Base.ViewIndex}, dims::NTuple{M,Integer}) where {T,N,M}
    offset = Base.compute_offset1(a, 1, I) * sizeof(T)

    Mem.retain(a.buf)
    b = ROCArray{T,M}(a.buf, dims, offset=a.offset+offset, own=false)
    finalizer(unsafe_free!, b)
    return b
end

@inline function unsafe_view(A, I, ::NonContiguous)
    Base.unsafe_view(Base._maybe_reshape_parent(A, Base.index_ndims(I...)), I...)
end

#FIXME: add pointer conversions

## reshape

# optimize reshape to return a ROCArray

function Base.reshape(a::ROCArray{T,M}, dims::NTuple{N,Int}) where {T,N,M}
    if prod(dims) != length(a)
        throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $len"))
    end

    if N == M && dims == size(a)
        return a
    end

    Mem.retain(a.buf)
    b = ROCArray{T,N}(a.buf, dims, offset=a.offset, own=false)
    finalizer(unsafe_free!, b)
    return b
end

# allow missing dimensions with Colon()
if VERSION < v"1.6.0-DEV.1358"
    Base.reshape(parent::ROCArray, dims::Tuple{Vararg{Union{Int,Colon}}}) =
    Base.reshape(parent, Base._reshape_uncolon(parent, dims))
end

## fft

#=
using AbstractFFTs

# defining our own plan type is the easiest way to pass around the plans in FFTW interface
# without ambiguities

struct FFTPlan{T}
    p::T
end

AbstractFFTs.plan_fft(A::ROCArray; kw_args...) = FFTPlan(plan_fft(A.data; kw_args...))
AbstractFFTs.plan_fft!(A::ROCArray; kw_args...) = FFTPlan(plan_fft!(A.data; kw_args...))
AbstractFFTs.plan_bfft!(A::ROCArray; kw_args...) = FFTPlan(plan_bfft!(A.data; kw_args...))
AbstractFFTs.plan_bfft(A::ROCArray; kw_args...) = FFTPlan(plan_bfft(A.data; kw_args...))
AbstractFFTs.plan_ifft!(A::ROCArray; kw_args...) = FFTPlan(plan_ifft!(A.data; kw_args...))
AbstractFFTs.plan_ifft(A::ROCArray; kw_args...) = FFTPlan(plan_ifft(A.data; kw_args...))

function Base.:(*)(plan::FFTPlan, A::ROCArray)
    x = plan.p * A.data
    ROCArray(x)
end
=#

## GPUArrays interfaces

GPUArrays.device(x::ROCArray) = x.buf.agent

GPUArrays.backend(::Type{<:ROCArray}) = ROCArrayBackend()

function Base.convert(::Type{ROCDeviceArray{T,N,AS.Global}}, a::ROCArray{T,N}) where {T,N}
    ptr = Base.unsafe_convert(Ptr{T}, a.buf)
    ROCDeviceArray{T,N,AS.Global}(a.dims, AMDGPU.LLVMPtr{T,AS.Global}(ptr+a.offset))
end
Adapt.adapt_storage(::AMDGPU.Adaptor, x::ROCArray{T,N}) where {T,N} =
    convert(ROCDeviceArray{T,N,AS.Global}, x)

## interop with CPU arrays

# We don't convert isbits types in `adapt`, since they are already
# considered GPU-compatible.

Adapt.adapt_storage(::Type{ROCArray}, xs::AT) where {AT<:AbstractArray} =
  isbitstype(AT) ? xs : convert(ROCArray, xs)

# if an element type is specified, convert to it
Adapt.adapt_storage(::Type{<:ROCArray{T}}, xs::AT) where {T, AT<:AbstractArray} =
  isbitstype(AT) ? xs : convert(ROCArray{T}, xs)

Adapt.adapt_storage(::Type{Array}, xs::ROCArray) = convert(Array, xs)

## Float32-preferring conversion

struct Float32Adaptor end

Adapt.adapt_storage(::Float32Adaptor, xs::AbstractArray) =
  isbits(xs) ? xs : convert(ROCArray, xs)

Adapt.adapt_storage(::Float32Adaptor, xs::AbstractArray{<:AbstractFloat}) =
  isbits(xs) ? xs : convert(ROCArray{Float32}, xs)

Adapt.adapt_storage(::Float32Adaptor, xs::AbstractArray{<:Complex{<:AbstractFloat}}) =
  isbits(xs) ? xs : convert(ROCArray{ComplexF32}, xs)

# not for Float16
Adapt.adapt_storage(::Float32Adaptor, xs::AbstractArray{Float16}) =
  isbits(xs) ? xs : convert(ROCArray, xs)

roc(xs) = adapt(Float32Adaptor(), xs)

Base.unsafe_convert(::Type{Ptr{T}}, x::ROCArray{T}) where T =
    Base.unsafe_convert(Ptr{T}, x.buf)


text:
module TestGraphics

import Modia3D
using  Modia3D.Test
using  Modia3D.StaticArrays

material1 = Modia3D.Material(;color="DarkRed", wireframe=true)
material2 = Modia3D.Material(;color=[0,0,255], shininess=0.2)

sphere    = Modia3D.Sphere(2.0, material=material1)
ellipsoid = Modia3D.Ellipsoid(1.0,2.0,3.0, material=material1)
box       = Modia3D.Box(1.0,2.0,3.0, material=material1)
cylinder  = Modia3D.Cylinder(3.0,1.0;Dy=0., material=material1)
capsule   = Modia3D.Capsule(3.0,1.0;Dy=0.5, material=material1)
beam      = Modia3D.Beam(1.0,2.0,3.0, material=material1)
cone      = Modia3D.Cone(3.0,1.0;Dy=0.5, relativeTipDiameter=0.2,
                         relativeInnerDiameter=0.5, material=material1)
pipe      = Modia3D.Pipe(3.0,1.0;Dy=0.5, relativeTipDiameter=0.2,
                         relativeInnerDiameter=0.5, material=material1)
spring    = Modia3D.Spring(3.0,1.0;Dy=0.5, windings=10,springRadius=0.1, material=material2)
gearWheel = Modia3D.GearWheel(3.0,1.0;Dy=0.5,relativeInnerDiameter=0.8,teeth=30,
                              angle=0.2, material=material2)
coordsys  = Modia3D.CoordinateSystem(2.0)
grid      = Modia3D.Grid(3.0,1.0, distance=0.1, lineWidth=2.0)

filename = joinpath(Modia3D.path, "objects", "engine", "crank", "crank.obj")
mesh1 = Modia3D.FileMesh(filename)
mesh2 = Modia3D.FileMesh(filename, 2.0; useMaterialColor=true, smoothNormals=true)
mesh3 = Modia3D.FileMesh(filename; scaleFactor=[2.0,3.0,4.0], useMaterialColor=true, smoothNormals=true)

text1 = Modia3D.TextShape("This is a text line")
font1 = Modia3D.Font()
font2 = Modia3D.Font(fontFamily="FreeSans", bold=false, italic=false, charSize=0.1,
                     color = "LightBlue", transparency=0.0)
text2 = Modia3D.TextShape("Another text line", font=font2, offset=[0.1,0.2,0.3],
                          axisAlignment=Modia3D.XZ_Plane,alignment=Modia3D.Left)

@testset "Modia3D.Graphics: Test geometry, color, fonts" begin
   @test material1.color     == MVector{3,Cint}(155,0,0)
   @test material1.wireframe == true
end

end


text:
#### DEFINE THE EVALUATION PARAMETERS ################

"""
Container type for the evaluation script parameters
#Fields
- `sim_dt` simulation time step
- rng::AbstractRNG
- n_episodes::Int64
- time_out::Int64
- callbacks::Tuple{Vararg{Any}}
"""
struct EvalConfig
    sim_dt::Float64
    rng::AbstractRNG
    n_episodes::Int64
    time_out::Int64
    callbacks::Tuple{Vararg{Any}}
    n_ped::Int64
    #verbose
    #render
end

function EvalConfig(; sim_dt::Float64 = 0.1,
                      rng::AbstractRNG = MersenneTwister(3),
                      n_episodes::Int64 = 1,
                      time_out::Int64 = 200,
                      callbacks::Tuple{Vararg{Any}} = (EgoCollisionCallback(),TerminalCallback()),
                      n_ped::Int64 = 2
                      )
    return EvalConfig(sim_dt, rng, n_episodes, time_out, callbacks, n_ped)
end

function Base.show(io::IO, conf::EvalConfig)
    @printf(io, "EvalConfig \n")
    @printf(io, "\t Sim step (s) %.2f \n", conf.sim_dt)
    @printf(io, "\t N episodes %4i \n", conf.n_episodes)
    @printf(io, "\t Time out %4i \n", conf.time_out)
    @printf(io, "\t Time out %4i \n", conf.time_out)
end

## Define Callbacks

struct TerminalCallback end

struct EgoCollisionCallback end

function AutomotiveDrivingModels.run_callback{S,D,I,M<:DriverModel}(
    callback::EgoCollisionCallback,
    rec::EntityQueueRecord{S,D,I},
    env::CrosswalkEnv,
    models::Dict{I,M},
    tick::Int,
    )
    return is_crash(rec[0])
end

"""
    is_crash(scene::Scene)
return true if the ego car is in collision in the given scene, do not check for collisions between
other participants
"""
function is_crash(scene::Scene)
    ego = scene[findfirst(EGO_ID, scene)]
    @assert ego.id == 1
    if ego.state.v ≈ 0
        return false
    end
    for veh in scene
        if veh.id != 1
            if is_colliding(ego, veh)
                return true
            end
        end
    end
    return false
end

function AutomotiveDrivingModels.run_callback{S,D,I,M<:DriverModel}(
    callback::TerminalCallback,
    rec::EntityQueueRecord{S,D,I},
    env::CrosswalkEnv,
    models::Dict{I,M},
    tick::Int,
    )

    return is_terminal(rec[0], env)
end

function is_terminal(scene::Scene, env::CrosswalkEnv)
    # check ego car position
    ego = scene[findfirst(EGO_ID, scene)]
    @assert ego.id == 1
    return ego.state.posG.x > env.params.end_pos
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs