text:
@testset "THREADPOOL" begin
    @test isconst(ThreadingUtilities, :THREADPOOL) # test that ThreadingUtilities.THREADPOOL is a constant
    @test ThreadingUtilities.THREADPOOL isa Vector{UInt}
    @test eltype(ThreadingUtilities.THREADPOOL) === UInt
    @test length(ThreadingUtilities.THREADPOOL) == (ThreadingUtilities.THREADBUFFERSIZE÷sizeof(UInt)) * (min(Threads.nthreads(),(Sys.CPU_THREADS)::Int) - 1) + (256 ÷ sizeof(UInt)) - 1
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).
@recipe function f(
    maybe_shaped_samples::DensitySampleVector,
    vsel::NTuple{2,Union{Symbol, Expr, Integer}};
    intervals = default_credibilities,
    interval_labels = [],
    colors = default_colors,
    mean = false,
    std = false,
    globalmode = false,
    marginalmode = true,
    smoothing = 0,
    diagonal = Dict(),
    upper = Dict(),
    right = Dict(),
    filter = false,
    closed = :left
)
    xindx = asindex(maybe_shaped_samples, vsel[1])
    yindx = asindex(maybe_shaped_samples, vsel[2])

    if length(xindx) > 1
        throw(ArgumentError("Symbol :$(vsel[1]) refers to a multivariate parameter. Use :($(vsel[1])[i]) instead."))
    elseif length(yindx) > 1
        throw(ArgumentError("Symbol :$(vsel[2]) refers to a multivariate parameter. Use :($(vsel[2])[i]) instead."))
    end

    samples = unshaped.(maybe_shaped_samples)
    filter ? samples = BAT.drop_low_weight_samples(samples) : nothing

    bins = get(plotattributes, :bins, 200)
    seriestype = get(plotattributes, :seriestype, :smallest_intervals)

    colorbar = false
    if seriestype == :histogram2d || seriestype == :histogram || seriestype == :hist
        colorbar = get(plotattributes, :colorbar, false)
    end

    xlabel, ylabel  = if !isshaped(maybe_shaped_samples)
        "v$xindx", "v$yindx"
    else
        getstring(maybe_shaped_samples, xindx), getstring(maybe_shaped_samples, yindx)
    end

    xguide := get(plotattributes, :xguide, xlabel)
    yguide := get(plotattributes, :yguide, ylabel)

    marg = MarginalDist(
        samples,
        (xindx, yindx),
        bins = bins,
        closed = closed,
        filter = filter
    )

    if seriestype == :scatter
        base_markersize = get(plotattributes, :markersize, 1.5)

        acc = findall(x -> x > 0, samples.weight)
        rej = findall(x -> x <= 0, samples.weight)

        color = parse(RGBA{Float64}, get(plotattributes, :seriescolor, :green))
        label = get(plotattributes, :label, isempty(rej) ? "samples" : "accepted")

        @series begin
            seriestype := :scatter
            label := label
            markersize := [w < 1 ? base_markersize : base_markersize * sqrt(w) for w in samples.weight[acc]]
            markerstrokewidth := 0
            markeralpha := 0.6
            markercolor := [w >= 1 ? color : RGBA(convert(RGB, color), color.alpha * w) for w in samples.weight[acc]]
            (flatview(samples.v)[xindx, acc], flatview(samples.v)[yindx, acc])
        end

        if !isempty(rej)
            @series begin
                seriestype := :scatter
                label := "rejected"
                markersize := base_markersize
                markerstrokewidth := 0
                markercolor := :red
                (flatview(samples.v)[xindx, rej], flatview(samples.v)[yindx, rej])
            end
        end

    else
        @series begin
            seriestype --> seriestype
            intervals --> intervals
            interval_labels --> interval_labels
            colors --> colors
            diagonal --> diagonal
            upper --> upper
            right --> right
            smoothing --> smoothing

            marg
        end
    end

    #------ stats ----------------------------
    stats = MCMCBasicStats(samples)

    mean_options = convert_to_options(mean)
    globalmode_options = convert_to_options(globalmode)
    marginalmode_options = convert_to_options(marginalmode)
    std_options = convert_to_options(std)

    if mean_options != ()
        mx= stats.param_stats.mean[xindx]
        my = stats.param_stats.mean[yindx]

        Σ_all = stats.param_stats.cov
        Σx = Σ_all[xindx, xindx]
        Σy = Σ_all[yindx, yindx]

        @series begin
            seriestype := :scatter
            label := get(mean_options, "label", "mean")
            seriestype==:marginal ? subplot := 3 :
            markeralpha := get(mean_options, "markeralpha", 1)
            markercolor := get(mean_options, "markercolor", :black)
            markersize := get(mean_options, "markersize", 4)
            markershape := get(mean_options, "markershape", :circle)
            markerstrokealpha := get(mean_options, "markerstrokealpha", 1)
            markerstrokecolor := get(mean_options, "markerstrokecolor", :black)
            markerstrokestyle := get(mean_options, "markerstrokestyle", :solid)
            markerstrokewidth := get(mean_options, "markerstrokewidth", 1)
            colorbar := colorbar

            if(std)
                xerror := sqrt(Σx)
                yerror := sqrt(Σy)
            end
           ([mx], [my])
        end
    end

   if globalmode_options != ()
        globalmode_x = stats.mode[xindx]
        globalmode_y = stats.mode[yindx]

        @series begin
            seriestype := :scatter
            label := get(globalmode_options, "label", "global mode")
            seriestype==:marginal ? subplot := 3 :
            markeralpha := get(globalmode_options, "markeralpha", 1)
            markercolor := get(globalmode_options, "markercolor", :black)
            markersize := get(globalmode_options, "markersize", 4)
            markershape := get(globalmode_options, "markershape", :rect)
            linealpha := 0
            markerstrokealpha := get(globalmode_options, "markerstrokealpha", 1)
            markerstrokecolor := get(globalmode_options, "markerstrokecolor", :black)
            markerstrokestyle := get(globalmode_options, "markerstrokestyle", :solid)
            markerstrokewidth := get(globalmode_options, "markerstrokewidth", 1)
            colorbar := colorbar

           ([globalmode_x], [globalmode_y])
        end
    end

    if marginalmode_options != ()
        marginalmode_values = find_marginalmodes(marg)
        for (i, l) in enumerate(marginalmode_values)
         @series begin
            seriestype := :scatter
            if i==1 && length(marginalmode_values)==1
                label := get(marginalmode_options, "label", "marginal mode")
            elseif i ==1
                label := get(marginalmode_options, "label", "marginal modes")
            else
                label :=""
            end

            seriestype == :marginal ? subplot := 3 :
            markercolor := get(marginalmode_options, "markercolor", :dimgrey)
            markersize := get(marginalmode_options, "markersize", 4)
            markershape := get(marginalmode_options, "markershape", :rect)
            markeralpha := get(marginalmode_options, "markeralpha", 1)
            linealpha := 0
            markerstrokealpha := get(marginalmode_options, "markerstrokealpha", 1)
            markerstrokecolor := get(marginalmode_options, "markerstrokecolor", :dimgrey)
            markerstrokestyle := get(marginalmode_options, "markerstrokestyle", :solid)
            markerstrokewidth := get(marginalmode_options, "markerstrokewidth", 1)
            colorbar := colorbar

            ([marginalmode_values[i][1]], [marginalmode_values[i][2]])
            end
        end
    end
end

ext: jl


text:
#=
This file is auto-generated. Do not edit.
=#

mutable struct HydroStorage <: HydroGen
    name::String
    available::Bool
    bus::Bus
    activepower::Float64
    reactivepower::Float64
    tech::TechHydro
    op_cost::TwoPartCost
    storagecapacity::Float64
    initial_storage::Float64
    internal::InfrastructureSystemsInternal
end

function HydroStorage(name, available, bus, activepower, reactivepower, tech, op_cost, storagecapacity, initial_storage, )
    HydroStorage(name, available, bus, activepower, reactivepower, tech, op_cost, storagecapacity, initial_storage, InfrastructureSystemsInternal())
end

function HydroStorage(; name, available, bus, activepower, reactivepower, tech, op_cost, storagecapacity, initial_storage, )
    HydroStorage(name, available, bus, activepower, reactivepower, tech, op_cost, storagecapacity, initial_storage, )
end

# Constructor for demo purposes; non-functional.

function HydroStorage(::Nothing)
    HydroStorage(;
        name="init",
        available=false,
        bus=Bus(nothing),
        activepower=0.0,
        reactivepower=0.0,
        tech=TechHydro(nothing),
        op_cost=TwoPartCost(nothing),
        storagecapacity=0.0,
        initial_storage=0.0,
    )
end

"""Get HydroStorage name."""
get_name(value::HydroStorage) = value.name
"""Get HydroStorage available."""
get_available(value::HydroStorage) = value.available
"""Get HydroStorage bus."""
get_bus(value::HydroStorage) = value.bus
"""Get HydroStorage activepower."""
get_activepower(value::HydroStorage) = value.activepower
"""Get HydroStorage reactivepower."""
get_reactivepower(value::HydroStorage) = value.reactivepower
"""Get HydroStorage tech."""
get_tech(value::HydroStorage) = value.tech
"""Get HydroStorage op_cost."""
get_op_cost(value::HydroStorage) = value.op_cost
"""Get HydroStorage storagecapacity."""
get_storagecapacity(value::HydroStorage) = value.storagecapacity
"""Get HydroStorage initial_storage."""
get_initial_storage(value::HydroStorage) = value.initial_storage
"""Get HydroStorage internal."""
get_internal(value::HydroStorage) = value.internal


text:
using oneAPI
using ExaTronKernels
using LinearAlgebra
using Random
using Test

try
    tmp = oneArray{Float64}(undef, 10)
catch e
    throw(e)
end

"""
Test ExaTron's internal routines written for GPU.

The current implementation assumes the following:
  - A thread block takes a matrix structure, (tx,ty), with
    n <= blockDim().x = blockDim().y <= 32 and blockDim().x is even.
  - Arguments passed on to these routines are assumed to be
    of size at least n. This is to prevent multiple thread
    divergence when we call a function with n_hat < n.
    Such a case occurs when we fix active variables.

We test the following routines, where [O] indicates if the routine
is checked if n < blockDim().x is OK.
  - dicf     [O][O]: this routine also tests dnsol and dtsol.
  - dicfs    [O][T]
  - dcauchy  [O][T]
  - dtrpcg   [O][T]
  - dprsrch  [O][T]
  - daxpy    [O][O]
  - dssyax   [O][O]: we do shuffle using blockDim().x.
  - dmid     [O][O]: we could use a single thread only to multiple divergences.
  - dgpstep  [O][O]
  - dbreakpt [O][O]: we use the existing ExaTron implementation as is.
  - dnrm2    [O][O]: we do shuffle using blockDim().x.
  - nrm2     [O][O]: we do shuffle using blockDim().x.
  - dcopy    [O][O]
  - ddot     [O][O]
  - dscal    [O][O]
  - dtrqsol  [O][O]
  - dspcg    [O][T]: we use a single thread to avoid multiple divergences.
  - dgpnorm  [O][O]
  - dtron    [O]
  - driver_kernel
"""

Random.seed!(0)
itermax = 10
n = 4
nblk = 4

@testset "dicf" begin
println("Testing dicf")
    function dicf_test(
        ::Val{n},
        d_in::oneDeviceArray{Float64, 2},
        d_out::oneDeviceArray{Float64, 2}
    ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        L = oneLocalArray(Float64, (n,n))

        for i in 1:n
            L[i,tx] = d_in[i,tx]
        end
        barrier()

        # Test Cholesky factorization.
        ExaTronKernels.dicf(n,L)
        if bx == 1
            for i in 1:n
                d_out[tx,i] = L[tx,i]
            end
        end
        barrier()
    end

    for i=1:itermax
        L = tril(rand(n,n))
        A = L*transpose(L)
        A .= tril(A) .+ (transpose(tril(A)) .- Diagonal(A))
        tron_A = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
        tron_A.vals .= A

        d_in = oneArray{Float64,2}(undef, (n,n))
        d_out = oneArray{Float64,2}(undef, (n,n))
        copyto!(d_in, tron_A.vals)
        @oneapi items=(n, nblk) groups=nblk dicf_test(Val{n}(), d_in, d_out)
        h_L = d_out |> Array

        tron_L = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
        tron_L.vals .= tron_A.vals
        indr = zeros(Int, n)
        indf = zeros(n)
        list = zeros(n)
        w = zeros(n)
        ExaTronKernels.dicf(n, n^2, tron_L, 5, indr, indf, list, w)

        @test norm(tron_A.vals .- tril(h_L)*transpose(tril(h_L))) <= 1e-10
        @test norm(tril(h_L) .- transpose(triu(h_L))) <= 1e-10
        @test norm(tril(tron_L.vals) .- tril(h_L)) <= 1e-9
    end
end

@testset "dicfs" begin
    println("Testing dicfs")
    function dicfs_test(::Val{n},
        alpha::Float64,
        dA::oneDeviceArray{Float64, 2},
        d_out::oneDeviceArray{Float64, 2}
    ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        wa1 = oneLocalArray(Float64, n)
        wa2 = oneLocalArray(Float64, n)
        A = oneLocalArray(Float64, (n,n))
        L = oneLocalArray(Float64, (n,n))

        for i in 1:n
            A[i,tx] = dA[i,tx]
        end
        barrier()

        ExaTronKernels.dicfs(n, alpha, A, L, wa1, wa2)
        if bx == 1
            for i in 1:n
                d_out[tx,i] = L[tx,i]
            end
        end
        barrier()
    end

    for i=1:itermax
        L = tril(rand(n,n))
        A = L*transpose(L)
        A .= tril(A) .+ (transpose(tril(A)) .- Diagonal(A))
        tron_A = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
        tron_L = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
        tron_A.vals .= A

        dA = oneArray{Float64,2}(undef, (n,n))
        d_out = oneArray{Float64,2}(undef, (n,n))
        alpha = 1.0
        copyto!(dA, tron_A.vals)
        @oneapi items=(n, nblk) groups=nblk dicfs_test(Val{n}(),alpha,dA,d_out)
        h_L = zeros(n,n)
        copyto!(h_L, d_out)
        iwa = zeros(Int, 3*n)
        wa1 = zeros(n)
        wa2 = zeros(n)
        ExaTronKernels.dicfs(n, n^2, tron_A, tron_L, 5, alpha, iwa, wa1, wa2)

        @test norm(tril(h_L) .- transpose(triu(h_L))) <= 1e-10
        @test norm(tril(tron_L.vals) .- tril(h_L)) <= 1e-9

        # Make it negative definite.
        for j=1:n
            tron_A.vals[j,j] = -tron_A.vals[j,j]
        end
        copyto!(dA, tron_A.vals)
        @oneapi items=(n, nblk) groups=nblk dicfs_test(Val{n}(),alpha,dA,d_out)
        copyto!(h_L, d_out)
        ExaTronKernels.dicfs(n, n^2, tron_A, tron_L, 5, alpha, iwa, wa1, wa2)

        @test norm(tril(h_L) .- transpose(triu(h_L))) <= 1e-10
        @test norm(tril(tron_L.vals) .- tril(h_L)) <= 1e-10
    end
end

@testset "dcauchy" begin
    println("Testing dcauchy")
    function dcauchy_test(::Val{n},
                            dx::oneDeviceArray{Float64, 1},
                            dl::oneDeviceArray{Float64, 1},
                            du::oneDeviceArray{Float64, 1},
                            dA::oneDeviceArray{Float64, 2},
                            dg::oneDeviceArray{Float64, 1},
                            delta::Float64,
                            alpha::Float64,
                            d_out1::oneDeviceArray{Float64, 1},
                            d_out2::oneDeviceArray{Float64, 1}
                            ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        xl = oneLocalArray(Float64, n)
        xu = oneLocalArray(Float64, n)
        g = oneLocalArray(Float64, n)
        s = oneLocalArray(Float64, n)
        wa = oneLocalArray(Float64, n)
        A = oneLocalArray(Float64, (n,n))

        for i in 1:n
            A[i,tx] = dA[i,tx]
        end
        x[tx] = dx[tx]
        xl[tx] = dl[tx]
        xu[tx] = du[tx]
        g[tx] = dg[tx]

        alpha = ExaTronKernels.dcauchy(n,x,xl,xu,A,g,delta,alpha,s,wa)
        if bx == 1
            d_out1[tx] = s[tx]
            d_out2[tx] = alpha
        end
        barrier()
    end

    for i=1:itermax
        L = tril(rand(n,n))
        A = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
        A.vals .= L*transpose(L)
        A.vals .= tril(A.vals) .+ (transpose(tril(A.vals)) .- Diagonal(A.vals))
        x = rand(n)
        xl = x .- abs.(rand(n))
        xu = x .+ abs.(rand(n))
        g = A.vals*x .+ rand(n)
        s = zeros(n)
        wa = zeros(n)
        alpha = 1.0
        delta = 2.0*norm(g)

        dx = oneArray{Float64}(undef, n)
        dl = oneArray{Float64}(undef, n)
        du = oneArray{Float64}(undef, n)
        dg = oneArray{Float64}(undef, n)
        dA = oneArray{Float64,2}(undef, (n,n))
        d_out1 = oneArray{Float64}(undef, n)
        d_out2 = oneArray{Float64}(undef, n)
        copyto!(dx, x)
        copyto!(dl, xl)
        copyto!(du, xu)
        copyto!(dg, g)
        copyto!(dA, A.vals)
        @oneapi items=(n, nblk) groups=nblk dcauchy_test(Val{n}(),dx,dl,du,dA,dg,delta,alpha,d_out1,d_out2)
        h_s = zeros(n)
        h_alpha = zeros(n)
        copyto!(h_s, d_out1)
        copyto!(h_alpha, d_out2)

        alpha = ExaTronKernels.dcauchy(n, x, xl, xu, A, g, delta, alpha, s, wa)

        @test norm(s .- h_s) <= 1e-10
        @test norm(alpha .- h_alpha) <= 1e-10
    end
end

@testset "dtrpcg" begin
    println("Testing dtrpcg")
    function dtrpcg_test(::Val{n}, delta::Float64, tol::Float64,
                            stol::Float64,
                            d_in::oneDeviceArray{Float64, 2},
                            d_g::oneDeviceArray{Float64, 1},
                            d_out_L::oneDeviceArray{Float64, 2},
                            d_out::oneDeviceArray{Float64, 1}
                            ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        A = oneLocalArray(Float64, (n,n))
        L = oneLocalArray(Float64, (n,n))

        g = oneLocalArray(Float64, n)
        w = oneLocalArray(Float64, n)
        p = oneLocalArray(Float64, n)
        q = oneLocalArray(Float64, n)
        r = oneLocalArray(Float64, n)
        t = oneLocalArray(Float64, n)
        z = oneLocalArray(Float64, n)

        for i in 1:n
            A[i,tx] = d_in[i,tx]
            L[i,tx] = d_in[i,tx]
        end
        g[tx] = d_g[tx]
        barrier()

        ExaTronKernels.dicf(n,L)
        info, iters = ExaTronKernels.dtrpcg(n,A,g,delta,L,tol,stol,n,w,p,q,r,t,z)
        if bx == 1
            d_out[tx] = w[tx]
            for i in 1:n
                d_out_L[i,tx] = L[i,tx]
            end
        end
        barrier()
    end

    delta = 100.0
    tol = 1e-6
    stol = 1e-6
    tron_A = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
    tron_L = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
    for i=1:itermax
        L = tril(rand(n,n))
        A = L*transpose(L)
        A .= tril(A) .+ (transpose(tril(A)) .- Diagonal(A))
        g = 0.1*ones(n)
        w = zeros(n)
        p = zeros(n)
        q = zeros(n)
        r = zeros(n)
        t = zeros(n)
        z = zeros(n)
        tron_A.vals .= A
        tron_L.vals .= A
        d_in = oneArray{Float64,2}(undef, (n,n))
        d_g = oneArray{Float64}(undef, n)
        d_out_L = oneArray{Float64,2}(undef, (n,n))
        d_out = oneArray{Float64}(undef, n)
        copyto!(d_in, A)
        copyto!(d_g, g)
        @oneapi items=(n, nblk) groups=nblk dtrpcg_test(Val{n}(),delta,tol,stol,d_in,d_g,d_out_L,d_out)
        oneAPI.synchronize()
        h_w = zeros(n)
        h_L = zeros(n,n)
        copyto!(h_L, d_out_L)
        copyto!(h_w, d_out)

        indr = zeros(Int, n)
        indf = zeros(n)
        list = zeros(n)
        ExaTronKernels.dicf(n, n^2, tron_L, 5, indr, indf, list, w)
        ExaTronKernels.dtrpcg(n, tron_A, g, delta, tron_L, tol, stol, n, w, p, q, r, t, z)

        @test norm(tril(h_L) .- tril(tron_L.vals)) <= tol
        @test norm(h_w .- w) <= tol
    end
end

@testset "dprsrch" begin
    println("Testing dprsrch")
    function dprsrch_test(::Val{n},
                            d_x::oneDeviceArray{Float64, 1},
                            d_xl::oneDeviceArray{Float64, 1},
                            d_xu::oneDeviceArray{Float64, 1},
                            d_g::oneDeviceArray{Float64, 1},
                            d_w::oneDeviceArray{Float64, 1},
                            d_A::oneDeviceArray{Float64, 2},
                            d_out1::oneDeviceArray{Float64, 1},
                            d_out2::oneDeviceArray{Float64, 1}
                            ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        xl = oneLocalArray(Float64, n)
        xu = oneLocalArray(Float64, n)
        g = oneLocalArray(Float64, n)
        w = oneLocalArray(Float64, n)
        wa1 = oneLocalArray(Float64, n)
        wa2 = oneLocalArray(Float64, n)
        A = oneLocalArray(Float64, (n,n))

        for i in 1:n
            A[i,tx] = d_A[i,tx]
        end
        x[tx] = d_x[tx]
        xl[tx] = d_xl[tx]
        xu[tx] = d_xu[tx]
        g[tx] = d_g[tx]
        w[tx] = d_w[tx]
        barrier()

        ExaTronKernels.dprsrch(n, x, xl, xu, A, g, w, wa1, wa2)
        if bx == 1
            d_out1[tx] = x[tx]
            d_out2[tx] = w[tx]
        end
        barrier()

    end

    for i=1:itermax
        L = tril(rand(n,n))
        A = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
        A.vals .= L*transpose(L)
        A.vals .= tril(A.vals) .+ (transpose(tril(A.vals)) .- Diagonal(A.vals))
        x = rand(n)
        xl = x .- abs.(rand(n))
        xu = x .+ abs.(rand(n))
        g = A.vals*x .+ rand(n)
        w = -g
        wa1 = zeros(n)
        wa2 = zeros(n)

        dx = oneArray{Float64}(undef, n)
        dl = oneArray{Float64}(undef, n)
        du = oneArray{Float64}(undef, n)
        dg = oneArray{Float64}(undef, n)
        dw = oneArray{Float64}(undef, n)
        dA = oneArray{Float64,2}(undef, (n,n))
        d_out1 = oneArray{Float64}(undef, n)
        d_out2 = oneArray{Float64}(undef, n)
        copyto!(dx, x)
        copyto!(dl, xl)
        copyto!(du, xu)
        copyto!(dg, g)
        copyto!(dw, w)
        copyto!(dA, A.vals)
        @oneapi items=(n, nblk) groups=nblk dprsrch_test(Val{n}(),dx,dl,du,dg,dw,dA,d_out1,d_out2)
        h_x = zeros(n)
        h_w = zeros(n)
        copyto!(h_x, d_out1)
        copyto!(h_w, d_out2)

        ExaTronKernels.dprsrch(n,x,xl,xu,A,g,w,wa1,wa2)

        @test norm(x .- h_x) <= 1e-10
        @test norm(w .- h_w) <= 1e-10
    end
end

@testset "daxpy" begin
    println("Testing daxpy")
    function daxpy_test(::Val{n},
                        da,
                        d_in::oneDeviceArray{Float64, 1},
                        d_out::oneDeviceArray{Float64, 1}
                        ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        y = oneLocalArray(Float64, n)

        x[tx] = d_in[tx]
        y[tx] = d_in[tx + n]
        barrier()

        ExaTronKernels.daxpy(n, da, x, 1, y, 1)
        if bx == 1
            d_out[tx] = y[tx]
        end
        barrier()

    end

    for i=1:itermax
        da = rand(1)[1]
        h_in = rand(2*n)
        h_out = zeros(n)
        d_in = oneArray{Float64}(undef, 2*n)
        d_out = oneArray{Float64}(undef, n)
        copyto!(d_in, h_in)
        @oneapi items=(n, nblk) groups=nblk daxpy_test(Val{n}(),da,d_in,d_out)
        copyto!(h_out, d_out)

        @test norm(h_out .- (h_in[n+1:2*n] .+ da.*h_in[1:n])) <= 1e-12
    end
end

@testset "dssyax" begin
    println("Testing dssyax")
    function dssyax_test(::Val{n},
                            d_z::oneDeviceArray{Float64, 1},
                            d_in::oneDeviceArray{Float64, 2},
                            d_out::oneDeviceArray{Float64, 1}
                        ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        z = oneLocalArray(Float64, n)
        q = oneLocalArray(Float64, n)
        A = oneLocalArray(Float64, (n,n))

        for i in 1:n
            A[i,tx] = d_in[i,tx]
        end
        z[tx] = d_z[tx]
        barrier()

        ExaTronKernels.dssyax(n, A, z, q)
        if bx == 1
            d_out[tx] = q[tx]
        end
        barrier()

    end

    for i=1:itermax
        z = rand(n)
        h_in = rand(n,n)
        h_out = zeros(n)
        d_z = oneArray{Float64}(undef, n)
        d_in = oneArray{Float64,2}(undef, (n,n))
        d_out = oneArray{Float64}(undef, n)
        copyto!(d_z, z)
        copyto!(d_in, h_in)
        @oneapi items=(n, nblk) groups=nblk dssyax_test(Val{n}(),d_z,d_in,d_out)
        copyto!(h_out, d_out)
        @test norm(h_out .- h_in*z) <= 1e-12
    end
end

@testset "dmid" begin
    println("Testing dmid")
    function dmid_test(::Val{n},
                        dx::oneDeviceArray{Float64, 1},
                        dl::oneDeviceArray{Float64, 1},
                        du::oneDeviceArray{Float64, 1},
                        d_out::oneDeviceArray{Float64, 1}
                    ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        xl = oneLocalArray(Float64, n)
        xu = oneLocalArray(Float64, n)

        x[tx] = dx[tx]
        xl[tx] = dl[tx]
        xu[tx] = du[tx]
        barrier()

        ExaTronKernels.dmid(n, x, xl, xu)
        if bx == 1
            d_out[tx] = x[tx]
        end
        barrier()

    end

    for i=1:itermax
        x = rand(n)
        xl = x .- abs.(rand(n))
        xu = x .+ abs.(rand(n))

        # Force some components to go below or above bounds
        # so that we can test all cases.
        for j=1:n
            k = rand(1:3)
            if k == 1
                x[j] = xl[j] - 0.1
            elseif k == 2
                x[j] = xu[j] + 0.1
            end
        end
        x_out = zeros(n)
        dx = oneArray{Float64}(undef, n)
        dl = oneArray{Float64}(undef, n)
        du = oneArray{Float64}(undef, n)
        d_out = oneArray{Float64}(undef, n)
        copyto!(dx, x)
        copyto!(dl, xl)
        copyto!(du, xu)
        @oneapi items=(n, nblk) groups=nblk dmid_test(Val{n}(),dx,dl,du,d_out)
        copyto!(x_out, d_out)

        ExaTronKernels.dmid(n, x, xl, xu)
        @test !(false in (x .== x_out))
    end
end

@testset "dgpstep" begin
    println("Testing dgpstep")
    function dgpstep_test(::Val{n},dx,
                            dl::oneDeviceArray{Float64, 1},
                            du::oneDeviceArray{Float64, 1},
                            alpha::Float64,
                            dw::oneDeviceArray{Float64, 1},
                            d_out::oneDeviceArray{Float64, 1}
                            ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        xl = oneLocalArray(Float64, n)
        xu = oneLocalArray(Float64, n)
        w = oneLocalArray(Float64, n)
        s = oneLocalArray(Float64, n)

        x[tx] = dx[tx]
        xl[tx] = dl[tx]
        xu[tx] = du[tx]
        w[tx] = dw[tx]
        barrier()

        ExaTronKernels.dgpstep(n, x, xl, xu, alpha, w, s)
        if bx == 1
            d_out[tx] = s[tx]
        end
        barrier()

    end

    for i=1:itermax
        x = rand(n)
        xl = x .- abs.(rand(n))
        xu = x .+ abs.(rand(n))
        w = rand(n)
        alpha = rand(1)[1]
        s = zeros(n)
        s_out = zeros(n)

        # Force some components to go below or above bounds
        # so that we can test all cases.
        for j=1:n
            k = rand(1:3)
            if k == 1
                if x[j] + alpha*w[j] >= xl[j]
                    w[j] = (xl[j] - x[j]) / alpha - 0.1
                end
            elseif k == 2
                if x[j] + alpha*w[j] <= xu[j]
                    w[j] = (xu[j] - x[j]) / alpha + 0.1
                end
            end
        end

        dx = oneArray{Float64}(undef, n)
        dl = oneArray{Float64}(undef, n)
        du = oneArray{Float64}(undef, n)
        dw = oneArray{Float64}(undef, n)
        d_out = oneArray{Float64}(undef, n)
        copyto!(dx, x)
        copyto!(dl, xl)
        copyto!(du, xu)
        copyto!(dw, w)
        @oneapi items=(n, nblk) groups=nblk dgpstep_test(Val{n}(),dx,dl,du,alpha,dw,d_out)
        copyto!(s_out, d_out)

        ExaTronKernels.dgpstep(n, x, xl, xu, alpha, w, s)
        @test !(false in (s .== s_out))
    end
end

@testset "dbreakpt" begin
    println("Testing dbreakpt")
    function dbreakpt_test(::Val{n},
                            dx::oneDeviceArray{Float64},
                            dl::oneDeviceArray{Float64},
                            du::oneDeviceArray{Float64},
                            dw::oneDeviceArray{Float64},
                            d_nbrpt::oneDeviceArray{Float64},
                            d_brptmin::oneDeviceArray{Float64},
                            d_brptmax::oneDeviceArray{Float64}
                            ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        xl = oneLocalArray(Float64, n)
        xu = oneLocalArray(Float64, n)
        w = oneLocalArray(Float64, n)
        x[tx] = dx[tx]
        xl[tx] = dl[tx]
        xu[tx] = du[tx]
        w[tx] = dw[tx]
        barrier()

        nbrpt, brptmin, brptmax = ExaTronKernels.dbreakpt(n, x, xl, xu, w)
        if bx == 1
            for i in 1:n
                d_nbrpt[i,tx] = nbrpt
                d_brptmin[i,tx] = brptmin
                d_brptmax[i,tx] = brptmax
            end
        end
        barrier()

    end

    for i=1:itermax
        x = rand(n)
        xl = x .- abs.(rand(n))
        xu = x .+ abs.(rand(n))
        w = 2.0*rand(n) .- 1.0     # (-1,1]
        h_nbrpt = zeros((n,n))
        h_brptmin = zeros((n,n))
        h_brptmax = zeros((n,n))

        dx = oneArray{Float64}(undef, n)
        dl = oneArray{Float64}(undef, n)
        du = oneArray{Float64}(undef, n)
        dw = oneArray{Float64}(undef, n)
        d_nbrpt = oneArray{Float64,2}(undef, (n,n))
        d_brptmin = oneArray{Float64,2}(undef, (n,n))
        d_brptmax = oneArray{Float64,2}(undef, (n,n))
        copyto!(dx, x)
        copyto!(dl, xl)
        copyto!(du, xu)
        copyto!(dw, w)
        @oneapi items=(n, nblk) groups=nblk dbreakpt_test(Val{n}(),dx,dl,du,dw,d_nbrpt,d_brptmin,d_brptmax)
        copyto!(h_nbrpt, d_nbrpt)
        copyto!(h_brptmin, d_brptmin)
        copyto!(h_brptmax, d_brptmax)

        nbrpt, brptmin, brptmax = ExaTronKernels.dbreakpt(n, x, xl, xu, w)
        @test !(false in (nbrpt .== h_nbrpt))
        @test !(false in (brptmin .== h_brptmin))
        @test !(false in (brptmax .== h_brptmax))
    end
end

@testset "dnrm2" begin
    println("Testing dnrm2")
    function dnrm2_test(::Val{n},
                            d_in::oneDeviceArray{Float64, 1},
                            d_out::oneDeviceArray{Float64, 1}
                        ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        x[tx] = d_in[tx]
        barrier()

        v = ExaTronKernels.dnrm2(n, x, 1)
        if bx == 1
            d_out[tx] = v
        end
        barrier()

    end

    @inbounds for i=1:itermax
        h_in = rand(n)
        h_out = zeros(n)
        d_in = oneArray{Float64}(undef, n)
        d_out = oneArray{Float64}(undef, n)
        copyto!(d_in, h_in)
        @oneapi items=(n, nblk) groups=nblk dnrm2_test(Val{n}(),d_in,d_out)
        copyto!(h_out, d_out)
        xnorm = norm(h_in, 2)

        @test norm(xnorm .- h_out) <= 1e-10
    end
end

@testset "nrm2" begin
    println("Testing nrm2")
    function nrm2_test(::Val{n}, d_A::oneDeviceMatrix{Float64}, d_out::oneDeviceVector{Float64}) where {n}
        tx = get_local_id()
        bx = get_group_id()

        wa = oneLocalArray(Float64, n)
        A = oneLocalArray(Float64, (n,n))

        for i in 1:n
            A[i,tx] = d_A[i,tx]
        end
        barrier()

        ExaTronKernels.nrm2!(wa, A, n)
        if bx == 1
            d_out[tx] = wa[tx]
        end
        barrier()

        return
    end

    @inbounds for i=1:itermax
        L = tril(rand(n,n))
        A = L*transpose(L)
        A .= tril(A) .+ (transpose(tril(A)) .- Diagonal(A))
        wa = zeros(n)
        tron_A = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
        tron_A.vals .= A
        ExaTronKernels.nrm2!(wa, tron_A, n)

        d_A = oneArray{Float64,2}(undef, (n,n))
        d_out = oneArray{Float64}(undef, n)
        h_wa = zeros(n)
        copyto!(d_A, A)
        @oneapi items=(n, nblk) groups=nblk nrm2_test(Val{n}(),d_A,d_out)
        copyto!(h_wa, d_out)

        @test norm(wa .- h_wa) <= 1e-10
    end
end

@testset "dcopy" begin
    println("Testing dcopy")
    function dcopy_test(::Val{n}, d_in::oneDeviceVector{Float64, 1},
                        d_out::oneDeviceVector{Float64, 1}) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        y = oneLocalArray(Float64, n)

        x[tx] = d_in[tx]
        barrier()

        ExaTronKernels.dcopy(n, x, 1, y, 1)

        if bx == 1
            d_out[tx] = y[tx]
        end
        barrier()

        return
    end

    @inbounds for i=1:itermax
        h_in = rand(n)
        h_out = zeros(n)
        d_in = oneArray{Float64}(undef, n)
        d_out = oneArray{Float64}(undef, n)
        copyto!(d_in, h_in)
        @oneapi items=(n,nblk) groups=nblk dcopy_test(Val{n}(),d_in,d_out)
        oneAPI.synchronize()
        copyto!(h_out, d_out)

        @test !(false in (h_in .== h_out))
    end
end

@testset "ddot" begin
    println("Testing ddot")
    function ddot_test(
                        ::Val{n},
                        d_in::oneDeviceArray{Float64, 1},
                        d_out::oneDeviceArray{Float64, 2}
                    ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        y = oneLocalArray(Float64, n)
        x[tx] = d_in[tx]
        y[tx] = d_in[tx]
        barrier()

        v = ExaTronKernels.ddot(n, x, 1, y, 1)

        if bx == 1
            for i in 1:n
                d_out[i,tx] = v
            end
        end
        barrier()

    end

    @inbounds for i=1:itermax
        h_in = rand(n)
        h_out = zeros((n,n))
        d_in = oneArray{Float64}(undef, n)
        d_out = oneArray{Float64,2}(undef, (n,n))
        copyto!(d_in, h_in)
        @oneapi items=(n, nblk) groups=nblk ddot_test(Val{n}(),d_in,d_out)
        copyto!(h_out, d_out)

        @test norm(dot(h_in,h_in) .- h_out, 2) <= 1e-10
    end
end

@testset "dscal" begin
    println("Testing dscal")
    function dscal_test(::Val{n}, da::Float64,
                        d_in::oneDeviceVector{Float64, 1},
                        d_out::oneDeviceVector{Float64, 1}
                        ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        x[tx] = d_in[tx]
        barrier()

        ExaTronKernels.dscal(n, da, x, 1)
        if bx == 1
            d_out[tx] = x[tx]
        end
        barrier()

    end

    for i=1:itermax
        h_in = rand(n)
        h_out = zeros(n)
        da = rand(1)[1]
        d_in = oneArray{Float64}(undef, n)
        d_out = oneArray{Float64}(undef, n)
        copyto!(d_in, h_in)
        @oneapi items=(n, nblk) groups=nblk dscal_test(Val{n}(),da,d_in,d_out)
        copyto!(h_out, d_out)

        @test norm(h_out .- (da.*h_in)) <= 1e-12
    end
end

@testset "dtrqsol" begin
    println("Testing dtrqsol")
    function dtrqsol_test(::Val{n}, d_x,
                            d_p,
                            d_out,
                            delta::Float64
                            ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        p = oneLocalArray(Float64, n)

        x[tx] = d_x[tx]
        p[tx] = d_p[tx]
        barrier()

        sigma = ExaTronKernels.dtrqsol(n, x, p, delta)
        if bx == 1
            for i in 1:n
                d_out[i,tx] = sigma
            end
        end
        barrier()
    end

    for i=1:itermax
        x = rand(n)
        p = rand(n)
        sigma = abs(rand(1)[1])
        delta = norm(x .+ sigma.*p)

        d_x = oneArray{Float64}(undef, n)
        d_p = oneArray{Float64}(undef, n)
        d_out = oneArray{Float64,2}(undef, (n,n))
        copyto!(d_x, x)
        copyto!(d_p, p)
        @oneapi items=(n, nblk) groups=nblk dtrqsol_test(Val{n}(),d_x,d_p,d_out,delta)

        d_out = d_out |> Array
        @test norm(sigma .- d_out) <= 1e-10
    end
end

@testset "dspcg" begin
    println("Testing dspcg")
    function dspcg_test(::Val{n}, delta::Float64, rtol::Float64,
                        cg_itermax::Int64,
                        dx::oneDeviceArray{Float64, 1},
                        dxl::oneDeviceArray{Float64, 1},
                        dxu::oneDeviceArray{Float64, 1},
                        dA::oneDeviceArray{Float64, 2},
                        dg::oneDeviceArray{Float64, 1},
                        ds::oneDeviceArray{Float64, 1},
                        d_out::oneDeviceArray{Float64, 1}
                        ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        xl = oneLocalArray(Float64, n)
        xu = oneLocalArray(Float64, n)
        g = oneLocalArray(Float64, n)
        s = oneLocalArray(Float64, n)
        w = oneLocalArray(Float64, n)
        wa1 = oneLocalArray(Float64, n)
        wa2 = oneLocalArray(Float64, n)
        wa3 = oneLocalArray(Float64, n)
        wa4 = oneLocalArray(Float64, n)
        wa5 = oneLocalArray(Float64, n)
        gfree = oneLocalArray(Float64, n)
        indfree = oneLocalArray(Int64, n)
        iwa = oneLocalArray(Int64, 2*n)

        A = oneLocalArray(Float64, (n,n))
        B = oneLocalArray(Float64, (n,n))
        L = oneLocalArray(Float64, (n,n))

        @inbounds for j=1:n
            A[j,tx] = dA[j,tx]
        end
        x[tx] = dx[tx]
        xl[tx] = dxl[tx]
        xu[tx] = dxu[tx]
        g[tx] = dg[tx]
        s[tx] = ds[tx]
        barrier()

        ExaTronKernels.dspcg(n, delta, rtol, cg_itermax, x, xl, xu,
                        A, g, s, B, L, indfree, gfree, w, iwa,
                        wa1, wa2, wa3, wa4, wa5)

        if bx == 1
            d_out[tx] = x[tx]
        end
        barrier()

    end

    for i=1:itermax
        L = tril(rand(n,n))
        A = L*transpose(L)
        A .= tril(A) .+ (transpose(tril(A)) .- Diagonal(A))
        tron_A = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
        tron_A.vals .= A
        tron_B = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
        tron_L = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
        x = rand(n)
        xl = x .- abs.(rand(n))
        xu = x .+ abs.(rand(n))
        g = A*x .+ rand(n)
        s = rand(n)
        delta = 2.0*norm(g)
        rtol = 1e-6
        cg_itermax = n
        w = zeros(n)
        wa = zeros(5*n)
        gfree = zeros(n)
        indfree = zeros(Int64, n)
        iwa = zeros(Int64, 3*n)

        dx = oneArray{Float64}(undef, n)
        dxl = oneArray{Float64}(undef, n)
        dxu = oneArray{Float64}(undef, n)
        dA = oneArray{Float64,2}(undef, (n,n))
        dg = oneArray{Float64}(undef, n)
        ds = oneArray{Float64}(undef, n)
        d_out = oneArray{Float64}(undef, n)

        copyto!(dx, x)
        copyto!(dxl, xl)
        copyto!(dxu, xu)
        copyto!(dA, tron_A.vals)
        copyto!(dg, g)
        copyto!(ds, s)

        @oneapi items=(n, nblk) groups=nblk dspcg_test(Val{n}(),delta,rtol,cg_itermax,dx,dxl,dxu,dA,dg,ds,d_out)
        h_x = zeros(n)
        copyto!(h_x, d_out)

        ExaTronKernels.dspcg(n, x, xl, xu, tron_A, g, delta, rtol, s, 5, cg_itermax,
                        tron_B, tron_L, indfree, gfree, w, wa, iwa)

        @test norm(x .- h_x) <= 1e-10
    end
end

@testset "dgpnorm" begin
    println("Testing dgpnorm")
    function dgpnorm_test(::Val{n},
        dx::oneDeviceArray{Float64, 1},
        dxl::oneDeviceArray{Float64, 1},
        dxu::oneDeviceArray{Float64, 1},
        dg::oneDeviceArray{Float64, 1},
        d_out::oneDeviceArray{Float64, 1}
    ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        xl = oneLocalArray(Float64, n)
        xu = oneLocalArray(Float64, n)
        g = oneLocalArray(Float64, n)

        x[tx] = dx[tx]
        xl[tx] = dxl[tx]
        xu[tx] = dxu[tx]
        g[tx] = dg[tx]
        barrier()

        v = ExaTronKernels.dgpnorm(n, x, xl, xu, g)
        if bx == 1
            d_out[tx] = v
        end
        barrier()

    end

    for i=1:itermax
        x = rand(n)
        xl = x .- abs.(rand(n))
        xu = x .+ abs.(rand(n))
        g = 2.0*rand(n) .- 1.0

        dx = oneArray{Float64}(undef, n)
        dxl = oneArray{Float64}(undef, n)
        dxu = oneArray{Float64}(undef, n)
        dg = oneArray{Float64}(undef, n)
        d_out = oneArray{Float64}(undef, n)

        copyto!(dx, x)
        copyto!(dxl, xl)
        copyto!(dxu, xu)
        copyto!(dg, g)

        @oneapi items=(n, nblk) groups=nblk dgpnorm_test(Val{n}(), dx, dxl, dxu, dg, d_out)
        h_v = zeros(n)
        copyto!(h_v, d_out)

        v = ExaTronKernels.dgpnorm(n, x, xl, xu, g)
        @test norm(h_v .- v) <= 1e-10
    end
end

@testset "dtron" begin
    println("Testing dtron")
    function dtron_test(::Val{n}, f::Float64, frtol::Float64, fatol::Float64, fmin::Float64,
                        cgtol::Float64, cg_itermax::Int, delta::Float64, task::Int,
                        disave::oneDeviceArray{Int}, ddsave::oneDeviceArray{Float64, 1},
                        dx::oneDeviceArray{Float64, 1}, dxl::oneDeviceArray{Float64, 1},
                        dxu::oneDeviceArray{Float64, 1}, dA::oneDeviceArray{Float64, 2},
                        dg::oneDeviceArray{Float64, 1}, d_out::oneDeviceArray{Float64, 1}
                        ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        xl = oneLocalArray(Float64, n)
        xu = oneLocalArray(Float64, n)
        g = oneLocalArray(Float64, n)
        xc = oneLocalArray(Float64, n)
        s = oneLocalArray(Float64, n)
        wa = oneLocalArray(Float64, n)
        wa1 = oneLocalArray(Float64, n)
        wa2 = oneLocalArray(Float64, n)
        wa3 = oneLocalArray(Float64, n)
        wa4 = oneLocalArray(Float64, n)
        wa5 = oneLocalArray(Float64, n)
        gfree = oneLocalArray(Float64, n)
        indfree = oneLocalArray(Int, n)
        iwa = oneLocalArray(Int, n)

        A = oneLocalArray(Float64, (n,n))
        B = oneLocalArray(Float64, (n,n))
        L = oneLocalArray(Float64, (n,n))

        @inbounds for j=1:n
            A[j,tx] = dA[j,tx]
        end
        x[tx] = dx[tx]
        xl[tx] = dxl[tx]
        xu[tx] = dxu[tx]
        g[tx] = dg[tx]
        barrier()

        ExaTronKernels.dtron(n, x, xl, xu, f, g, A, frtol, fatol, fmin, cgtol,
                        cg_itermax, delta, task, B, L, xc, s, indfree, gfree,
                        disave, ddsave, wa, iwa, wa1, wa2, wa3, wa4, wa5)
        if bx == 1
            d_out[tx] = x[tx]
        end
        barrier()

        return nothing
    end

    for i=1:itermax
        L = tril(rand(n,n))
        A = L*transpose(L)
        A .= tril(A) .+ (transpose(tril(A)) .- Diagonal(A))
        x = rand(n)
        xl = x .- abs.(rand(n))
        xu = x .+ abs.(rand(n))
        c = rand(n)
        g = A*x .+ c
        xc = zeros(n)
        s = zeros(n)
        wa = zeros(7*n)
        gfree = zeros(n)
        indfree = zeros(Int, n)
        iwa = zeros(Int, 3*n)
        isave = zeros(Int, 3)
        dsave = zeros(3)
        task = 0
        fatol = 0.0
        frtol = 1e-12
        fmin = -1e32
        cgtol = 0.1
        cg_itermax = n
        delta = 2.0*norm(g)
        f = 0.5*transpose(x)*A*x .+ transpose(x)*c

        tron_A = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
        tron_A.vals .= A
        tron_B = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
        tron_L = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)

        dx = oneArray{Float64}(undef, n)
        dxl = oneArray{Float64}(undef, n)
        dxu = oneArray{Float64}(undef, n)
        dA = oneArray{Float64,2}(undef, (n,n))
        dg = oneArray{Float64}(undef, n)
        disave = oneArray{Int}(undef, n)
        ddsave = oneArray{Float64}(undef, n)
        d_out = oneArray{Float64}(undef, n)

        copyto!(dx, x)
        copyto!(dxl, xl)
        copyto!(dxu, xu)
        copyto!(dA, tron_A.vals)
        copyto!(dg, g)

        @oneapi items=(n, nblk) groups=nblk dtron_test(Val{n}(),f,frtol,fatol,fmin,cgtol,cg_itermax,delta,task,disave,ddsave,dx,dxl,dxu,dA,dg,d_out)
        h_x = zeros(n)
        copyto!(h_x, d_out)

        task_str = Vector{UInt8}(undef, 60)
        for (i,s) in enumerate("START")
            task_str[i] = UInt8(s)
        end

        ExaTronKernels.dtron(n, x, xl, xu, f, g, tron_A, frtol, fatol, fmin, cgtol,
                        cg_itermax, delta, task_str, tron_B, tron_L, xc, s, indfree,
                        isave, dsave, wa, iwa)
        @test norm(x .- h_x) <= 1e-10
    end
end

@testset "driver_kernel" begin
    println("Testing driver")
    @inline function eval_f(n, x, dA, dc)
        f = 0.0
        @inbounds for i=1:n
            @inbounds for j=1:n
                f += x[i]*dA[i,j]*x[j]
            end
        end
        f = 0.5*f
        @inbounds for i=1:n
            f += x[i]*dc[i]
        end
        barrier()
        return f
    end

    @inline function eval_g(n, x, g, dA, dc)
        @inbounds for i=1:n
            gval = 0.0
            @inbounds for j=1:n
                gval += dA[i,j]*x[j]
            end
            g[i] = gval + dc[i]
        end
        barrier()
        return
    end

    @inline function eval_h(n, scale, x, A, dA)
        tx = get_local_id()

        @inbounds for j=1:n
            A[j,tx] = dA[j,tx]
        end
        barrier()
        return
    end

    @inline function driver_kernel(N::Val{n}, max_feval::Int, max_minor::Int,
                           x::oneDeviceArray{Float64,1}, xl::oneDeviceArray{Float64,1},
                           xu::oneDeviceArray{Float64,1}, dA::oneDeviceArray{Float64,2},
                           dc::oneDeviceArray{Float64,1}) where {n}

        # Problem with dynamic memory
        g = oneLocalArray(Float64, n)
        xc = oneLocalArray(Float64, n)
        s = oneLocalArray(Float64, n)
        wa = oneLocalArray(Float64, n)
        wa1 = oneLocalArray(Float64, n)
        wa2 = oneLocalArray(Float64, n)
        wa3 = oneLocalArray(Float64, n)
        wa4 = oneLocalArray(Float64, n)
        wa5 = oneLocalArray(Float64, n)
        gfree = oneLocalArray(Float64, n)
        dsave = oneLocalArray(Float64, n)
        indfree = oneLocalArray(Int, n)
        iwa = oneLocalArray(Int, 2*n)
        isave = oneLocalArray(Int, n)

        A = oneLocalArray(Float64, (n,n))
        B = oneLocalArray(Float64, (n,n))
        L = oneLocalArray(Float64, (n,n))

        task = 0
        status = 0

        delta = 0.0
        fatol = 0.0
        frtol = 1e-12
        fmin = -1e32
        gtol = 1e-6
        cgtol = 0.1
        cg_itermax = n

        f = 0.0
        nfev = 0
        ngev = 0
        nhev = 0
        minor_iter = 0
        search = true

        while search

            # [0|1]: Evaluate function.

            if task == 0 || task == 1
                f = eval_f(n, x, dA, dc)
                nfev += 1
                if nfev >= max_feval
                    search = false
                end
            end

            # [2] G or H: Evaluate gradient and Hessian.

            if task == 0 || task == 2
                eval_g(n, x, g, dA, dc)
                eval_h(n, 1.0, x, A, dA)
                ngev += 1
                nhev += 1
                minor_iter += 1
            end

            # Initialize the trust region bound.

            if task == 0
                gnorm0 = ExaTronKernels.dnrm2(n, g, 1)
                delta = gnorm0
            end

            # Call Tron.

            if search
                delta, task = ExaTronKernels.dtron(n, x, xl, xu, f, g, A, frtol, fatol, fmin, cgtol,
                                            cg_itermax, delta, task, B, L, xc, s, indfree, gfree,
                                            isave, dsave, wa, iwa, wa1, wa2, wa3, wa4, wa5)
            end

            # [3] NEWX: a new point was computed.

            if task == 3
                gnorm_inf = ExaTronKernels.dgpnorm(n, x, xl, xu, g)
                if gnorm_inf <= gtol
                    task = 4
                end

                if minor_iter >= max_minor
                    status = 1
                    search = false
                end
            end

            # [4] CONV: convergence was achieved.

            if task == 4
                search = false
            end
        end

        return status, minor_iter
    end

    function driver_kernel_test(N::Val{n}, max_feval::Int, max_minor::Int,
                                dx::oneDeviceArray{Float64,1},
                                dxl::oneDeviceArray{Float64,1},
                                dxu::oneDeviceArray{Float64,1},
                                dA::oneDeviceArray{Float64,2},
                                dc::oneDeviceArray{Float64,1},
                                d_out::oneDeviceArray{Float64,1}
                                ) where {n}
        tx = get_local_id()
        bx = get_group_id()

        x = oneLocalArray(Float64, n)
        xl = oneLocalArray(Float64, n)
        xu = oneLocalArray(Float64, n)

        x[tx] = dx[tx]
        xl[tx] = dxl[tx]
        xu[tx] = dxu[tx]
        barrier()

        status, minor_iter = driver_kernel(N, max_feval, max_minor, x, xl, xu, dA, dc)

        if bx == 1
            d_out[tx] = x[tx]
        end
        barrier()

    end

    max_feval = 500
    max_minor = 100

    tron_A = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
    tron_B = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)
    tron_L = ExaTronKernels.TronDenseMatrix{Array{Float64,2}}(n)

    dx = oneArray{Float64}(undef, n)
    dxl = oneArray{Float64}(undef, n)
    dxu = oneArray{Float64}(undef, n)
    dA = oneArray{Float64,2}(undef, (n,n))
    dc = oneArray{Float64}(undef, n)
    d_out = oneArray{Float64}(undef, n)

    for i=1:itermax
        L = tril(rand(n,n))
        A = L*transpose(L)
        A .= tril(A) .+ (transpose(tril(A)) .- Diagonal(A))
        x = rand(n)
        xl = x .- abs.(rand(n))
        xu = x .+ abs.(rand(n))
        c = rand(n)

        tron_A.vals .= A

        copyto!(dx, x)
        copyto!(dxl, xl)
        copyto!(dxu, xu)
        copyto!(dA, tron_A.vals)
        copyto!(dc, c)

        function eval_f_cb(x)
            f = 0.5*(transpose(x)*A*x) + transpose(c)*x
            return f
        end

        function eval_g_cb(x, g)
            g .= A*x .+ c
        end

        function eval_h_cb(x, mode, rows, cols, scale, lambda, values)
            nz = 1
            if mode == :Structure
                @inbounds for j=1:n
                    @inbounds for i=j:n
                        rows[nz] = i
                        cols[nz] = j
                        nz += 1
                    end
                end
            else
                @inbounds for j=1:n
                    @inbounds for i=j:n
                        values[nz] = A[i,j]
                        nz += 1
                    end
                end
            end
        end

        nele_hess = div(n*(n+1), 2)
        tron = ExaTronKernels.createProblem(n, xl, xu, nele_hess, eval_f_cb, eval_g_cb, eval_h_cb; :matrix_type=>:Dense, :max_minor=>max_minor)
        copyto!(tron.x, x)
        status = ExaTronKernels.solveProblem(tron)
        @oneapi items=(n, nblk) groups=nblk driver_kernel_test(Val(n),max_feval,max_minor,dx,dxl,dxu,dA,dc,d_out)
        h_x = zeros(n)
        copyto!(h_x, d_out)

        @test norm(h_x .- tron.x) <= 1e-10
    end
end


text:
push!(ARGS, "../../input_files/dynamic/basin_switching/16GPa/switch1.dat")
include("../../Basin.jl")
pop!(ARGS)
push!(ARGS, "../../input_files/dynamic/basin_switching/16GPa/switch2.dat")
include("../../Basin.jl")
pop!(ARGS)
push!(ARGS, "../../input_files/dynamic/basin_switching/16GPa/switch3.dat")
include("../../Basin.jl")
pop!(ARGS)
push!(ARGS, "../../input_files/dynamic/basin_switching/16GPa/switch4.dat")
include("../../Basin.jl")
pop!(ARGS)
push!(ARGS, "../../input_files/dynamic/basin_switching/16GPa/switch5.dat")
include("../../Basin.jl")
pop!(ARGS)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module PolicyAssignmentsClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_PolicyAssignmentsApi.jl")

# export models
export convert, CloudError
export convert, ErrorAdditionalInfo
export convert, ErrorResponse
export convert, Identity
export convert, ParameterValues
export convert, ParameterValuesValue
export convert, PolicyAssignment
export convert, PolicyAssignmentListResult
export convert, PolicyAssignmentProperties
export convert, PolicySku

# export operations
export convert, PolicyAssignmentsApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
export proxTVCondat!

function proxTVCondat!(x::Vector{T}, λ::Float64; shape=[], kargs...) where T
  x_ = reshape(x, shape...)
  nhood = Int64[1 0 0;0 1 0; 0 0 1]
  omega = T[1,1,1]

  y = copy(x_)
  y[:] .= 0.0

  for d=1:length(omega)
    y_ = copy(x_)
    tv_denoise_3d_condat!(y_, nhood[d,:], λ*omega[d])
    y .+= y_ ./ length(omega)
  end
  x[:] = y[:]
  return y
end

mutable struct StartRange <: AbstractArray{CartesianIndices,3}
  x::CartesianIndices
  y::CartesianIndices
  z::CartesianIndices
end

Base.size(R::StartRange) = (Int(3),)
Base.IndexStyle(::Type{StartRange}) = IndexLinear()
Base.getindex(R::StartRange,i::Int) = if i==1 R.x elseif i==2 R.y elseif i==3 R.z end;

"""
This function checks if the cartesian index exceeds size.
"""
function inrange(size::Tuple{Int64,Int64,Int64},range::CartesianIndex{3})

    if range.I[1] > size[1] || range.I[1] < 1 || range.I[2] > size[2] ||
       range.I[2] < 1 || range.I[3] > size[3] || range.I[3] < 1
      return false
    else
      return true
    end

end

"""
This function returns a StartRange variable, which contains the start planes for
the 1d tv extraction.
"""
function get_startrange(size::Tuple{Int64,Int64,Int64},step::Array{T,1}) where {T<:Real}

  output = StartRange(
    CartesianIndices(CartesianIndex((0,0,0))),
    CartesianIndices(CartesianIndex((0,0,0))),
    CartesianIndices(CartesianIndex((0,0,0)))
  )

  output.x = CartesianIndices((1:step[1],1:size[2],1:size[3]))
  output.y = CartesianIndices((1:size[1],1:step[2],1:size[3]))
  output.z = CartesianIndices((1:size[1],1:size[2],1:step[3]))

  if step[1] < 0
    output.x = CartesianIndices(((size[1]+step[1]):size[1],1:size[2],1:size[3]))
  end

  if step[2] < 0
    output.y = CartesianIndices((1:size[1],(size[2]+step[2]):size[2],1:size[3]))
  end

  if step[3] < 0
    output.z = CartesianIndices((1:size[1],1:size[2],(size[3]+step[3]):size[3]))
  end

  return output
end

"""
This function returns the one dimensional tv problem depending on a start pixel
neighbor and a direction increment.
"""
function tv_get_onedim_data!(tvData::Array{T,3}, tvOneDim::Array{T,1}, neighbor,
    increment, arrayCount::Array{Int64,1}) where {T<:Real}

  tvSize = size(tvData)

  arrayCount[1] = 0

  # While neighbor does not exceeds tvSize, add data to 1d problem
  while true
    if inrange(tvSize,neighbor)
      arrayCount[1] = arrayCount[1]+1
      @inbounds tvOneDim[arrayCount[1]] = tvData[neighbor]
      neighbor = neighbor + increment
    else
      break
    end
  end

end

"""
This function sorts the 1d tv result back into the 3d data.
"""
function tv_push_onedim_data!(tvData::Array{T,3},tvOneDim::Array{T,1},
     arrayCount::Int64,neighbor,increment) where {T<:Real}

  for i=1:arrayCount
    @inbounds tvData[neighbor] = tvOneDim[i]
    neighbor = neighbor + increment
  end

end

"""
This function extracts 1d problems from the 3d data and starts the 1d tv function.
"""
function tv_denoise_3d_condat!(tvData::Array{T,3}, nhood::Array{Int64,1}, lambda) where {T<:Real}

  tvSize = size(tvData)
  cartRange = get_startrange(tvSize,nhood[:])
  increment = CartesianIndex((nhood[1],nhood[2],nhood[3]));
  tvOneDim = Array{eltype(tvData)}(undef, Int64(ceil(sqrt(tvSize[1]*tvSize[2]*tvSize[3]))))
  arrayCount = Array{Int64}(undef, 1)
  for R in cartRange

    for k in R
      neighbor = k;
      tv_get_onedim_data!(tvData,tvOneDim,neighbor,increment,arrayCount)

      tv_denoise_1d_condat!(tvOneDim,arrayCount[1],lambda)

      neighbor = k

      tv_push_onedim_data!(tvData,tvOneDim,arrayCount[1],neighbor,increment)
    end
  end
end

"""
This function performs the 1d tv algorithm.
"""
function tv_denoise_1d_condat!(c::Array{T,1},width::Int64,lambda) where {T<:Real}

  cLength = width

  k = 1
  k0 = 1
  umin = lambda
  umax = -lambda
  vmin = c[1] - lambda
  vmax = c[1] + lambda
  kplus = 1
  kminus = 1

  twolambda = 2*lambda
  minlambda = -lambda

  while true

    while k == cLength
      if umin < 0
        while true
          c[k0] = vmin
          k0 = k0+1
          !(k0<=kminus) && break
        end

        k=k0
        kminus=k
        vmin=c[kminus]
        umin=lambda

        umax = vmin + umin - vmax
      elseif umax > 0
        while true
          c[k0] = vmax
          k0 = k0+1
          !(k0<=kplus) && break
        end

        k=k0
        kplus=k
        vmax=c[kplus]
        umax=minlambda

        umin = vmax + umax - vmin
      else
          vmin = vmin + umin/(k-k0+1)
          while true
            c[k0] = vmin
            k0 = k0 + 1
            !(k0<=k) && break
          end
          return
      end
    end

    umin = umin + c[k+1] - vmin
    if umin < minlambda
      # Inplace soft thresholding
      #vmin = vmin>mu ? vmin-mu : vmin<-mu ? vmin+mu : 0.0
      while true
        c[k0] = vmin
        k0 = k0 + 1
        !(k0<=kminus) && break
      end
      k=k0
      kminus=k
      kplus=kminus
      vmin = c[kplus]
      vmax = vmin + twolambda
      umin = lambda
      umax = minlambda
    else
      umax = umax + c[k+1] - vmax
      if umax > lambda
        # Inplace soft thresholding
        #vmax = vmax>mu ? vmax-mu : vmax<-mu ? vmax+mu : 0.0;
        while true
          c[k0]=vmax
          k0 = k0 + 1
          !(k0<=kplus) && break
        end
        k = k0
        kminus = k
        kplus = kminus
        vmax = c[kplus]
        vmin = vmax - twolambda
        umin = lambda
        umax = minlambda
      else
        k = k + 1
        if umin >= lambda
          kminus = k
          vmin = vmin + (umin-lambda)/(kminus-k0+1)
          umin = lambda
        end
        if umax <= minlambda
          kplus = k
          vmax = vmax + (umax+lambda)/(kplus -k0 +1)
          umax = minlambda
        end
      end
    end
  end
end


text:
"""
Buffer one can easily insert/remove line segments from.
This is great for working with line data that is changing a lot with minimal
performanc impact. It's also great for batching draw calls, since you can
insert lines in different points of the program, and then draw it with a single
opengl draw call.
"""
struct LinesegmentBuffer{N}
    positions::GPUVector{Point{N, Float32}}
    colors::GPUVector{RGBAf0}
    thickness::GPUVector{Float32}
    robj::RenderObject
    range::Signal{Int}
end

function LinesegmentBuffer(pos::Point{N, <: AbstractFloat} = Point3f0(0)) where N
    positions = gpuvec(Point{N, Float32}[])
    colors = gpuvec(RGBAf0[])
    thickness = gpuvec(Float32[])
    range = Signal(0)
    robj = visualize(
        positions.buffer, :linesegment,
        color = colors.buffer,
        thickness = thickness.buffer,
        indices = range,
        boundingbox = Signal(AABB{Float32}())
    ).children[]
    robj.boundingbox = Signal(AABB{Float32}())
    LinesegmentBuffer{N}(
        positions,
        colors,
        thickness,
        robj,
        range
    )
end

same_length_array(array, value) = fill(value, length(array))
same_length_array(arr, value::Vector) = value

function Base.append!(lsb::LinesegmentBuffer{N}, pos::Vector{Point{N, Float32}}, color, thickness) where N
    thickv = Float32.(same_length_array(pos, thickness))
    append!(lsb.positions, pos)
    append!(lsb.colors, same_length_array(pos, to_color(color)))
    append!(lsb.thickness, thickv)
    bb = value(lsb.robj.boundingbox)
    for (s, pos) in zip(thickv, pos)
        pos3d = Vec{3, Float32}(to_nd(pos, Val{3}, 0))
        bb = update(bb, pos3d)
        bb = update(bb, pos3d .+ s)
    end
    push!(lsb.robj.boundingbox, bb)
    push!(lsb.range, length(lsb.positions))
    return
end

function Base.empty!(lsb::LinesegmentBuffer)
    resize!(lsb.positions, 0)
    resize!(lsb.colors, 0)
    resize!(lsb.thickness, 0)
    push!(lsb.range, 0)
    push!(lsb.robj.boundingbox, AABB{Float32}())
    return
end

function extract_view(x::SubArray)
    x.parent, x.indexes[1]
end
function extract_view(x::ArrayNode)
    p = if isa(to_value(x).parent, ArrayNode)
        to_value(x).parent
    else
        lift_node(x-> x.parent, x)
    end
    idx = if isa(to_value(x).parent, ArrayNode)
        to_value(x).indexes[1]
    else
        lift_node(x-> x.indexes[1], x)
    end
    p, idx
end

function lines_2glvisualize(kw_args)
    result = Dict{Symbol, Any}()
    for (k, v) in kw_args
        k in (:x, :y, :z, :scale, :rotation, :offset, :camera) && continue
        if k == :linestyle
            # TODO implement pattern as signal
            result[:pattern] = to_value(v)
            continue
        end
        if k == :drawover
            if v == true
                result[:prerender] = ()-> glDisable(GL_DEPTH_TEST)
            end
            continue
        end
        if k == :colornorm
            k = :color_norm
        end
        if k == :colormap
            k = :color_map
        end
        if k == :linewidth
            k = :thickness
        end
        if k == :positions
            k = :vertex
            if isa(to_value(v), SubArray)
                v, idx = extract_view(v)
                result[:indices] = to_signal(to_index_buffer((), idx))
            end
        end
        result[k] = to_signal(v)
    end
    result[:fxaa] = false
    result
end

function _lines(scene, style, attributes)
    attributes = lines_defaults(scene, attributes)
    data = lines_2glvisualize(attributes)
    pos = data[:vertex]
    delete!(data, :vertex)
    viz = GLVisualize._default(pos, Style(style), data)
    viz = GLVisualize.assemble_shader(viz).children[]
    insert_scene!(scene, style, viz, attributes)
end

for arg in ((:x, :y), (:x, :y, :z), (:positions,))
    insert_expr = map(arg) do elem
        :(attributes[$(QuoteNode(elem))] = $elem)
    end
    @eval begin
        function lines(scene::makie, $(arg...), attributes::Dict)
            $(insert_expr...)
            _lines(scene, :lines, attributes)
        end
        function linesegment(scene::makie, $(arg...), attributes::Dict)
            $(insert_expr...)
            _lines(scene, :linesegment, attributes)
        end
    end
end

function linesegment(scene::makie, pos::AbstractVector{<: Union{Tuple{P, P}, Pair{P, P}}}, attributes::Dict) where P <: Point
    positions = lift_node(x->reinterpret(P, x), to_node(pos))
    linesegment(scene, positions, attributes)
end

function arc(pmin, pmax, a1, a2)

    xy = Vector{Point2f0}(361)

    xcenter = (x_lin(xmin) + x_lin(xmax)) / 2.0;
    ycenter = (y_lin(ymin) + y_lin(ymax)) / 2.0;
    width = abs(x_lin(xmax) - x_lin(xmin)) / 2.0;
    height = abs(y_lin(ymax) - y_lin(ymin)) / 2.0;

    start = min(a1, a2);
    stop = max(a1, a2);
    start += (stop - start) / 360 * 360;

    n = 0;
    for a in start:stop
        x[n] = x_log(xcenter + width  * cos(a * M_PI / 180));
        y[n] = y_log(ycenter + height * sin(a * M_PI / 180));
        n += 1
    end
    if (n > 1)
        lines(x, y)
    end

end


text:
using Printf

println("Classification Tests:")
tests = [1:12, 28, 496, 220, 1184, 12496, 1264460, 790, 909, 562, 1064, 1488]
for i in tests
    (class, a) = aliquotclassifier(i)
    println(@sprintf("%8d => ", i), @sprintf("%16s, ", class), a)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ d4ecf168-0fb9-11ec-374d-af81d65de75e
md"""# Welcome to Week 2

## Learning Objectives

This week we will be laying the conceptual framework that will form the infrastructure for understanding and working on the Learning Health System cycle. What is this cycle?

- Data to Knowledge
- Knowledge to Practice
- Practice to Data

Our cycle has three concepts:

1. Data
2. Knowledge
3. Practice

In this week we will focus on three concepts, but we are going to ignore _Practice_ and we are going to insert the concept of _Information_.

1. Data
1. Information
1. Knowledge

(DIK)

All elements of our DIK trio our commonplace in our lives. In fact, I think it would be hard to go a day outside of an off-the-grid vacation where we didn't read or hear about one of the three. Nonetheless, we hope to introduce you to a new way of thinking about Data, Information, and Knowledge, a new way that will help you be more insightful and productive in our current era of healthcare.

The key to success for Week 2 is to slow down and think twice about the topics we are discussing, even if they seem commonplace.

"""

# ╔═╡ Cell order:
# ╠═d4ecf168-0fb9-11ec-374d-af81d65de75e


text:
#module CustomStencil
include("types.jl")
include("parsing.jl")
include("stencil_creation.jl")
include("kernel_gen.jl")
include("mgpu.jl")
include("fft_2d.jl")
include("fft_3d.jl")

using TimerOutputs
using StaticArrays
timers = Vector{TimerOutput}(undef,0)
using DSP

export @def_stencil_expression
export CreateStencilDefinition
export NewStencilInstance
export StencilInstance
export StencilDefinition
export ApplyStencil
export ApplyMultiGPU
@info collect(devices())
@info CUDA.version()
"""
    macro def_stencil_expression(expr::Expr)

Returns an expression with @sum macro and :abs :min: max functions
evaluated and replaced.
"""
macro def_stencil_expression(expr::Expr)#, data, current, t)
    expr = Meta.quot(expr)
    expr = eval_macro(expr, Symbol("@sum"))
    for i in [:abs, :max, :min]
        expr = eval_func!(expr, i)
    end
    return expr
end

function def_stencil_expression(expr::Expr)
    expr = eval_macro(expr, Symbol("@sum"))
    for i in [:abs, :max, :min]
        expr = eval_func!(expr, i)
    end
    return expr
end
```
    Used by eval to parse an expression to symbolic math
```
D = Array{Basic,1}(undef,1)
c = Array{Float64,1}(undef,1)
v = symbols("vsq_sym")

"""
    CreateStencilDefinition(expr::Expr, coefs::Union{Array{Float64,1}, Nothing})

Creates a Stencil definition struct given an expression and optinally a
coefficient matrix.
"""
function CreateStencilDefinition(expr::Expr, coefs::Union{Vector, Nothing})
    #println(expr)
    if coefs != nothing
        coefs = Float64.(coefs)
    end
    @show mn,mx = minmax_const(expr)
    m = mx > abs(mn) ? mx : abs(mn)
    if !(coefs isa Nothing)
        if (length(coefs)) == m + 1
            nothing
        elseif (length(coefs)) == 1000 + 2m + 1
            nothing
        else
            error("Coefficient matrix should have length of max radius + 1. is " ,length(coefs), " should be ", m +1)
        end
    end
    return StencilDefinition(expr, coefs, m, false)
end

"""
    CreateStencilDefinition(coefs::Array{Float64, 3}; uses_vsq=false)

Creates a Stencil definition struct given a coefficient array and whether
it should be multiplied by the vsq array (defaults to false.)
"""
function CreateStencilDefinition(coefs::Array{T, 3}; uses_vsq=false) where T <: Number
    # check coefs is cube

    coefs = Float64.(coefs)

    sz = size(coefs)
    !(sz[1] == sz[2] == sz[3]) &&  error("coefs array must have same dimensions size")
    max_radius = Int((sz[1] - 1)/2)

    return StencilDefinition(nothing, coefs, max_radius, uses_vsq)
end

function CreateStencilDefinition(coefs::Array{T, 2}; uses_vsq=false) where T <: Number
    # check coefs is cube

    coefs = Float64.(coefs)

    sz = size(coefs)
    !(sz[1] == sz[2]) &&  error("coefs array must have same dimensions size")
    max_radius = Int((sz[1] - 1)/2)

    return StencilDefinition(nothing, coefs, max_radius, uses_vsq)
end

"""
    function NewStencilInstance(varstencil::StencilDefinition...)

"""
function NewStencilInstance(varstencil::StencilDefinition...;prev_time_coeff=0, m_step=false, bdim=32)
    st_sym,max_radius = CombineStencils(varstencil...)
    fex,uses_vsq,fz,bz = GenStencilKernel(st_sym,max_radius,prev_time_coeff, bdim)
    kernel = ConvertToFunction(fex, uses_vsq)
    m_kernel = nothing
    m_rad = false
    m_fz = false
    m_bz = false
    fex2 = nothing
    new_sym = nothing
    if !uses_vsq && m_step != false && prev_time_coeff==0
        MAX_RAD = 16
        if max_radius*m_step > MAX_RAD
            m_step = false
            @warn "Combining timesteps results in too big stencil. Time combination aborted"
        else
            r = max_radius
            new_sym = combine_more_time_steps(st_sym, r, m_step)
            # for i = 2:m_step
            #     new_sym = CombineTimeSteps(st_sym, r, new_sym, r*i)
            # end
            m_rad = max_radius*m_step
            #println(new_sym)
            fex2,uses_vsq2,m_fz,m_bz = GenStencilKernel(new_sym,m_rad,prev_time_coeff, bdim)
            m_kernel = ConvertToFunction(fex2, uses_vsq2)
        end
    end
    return StencilInstance(varstencil, st_sym, max_radius,fz, bz, bdim, kernel, fex, uses_vsq,
                        m_step, m_kernel, fex2,new_sym, m_rad, m_fz, m_bz)
end

function ApplyStencil(st_inst::StencilInstance, org_data, t_steps::Integer; vsq=nothing)
    # Fix padding for multi time step
    #global timers
    to = TimerOutput()
    radius = st_inst.max_radius
    pad_radius = st_inst.max_radius
    offx = 0
    offy = 0
    if st_inst.combined_time_step != false
        pad_radius = st_inst.m_max_radius
        offx = st_inst.m_max_radius
        offy = st_inst.m_max_radius
    end

    # -------------------------------
    #radius = st_inst.max_radius
    data = PadData(pad_radius, org_data)
    bdimx = st_inst.bdim
    bdimy = st_inst.bdim
    dx = size(data,1)
    dy = size(data,2)
    dev_data = CuArray(data)
    dev_out = CUDA.zeros(Float32, size(data))
    dev_vsq = nothing
    if st_inst.uses_vsq
        if vsq isa Nothing
            error("vsq array not provided")
        end
        dev_vsq = CuArray(vsq)
    end
    bx = Int(floor((dx - 2*pad_radius)/bdimx))
    by = Int(floor((dy - 2*pad_radius)/bdimy))
    i = 0
    #st_inst.combined_time_step = false
    #println(st_inst.combined_time_step)
    @time begin
    while i < t_steps
        if st_inst.combined_time_step != false
            if i + st_inst.combined_time_step <= t_steps
                args = (dev_data, dev_out, 0, 0, 0, 0)
                if st_inst.uses_vsq
                    args = (args..., dev_vsq)
                end
                #println("Running combined $i")
                @timeit to "combined time step $(i)" begin
                @cuda(blocks=(bx,by,1), threads=(bdimx,bdimy),
                            shmem=((bdimx+2*st_inst.m_max_radius)*
                                (bdimy+2*st_inst.m_max_radius))*sizeof(Float32),
                            st_inst.m_kernel(args...))
                end # Timer end
                i += st_inst.combined_time_step
                dev_data,dev_out = dev_out,dev_data
                continue
            end
        end
        args = (dev_data, dev_out, 0, 0, offx, offy)
        if st_inst.uses_vsq
            args = (args..., dev_vsq)
        end
        #println("t = $(i), (bx,by) = $((bx,by))")

        @cuda(blocks=(bx,by,1), threads=(bdimx,bdimy),
                    shmem=((bdimx+2*radius)*(bdimy+2*radius))*sizeof(Float32),
                    st_inst.kernel(args...))

        dev_data,dev_out = dev_out,dev_data
        i += 1

    end # Timer end

    CUDA.synchronize(CuDefaultStream())
    end
    return view(Array(dev_data), (pad_radius+1):(dx-pad_radius), (pad_radius+1):(dy-pad_radius),:)

end

function ApplyFFTstencil(st_inst::StencilInstance, org_data, t_steps::Integer, multi=1)
    @show z = size(org_data, 3)
    template_single = matrix_from_expression(st_inst.stencil_sym, st_inst.max_radius)
    template_multi = nothing
    if multi > 1
        template_multi = DSP.conv(template_single, template_single)
        for i = 2:multi
            template_multi = DSP.conv(template_multi, template_single)
        end
    end

    ## 2d STENCIL
    if z == 1
        template_single = template_single[:,:,st_inst.max_radius+1]
        if template_multi != nothing
            template_single = template_multi[:,:,multi*st_inst.max_radius+1]
        end
        org_data = org_data[:,:,1]
        step = 1
        if multi > 1
            step = multi
        end
        t_steps = t_steps÷step
        return fft_stencil_2d(org_data, template_single, t_steps)
    end
    if template_multi != nothing
        template_single = template_multi
    end
    step = 1
    if multi > 1
        step = multi
    end
    t_steps = t_steps÷step
    return fft_stencil_3d(org_data, template_single, t_steps)

end

# function ApplyStencil(st_inst::StencilInstance, data, t_steps::Integer; vsq=nothing)
#     radius = st_inst.max_radius
#     bdimx = 32
#     bdimy = 32
#     dx = size(data,1)
#     dy = size(data,2)
#     dev_data = CuArray(data)
#     dev_out = CUDA.zeros(Float32, size(data))
#     dev_vsq = nothing
#     if st_inst.uses_vsq
#         if vsq isa Nothing
#             error("Vsq array not provided")
#         end
#         dev_vsq = CuArray(vsq)
#     end
#     bx = Int(floor((dx - 2*radius)/bdimx))
#     by = Int(floor((dy - 2*radius)/bdimy))
#     for i = 1:t_steps
#         if i%2 == 1
#             if st_inst.uses_vsq
#                 @cuda(blocks=(bx,by,1), threads=(bdimx,bdimy),
#                             shmem=((bdimx+2*radius)*(bdimy+2*radius))*sizeof(Float32),
#                             st_inst.kernel(dev_data,dev_out,0,dev_vsq))
#             else
#                 @cuda(blocks=(bx,by,1), threads=(bdimx,bdimy),
#                             shmem=((bdimx+2*radius)*(bdimy+2*radius))*sizeof(Float32),
#                             st_inst.kernel(dev_data,dev_out,0))
#             end
#         else
#             if st_inst.uses_vsq
#                 @cuda(blocks=(bx,by,1), threads=(bdimx,bdimy),
#                             shmem=((bdimx+2*radius)*(bdimy+2*radius))*sizeof(Float32),
#                             st_inst.kernel(dev_out,dev_data,0,dev_vsq))
#             else
#                 @cuda(blocks=(bx,by,1), threads=(bdimx,bdimy),
#                             shmem=((bdimx+2*radius)*(bdimy+2*radius))*sizeof(Float32),
#                             st_inst.kernel(dev_out,dev_data,0))
#             end
#         end
#         println("t = $(i), (bx,by) = $((bx,by))")
#     end
#
#     (t_steps%2 == 1) ? (return Array(dev_out)) : (return Array(dev_data))
#
# end

#end


text:
@for_libpetsc begin
    function initialized(::Type{$PetscScalar})
        r_flag = Ref{PetscBool}()
        @chk ccall((:PetscInitialized, $libpetsc), PetscErrorCode, (Ptr{PetscBool},), r_flag)
        return r_flag[] == PETSC_TRUE
    end
    function initialize(::Type{$PetscScalar})
        if !initialized($PetscScalar)
            MPI.Initialized() || MPI.Init()
            @chk ccall((:PetscInitializeNoArguments, $libpetsc), PetscErrorCode, ())

            # disable signal handler
            @chk ccall((:PetscPopSignalHandler, $libpetsc), PetscErrorCode, ())

            atexit(() -> finalize($PetscScalar))
        end
        return nothing
    end
    function finalized(::Type{$PetscScalar})
        r_flag = Ref{PetscBool}()
        @chk ccall((:PetscFinalized, $libpetsc), PetscErrorCode, (Ptr{PetscBool},), r_flag)
        return r_flag[] == PETSC_TRUE
    end
    function finalize(::Type{$PetscScalar} )
        if !finalized($PetscScalar)
            @chk ccall((:PetscFinalize, $libpetsc), PetscErrorCode, ())
        end
        return nothing
    end
end

"""
    PETSc.initialize([T])

Initialize PETSc for scalar type `T`, if not already initialized. If no `T` is
provided, then it will be initialized for all supported scalar types.

Additionally:
 - This will initialize MPI if it has not already been initialized.
 - It will disable the PETSc signal handler (via
   [`PetscPopSignalHandler`](https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscPopSignalHandler.html))
 - Add an [`atexit`](https://docs.julialang.org/en/v1/base/base/#Base.atexit)
   hook to call [`PETSc.finalize`](@ref).
"""
function initialize()
    map(initialize, scalar_types)
    return nothing
end

"""
    PETSc.finalize([T])

Finalize PETSc for scalar type `T`. If no `T` is provided, then it will be finalized for all supported scalar types.

It is generally not necessary to call this function directly, as it is added as an [`atexit`](https://docs.julialang.org/en/v1/base/base/#Base.atexit) hook in [`PETSc.initialize`](@ref).
"""
function finalize()
    map(finalize, scalar_types)
    return nothing
end


text:
module OptionPricing

using Reexport
@reexport using Dates
import Statistics: mean

export AmericanOption, Put, simulate_price, price_option, SwingOption,
SwingContract, Valuation, PriceProcess

abstract type AbstractOption end
abstract type RealOption <:AbstractOption end
abstract type FinancialOption <: AbstractOption end

# Option Types may use Traits
abstract type OptionStyle end
struct Put <: OptionStyle end
struct Call <: OptionStyle end

include("american.jl")
include("swing.jl")

# power_series takes an an array and returns an array with each each
# value raised to 0 up to order.
function power_series(price::AbstractArray, order::Int)
    return price .^ (0:order)'
end

end # module


text:
"""
    find_feg_separator_lt(skeleton, faces, root)

Finds the a valid separator for a finite element graph given its skeleton and
faces of its planar embedding represented as a set of sets of edges. The search
will initiate on a bfs tree rooted at the given vertex.
"""
function find_feg_separator_lt(skeleton::LabeledGraph{T}, faces::Set{Set{Pair{T, T}}}, root::T)::Tuple{Set{T}, Set{T}, Set{T}} where {T}
    g = _find_finite_element_graph(skeleton, faces)

    g_star_star = copy(skeleton)
    new_vertices = Dict{Number, Set{Pair{T, T}}}()
    # Only need to add vertex to faces with more than 3 vertices
    faces = filter(f -> length(f) > 3, faces)
    for face in faces
        new_vertex = rand(T)
        push!(new_vertices, new_vertex => face)
        lg_add_vertex!(g_star_star, new_vertex)
        face_vertices = reduce((x, y) -> union!(x, Set([y.first, y.second])), face, init=Set{Number}())
        for vertex in face_vertices
            lg_add_edge!(g_star_star, vertex, new_vertex)
        end
    end

    (c_star_star, a_star_star, b_star_star) = find_separator_lt(g_star_star, root)

    a = filter(vertex -> !(vertex in keys(new_vertices)), a_star_star)
    b = filter(vertex -> !(vertex in keys(new_vertices)), b_star_star)

    for pair in new_vertices
        if (pair.first in c_star_star)
            bad_face_vertecies =  reduce((x, y) -> union!(x, Set([y.first, y.second])), pair.second, init=Set{Number}())

            a_star_star_boundry = filter(vertex -> vertex in a_star_star, bad_face_vertecies)
            b_star_star_boundry = filter(vertex -> vertex in b_star_star, bad_face_vertecies)

            if length(a_star_star_boundry) < length(b_star_star_boundry)
                setdiff!(a, a_star_star_boundry)
                setdiff!(a_star_star, a_star_star_boundry)
            else
                setdiff!(b, b_star_star_boundry)
                setdiff!(b_star_star, b_star_star_boundry)
            end
        end
    end

    return (setdiff(lg_vertices(g), a, b), a, b)
end

"""
    find_feg_separator_lt_best(skeleton, faces)

Finds the a valid separator for a finite element graph given its skeleton and
faces of its planar embedding represented as a set of sets of edges. Repeats
the search from every possible BFS root and returns the most balanced
separator.
"""
function find_feg_separator_lt_best(skeleton::LabeledGraph{T}, faces::Set{Set{Pair{T, T}}})::Tuple{Set{T}, Set{T}, Set{T}} where {T}
    separators = map(root -> find_feg_separator_lt(skeleton, faces, root), collect(keys(skeleton.labels)))

    balances = map(separator -> min(length(separator[2]) / length(separator[3]), length(separator[3]) / length(separator[2])), separators)
    ######## See _find_feg_separator_lt_no_empty in biclique.jl
    # balances = map(separator -> begin
    #     if !isempty(separator[2]) && !isempty(separator[3])
    #         return min(length(separator[2]) / length(separator[3]), length(separator[3]) / length(separator[2]))
    #     else
    #         return -1.0
    #     end
    # end, separators)

    max_balance = max(balances...)
    for i in 1:length(balances)
        if balances[i] == max_balance return separators[i] end
    end
    return separators[1]

end

"""
    find_separator_lt(lg, root)

Finds a valid separator for a planar graph using the Lipton-Tarjan algorithm.
The search will initiate on a bfs tree rooted at the given vertex.
"""
function find_separator_lt(lg::LabeledGraph{T}, root::T)::Tuple{Set{T}, Set{T}, Set{T}} where {T}
    levels = @pipe LightGraphs.bfs_tree(lg.graph, lg.labels[root]) |> _find_bfs_levels(_, lg.labels[root]) |> map(level -> Set(convert_vertices(lg.labels, collect(level))), _)

    # Phase I
    middle_index = 1
    for i in 2:length(levels)
        if sum(map(j -> length(levels[j]), 1:i)) > (LightGraphs.nv(lg.graph) / 2)
            middle_index = i
            break
        end
    end
    middle = levels[middle_index]

    if length(middle) < (2 * sqrt(2 * LightGraphs.nv(lg.graph)))
        above = middle_index == 1 ? Set() : @pipe map(i -> levels[i], 1:(middle_index - 1)) |> reduce(union, _)
        below = middle_index < length(levels) ? (@pipe map(i -> levels[i], (middle_index + 1):length(levels)) |> reduce(union, _)) : Set()
        a = length(above) < length(below) ? above : below;
        b = length(above) < length(below) ? below : above;
        return (middle, a, b)
    end

    # Phase II
    upper_index = -1;
    for i in middle_index:-1:1
        distance = middle_index - i
        if length(levels[i]) < 2 * (sqrt(LightGraphs.nv(lg.graph)) - distance)
            upper_index = i
            break
        end
    end
    if sum(map(i -> length(levels[i]), 1:upper_index)) > 1/3
        # upper is separator
        above = upper_index > 1 ? (@pipe map(i -> levels[i], 1:(upper_index - 1)) |> reduce(union, _)) : Set()
        below = upper_index < length(levels) ? (@pipe map(i -> levels[i], (upper_index + 1):length(levels)) |> reduce(union, _)) : Set()
        a = length(above) < length(below) ? above : below;
        b = length(above) < length(below) ? below : above;
        return (levels[upper_index], a, b)
    end

    lower_index = -1;
    for i in middle_index:length(levels)
        distance = i - middle_index
        if length(levels[i]) < 2 * (sqrt(LightGraphs.nv(lg.graph)) - distance)
            lower_index = i
            break
        end
    end
    if sum(map(i -> length(levels[i]), lower_index:length(levels))) > 1/3
        # lower is separator
        above = lower_index > 1 ? (@pipe map(i -> levels[i], 1:(lower_index - 1)) |> reduce(union, _)) : Set()
        below = lower_index < length(levels) ? (@pipe map(i -> levels[i], (lower_index + 1):length(levels)) |> reduce(union, _)) : Set()
        a = length(above) < length(below) ? above : below;
        b = length(above) < length(below) ? below : above;
        return (levels[lower_index], a, b)
    end

    # Phase III
    return find_separator_fcs(lg, root)

end

"""
    find_separator_fcs(lg, root)

Finds a valid separator for a planar graph using finite cycles as separators.
The search will initiate on a bfs tree rooted at the given vertex.
"""
function find_separator_fcs(lg::LabeledGraph{T}, root::T)::Tuple{Set{T}, Set{T}, Set{T}} where {T}
    parents = LightGraphs.bfs_parents(lg.graph, lg.labels[root])

    tree_edges = map(tup -> LightGraphs.Edge(tup[1] => tup[2]), zip(parents, 1:length(parents)))
    non_tree_edge = filter(edge -> !(in(edge, tree_edges) || in(reverse(edge), tree_edges)), collect(LightGraphs.edges(lg.graph)))[1]

    fundamental_cycle = @pipe _find_fundamental_cycle(parents, non_tree_edge) |> Set(convert_vertices(lg.labels, collect(_)))
    return (fundamental_cycle, _find_partitions(lg, fundamental_cycle)...)
end

"""
    find_separator_fcs_best(lg, root)

Finds a valid separator for a planar graph using finite cycles as separators.
Repeats the search from every possible BFS root and returns the most balanced
separator.
"""
function find_separator_fcs_best(lg::LabeledGraph{T}, root::T)::Tuple{Set{T}, Set{T}, Set{T}} where {T}
    # THIS ALGORITHM IS O(n^2) BECAUSE IT CHOOSES THE MOST BALANCED
    parents = LightGraphs.bfs_parents(lg.graph, lg.labels[root])

    tree_edges = map(tup -> LightGraphs.Edge(tup[1] => tup[2]), zip(parents, 1:length(parents)))
    non_tree_edges = filter(edge -> !(in(edge, tree_edges) || in(reverse(edge), tree_edges)), collect(LightGraphs.edges(lg.graph)))

    fundamental_cycles = @pipe map(non_tree_edge -> _find_fundamental_cycle(parents, non_tree_edge), non_tree_edges) |> map(cycle -> Set(convert_vertices(lg.labels, collect(cycle))), _)
    partitions = map(cycle -> (cycle, _find_partitions(lg, cycle)...), fundamental_cycles)
    # TODO: What happens when sets are empty? Can get Inf
    ######## See _find_feg_separator_lt_no_empty in biclique.jl
    balances = map(partition -> min(length(partition[2]) / length(partition[3]), length(partition[3]) / length(partition[2])), partitions)

    max_balance = max(balances...)
    for i in 1:length(balances)
        if balances[i] == max_balance return partitions[i] end
    end
    return partitions[1]
end

"""
    _find_partitions(lg, separator)

Given a graph and separator find the sets of vertices separated.
"""
function _find_partitions(lg::LabeledGraph{T}, separator::Set{T})::Tuple{Set{T}, Set{T}} where {T}
    lg_without_separator = copy(lg)

    for vertex in separator
        lg_rem_vertex!(lg_without_separator, vertex)
    end

    components = @pipe LightGraphs.connected_components(lg_without_separator.graph) |>
        map(component -> convert_vertices(lg_without_separator.labels, component), _) |>
        sort(_, by=length, rev=true)
    a = Array{T, 1}()
    b = Array{T, 1}()
    for component in components
        append!(sum(map(length, a)) < sum(map(length, b)) ? a : b, component)
    end

    return (Set(a), Set(b))
end

"""
    _find_fundamental_cycle(parents, non_tree_edge)

Given a bfs tree (in array form as provided by LightGraphs) and a non tree edge
return the fundamtal cycle as a set of vertices.
"""
function _find_fundamental_cycle(parents::Vector{Int}, non_tree_edge::LightGraphs.Edge)::Set{Int}
    left_vertex = LightGraphs.src(non_tree_edge)
    left_path = [left_vertex]
    while left_vertex != parents[left_vertex]
        left_vertex = parents[left_vertex]
        pushfirst!(left_path, left_vertex)
    end

    right_vertex = LightGraphs.dst(non_tree_edge)
    right_path = [right_vertex]
    while right_vertex != parents[right_vertex]
        right_vertex = parents[right_vertex]
        pushfirst!(right_path, right_vertex)
    end

    for i in 1:max(length(left_path), length(right_path))
        if left_path[i] != right_path[i]
            return Set(union(left_path[i:end], right_path[i:end], left_path[i - 1]))
        end
    end

    return Set(union(left_path, right_path))
end

"""
    _find_bfs_levels(tree, root)

Run breadth first search on a graph and return its BFS tree as an array of sets
representing levels of the tree.
"""
function _find_bfs_levels(tree::LightGraphs.AbstractGraph, root::Int)::Array{Set{Int}, 1}
    levels = Dict{Int, Set{Int}}()

    levels[1] = Set([root])
    i = 2
    while true
        levels[i] = Set(reduce(union, map(vertex -> LightGraphs.outneighbors(tree, vertex), collect(levels[i - 1]))))

        if isempty(levels[i]) break end

        i += 1
    end

    res = Array{Set{Int}}(undef, i - 1)
    for i in 1:length(res)
        res[i] = levels[i]
    end

    return res
end

"""
    _is_valid_separator(lg, separator, a, b)

Checks if a separator is valid on graph lg and separates A and B
"""
function _is_valid_separator(lg::LabeledGraph{T}, separator::Set{T}, a::Set{T}, b::Set{T})::Bool where T
    is_separating = true
    for source in a
        break_early = false
        for destination in b
            is_path = ClutteredEnvPathOpt.LightGraphs.has_path(lg.graph, source, destination, exclude_vertices=collect(separator))
            is_valid = is_separating && !is_path
            if !is_valid
                break_early = true
                break
            end
        end

        if break_early break end
    end

   is_covering =  all(map(vertex -> in(vertex, union(separator, a, b)), collect(keys(lg.labels))))

   is_disjoint =
    all(map(vertex -> !in(vertex, union(a, b)), collect(separator))) &&
    all(map(vertex -> !in(vertex, union(b, separator)), collect(a))) &&
    all(map(vertex -> !in(vertex, union(a, separator)), collect(b)))

    return is_separating && is_covering && is_disjoint
end

"""
    _find_finite_element_graph(skeleton, faces)

Creates a finite element graph from a skeleton and set of sets of edges
representing faces.
"""
function _find_finite_element_graph(skeleton::LabeledGraph{T}, faces::Set{Set{Pair{T, T}}})::LabeledGraph{T} where {T}
    g = copy(skeleton)

    for face in faces
        face_vertices = reduce((x, y) -> union!(x, Set([y.first, y.second])), face, init=Set{Number}())
        for edge in face
            for dest in face_vertices
                if (edge.first != dest)
                    lg_add_edge!(g, edge.first, dest)
                end
            end
        end
    end

    return g
end

"""
    pp_expell(lg, separator, a, b)

A postprocessing algorithm to shrink the size of a separator by expelling
vertices in the separator not connected to both A and B.
"""
function pp_expell(lg::LabeledGraph{T}, separator::Set{T}, a::Set{T}, b::Set{T})::Tuple{Set{T}, Set{T}, Set{T}} where {T}
    pp_separator = copy(separator)
    pp_a = copy(a)
    pp_b = copy(b)

    for separator_vertex in separator
        source = lg.labels[separator_vertex]

        is_connected_a = any(
            map(
                destination -> LightGraphs.has_path(lg.graph, source, destination, exclude_vertices=collect(convert_vertices_rev(lg.labels, setdiff(pp_separator, Set([separator_vertex]))))),
                collect(convert_vertices_rev(lg.labels, pp_a))
            )
        )

        is_connected_b = any(
            map(
                destination -> LightGraphs.has_path(lg.graph, source, destination, exclude_vertices=collect(convert_vertices_rev(lg.labels, setdiff(pp_separator, Set([separator_vertex]))))),
                collect(convert_vertices_rev(lg.labels, pp_b))
            )
        )

        if !is_connected_a && !is_connected_b
            push!(length(pp_a) < length(pp_b) ? pp_a : pp_b, separator_vertex)
            delete!(pp_separator, separator_vertex)
        elseif is_connected_a && !is_connected_b
            push!(pp_a, separator_vertex)
            delete!(pp_separator, separator_vertex)
        elseif !is_connected_a && is_connected_b
            push!(pp_b, separator_vertex)
            delete!(pp_separator, separator_vertex)
        end
    end

    # Assure A and B are nonempty

    # triggered = false
    # if isempty(pp_a) || isempty(pp_b)
    #     println("\n#########\nA: $pp_a\nB:$pp_b\nC:$pp_separator")
    #     triggered = true
    # end

    if isempty(pp_a)
        # Move vertices from B into A
        # First, need to move vertices in B connected to every other vertex in B into C
        for b in pp_b
            if length(pp_b) <= 1
                #println("Don't want an empty B")
                break
            end
            source = lg.labels[b]
            num_edges_with_pp_b = sum(
                map(destination -> begin
                    if LightGraphs.has_edge(lg.graph, source, destination)
                        return 1
                    else
                        return 0
                    end
                end,
                    collect(convert_vertices_rev(lg.labels, setdiff(pp_b, Set(b))))
                )
            )
            if num_edges_with_pp_b == length(setdiff(pp_b, b))
                #println("Moving $b from B into C")
                push!(pp_separator, b)
                delete!(pp_b, b)
            end
        end

        # Next we move into A either
        # 1) vertices in B that aren't adjacent to any of the other vertices in B, or
        # 2) a vertex in B that is adjacent to the least amount of other vertices in B, also moving
        # those vertices it is adjacent to into C
        size_pp_b = length(pp_b)
        added = 0
        nv = length(pp_separator) + length(pp_b) + length(pp_a)
        min = nv * (nv - 1) / 2
        min_vertex = -1

        # Case 1
        for b in pp_b
            if added <= floor(size_pp_b/2) && length(pp_b) > 1 # don't want to add more than half into A and don't want B to be empty
                source = lg.labels[b]
                num_edges_with_pp_b = sum(
                    map(destination -> begin
                        if LightGraphs.has_edge(lg.graph, source, destination)
                            return 1
                        else
                            return 0
                        end
                    end,
                        collect(convert_vertices_rev(lg.labels, setdiff(pp_b, Set(b))))
                    )
                )
                if num_edges_with_pp_b == 0
                    #println("Moving $b from B into A")
                    push!(pp_a, b)
                    delete!(pp_b, b)
                    added += 1
                elseif num_edges_with_pp_b < min
                    min = num_edges_with_pp_b
                    min_vertex = b
                end
            end
        end

        # Case 2
        if added == 0
            source = lg.labels[min_vertex]
            rev = _reverse_labels(lg.labels)
            connections = map(destination -> begin
                if LightGraphs.has_edge(lg.graph, source, destination)
                    return rev[destination]
                else
                    return 0
                end
            end,
                collect(convert_vertices_rev(lg.labels, setdiff(pp_b, Set(min_vertex))))
            )

            connections = filter(c -> c != 0, connections)

            for v in connections
                push!(pp_separator, v)
                delete!(pp_b, v)
            end

            push!(pp_a, min_vertex)
            delete!(pp_b, min_vertex)
        end

    elseif isempty(pp_b)
        # Move vertices from A into B
        # First, need to move vertices in A connected to every other vertex in A into C
        for a in pp_a
            if length(pp_a) <= 1
                #println("Don't want an empty A")
                break
            end
            source = lg.labels[a]
            num_edges_with_pp_a = sum(
                map(destination -> begin
                    if LightGraphs.has_edge(lg.graph, source, destination)
                        return 1
                    else
                        return 0
                    end
                end,
                    collect(convert_vertices_rev(lg.labels, setdiff(pp_a, Set(a))))
                )
            )
            if num_edges_with_pp_a == length(setdiff(pp_a, a))
                #println("Moving $a from A into C")
                push!(pp_separator, a)
                delete!(pp_a, a)
            end
        end

        # Next we move into B either
        # 1) vertices in A that aren't adjacent to any of the other vertices in A, or
        # 2) a vertex in A that is adjacent to the least amount of other vertices in A, also moving
        # those vertices it is adjacent to into C
        size_pp_a = length(pp_a)
        added = 0
        nv = length(pp_separator) + length(pp_b) + length(pp_a)
        min = nv * (nv - 1) / 2
        min_vertex = -1

        # Case 1
        for a in pp_a
            if added <= floor(size_pp_a/2) && length(pp_a) > 1 # don't want to add more than half into B and don't want A to be empty
                source = lg.labels[a]
                num_edges_with_pp_a = sum(
                    map(destination -> begin
                        if LightGraphs.has_edge(lg.graph, source, destination)
                            return 1
                        else
                            return 0
                        end
                    end,
                        collect(convert_vertices_rev(lg.labels, setdiff(pp_a, Set(a))))
                    )
                )
                if num_edges_with_pp_a == 0
                    #println("Moving $a from A into B")
                    push!(pp_b, a)
                    delete!(pp_a, a)
                    added += 1
                elseif num_edges_with_pp_a < min
                    min = num_edges_with_pp_a
                    min_vertex = a
                end
            end
        end

        # Case 2
        if added == 0
            source = lg.labels[min_vertex]
            rev = _reverse_labels(lg.labels)
            connections = map(destination -> begin
                if LightGraphs.has_edge(lg.graph, source, destination)
                    return rev[destination]
                else
                    return 0
                end
            end,
                collect(convert_vertices_rev(lg.labels, setdiff(pp_a, Set(min_vertex))))
            )

            connections = filter(c -> c != 0, connections)

            for v in connections
                push!(pp_separator, v)
                delete!(pp_a, v)
            end

            push!(pp_b, min_vertex)
            delete!(pp_a, min_vertex)
        end

    end

    # if triggered
    #     println("\nThe new sets are \nA: $pp_a\nB:$pp_b\nC:$pp_separator\n")
    # end

    return (pp_separator, pp_a, pp_b)
end