text:
using CUDA
using ExaTronKernels
using LinearAlgebra
using Random
using Test

try
    tmp = CuArray{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
    function dicf_test(n::Int, d_in::CuDeviceArray{Float64},
                        d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        L = CuDynamicSharedArray(Float64, (n,n))
        for i in 1:n
            L[i,tx] = d_in[i,tx]
        end
        CUDA.sync_threads()

        # Test Cholesky factorization.
        ExaTronKernels.dicf(n,L)

        if bx == 1
            for i in 1:n
                d_out[i,tx] = L[i,tx]
            end
        end
        CUDA.sync_threads()
        return
    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 = CuArray{Float64,2}(undef, (n,n))
        d_out = CuArray{Float64,2}(undef, (n,n))
        copyto!(d_in, tron_A.vals)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=(n^2*sizeof(Float64)) dicf_test(n,d_in,d_out)
        h_L = zeros(n,n)
        copyto!(h_L, d_out)

        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-10
    end
end

@testset "dicfs" begin
    function dicfs_test(n::Int, alpha::Float64,
                        dA::CuDeviceArray{Float64},
                        d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        wa1 = CuDynamicSharedArray(Float64, n)
        wa2 = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        A = CuDynamicSharedArray(Float64, (n,n), (2*n)*sizeof(Float64))
        L = CuDynamicSharedArray(Float64, (n,n), (2*n+n^2)*sizeof(Float64))

        @inbounds for j=1:n
            A[j,tx] = dA[j,tx]
        end
        CUDA.sync_threads()

        ExaTronKernels.dicfs(n, alpha, A, L, wa1, wa2)
        if bx == 1
            @inbounds for j=1:n
                d_out[j,tx] = L[j,tx]
            end
        end
        CUDA.sync_threads()

        return
    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 = CuArray{Float64,2}(undef, (n,n))
        d_out = CuArray{Float64,2}(undef, (n,n))
        alpha = 1.0
        copyto!(dA, tron_A.vals)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((2*n+2*n^2)*sizeof(Float64)) dicfs_test(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)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((2*n+2*n^2)*sizeof(Float64)) dicfs_test(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
    function dcauchy_test(n::Int,dx::CuDeviceArray{Float64},
                            dl::CuDeviceArray{Float64},
                            du::CuDeviceArray{Float64},
                            dA::CuDeviceArray{Float64},
                            dg::CuDeviceArray{Float64},
                            delta::Float64,
                            alpha::Float64,
                            d_out1::CuDeviceArray{Float64},
                            d_out2::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        xl = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        xu = CuDynamicSharedArray(Float64, n, (2*n)*sizeof(Float64))
        g =  CuDynamicSharedArray(Float64, n, (3*n)*sizeof(Float64))
        s =  CuDynamicSharedArray(Float64, n, (4*n)*sizeof(Float64))
        wa = CuDynamicSharedArray(Float64, n, (5*n)*sizeof(Float64))
        A =  CuDynamicSharedArray(Float64, (n,n), (6*n)*sizeof(Float64))

        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
        CUDA.sync_threads()

        return
    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 = CuArray{Float64}(undef, n)
        dl = CuArray{Float64}(undef, n)
        du = CuArray{Float64}(undef, n)
        dg = CuArray{Float64}(undef, n)
        dA = CuArray{Float64,2}(undef, (n,n))
        d_out1 = CuArray{Float64}(undef, n)
        d_out2 = CuArray{Float64}(undef, n)
        copyto!(dx, x)
        copyto!(dl, xl)
        copyto!(du, xu)
        copyto!(dg, g)
        copyto!(dA, A.vals)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((6*n+n^2)*sizeof(Float64)) dcauchy_test(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
    function dtrpcg_test(n::Int, delta::Float64, tol::Float64,
                            stol::Float64, d_in::CuDeviceArray{Float64},
                            d_g::CuDeviceArray{Float64},
                            d_out_L::CuDeviceArray{Float64},
                            d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        A = CuDynamicSharedArray(Float64, (n,n))
        L = CuDynamicSharedArray(Float64, (n,n), (n^2)*sizeof(Float64))

        g = CuDynamicSharedArray(Float64, n, (2*n^2)*sizeof(Float64))
        w = CuDynamicSharedArray(Float64, n, (2*n^2 + n)*sizeof(Float64))
        p = CuDynamicSharedArray(Float64, n, (2*n^2 + 2*n)*sizeof(Float64))
        q = CuDynamicSharedArray(Float64, n, (2*n^2 + 3*n)*sizeof(Float64))
        r = CuDynamicSharedArray(Float64, n, (2*n^2 + 4*n)*sizeof(Float64))
        t = CuDynamicSharedArray(Float64, n, (2*n^2 + 5*n)*sizeof(Float64))
        z = CuDynamicSharedArray(Float64, n, (2*n^2 + 6*n)*sizeof(Float64))

        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]
        CUDA.sync_threads()

        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
        CUDA.sync_threads()

        return
    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 = CuArray{Float64,2}(undef, (n,n))
        d_g = CuArray{Float64}(undef, n)
        d_out_L = CuArray{Float64,2}(undef, (n,n))
        d_out = CuArray{Float64}(undef, n)
        copyto!(d_in, A)
        copyto!(d_g, g)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((2*n^2+7*n)*sizeof(Float64)) dtrpcg_test(n,delta,tol,stol,d_in,d_g,d_out_L,d_out)
        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
    function dprsrch_test(n::Int,d_x::CuDeviceArray{Float64},
                            d_xl::CuDeviceArray{Float64},
                            d_xu::CuDeviceArray{Float64},
                            d_g::CuDeviceArray{Float64},
                            d_w::CuDeviceArray{Float64},
                            d_A::CuDeviceArray{Float64},
                            d_out1::CuDeviceArray{Float64},
                            d_out2::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        xl = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        xu = CuDynamicSharedArray(Float64, n, (2*n)*sizeof(Float64))
        g = CuDynamicSharedArray(Float64, n, (3*n)*sizeof(Float64))
        w = CuDynamicSharedArray(Float64, n, (4*n)*sizeof(Float64))
        wa1 = CuDynamicSharedArray(Float64, n, (5*n)*sizeof(Float64))
        wa2 = CuDynamicSharedArray(Float64, n, (6*n)*sizeof(Float64))
        A = CuDynamicSharedArray(Float64, (n,n), (7*n)*sizeof(Float64))
        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]
        CUDA.sync_threads()

        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
        CUDA.sync_threads()

        return
    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 = CuArray{Float64}(undef, n)
        dl = CuArray{Float64}(undef, n)
        du = CuArray{Float64}(undef, n)
        dg = CuArray{Float64}(undef, n)
        dw = CuArray{Float64}(undef, n)
        dA = CuArray{Float64,2}(undef, (n,n))
        d_out1 = CuArray{Float64}(undef, n)
        d_out2 = CuArray{Float64}(undef, n)
        copyto!(dx, x)
        copyto!(dl, xl)
        copyto!(du, xu)
        copyto!(dg, g)
        copyto!(dw, w)
        copyto!(dA, A.vals)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((7*n+n^2)*sizeof(Float64)) dprsrch_test(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
    function daxpy_test(n::Int, da, d_in::CuDeviceArray{Float64},
                        d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        y = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        x[tx] = d_in[tx]
        y[tx] = d_in[tx + n]
        CUDA.sync_threads()

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

        return
    end

    for i=1:itermax
        da = rand(1)[1]
        h_in = rand(2*n)
        h_out = zeros(n)
        d_in = CuArray{Float64}(undef, 2*n)
        d_out = CuArray{Float64}(undef, n)
        copyto!(d_in, h_in)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((2*n)*sizeof(Float64)) daxpy_test(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
    function dssyax_test(n::Int,d_z::CuDeviceArray{Float64},
                            d_in::CuDeviceArray{Float64},
                            d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        z = CuDynamicSharedArray(Float64, n)
        q = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        A = CuDynamicSharedArray(Float64, (n,n), (2*n)*sizeof(Float64))
        for i in 1:n
            A[i,tx] = d_in[i,tx]
        end
        z[tx] = d_z[tx]
        CUDA.sync_threads()

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

        return
    end

    for i=1:itermax
        z = rand(n)
        h_in = rand(n,n)
        h_out = zeros(n)
        d_z = CuArray{Float64}(undef, n)
        d_in = CuArray{Float64,2}(undef, (n,n))
        d_out = CuArray{Float64}(undef, n)
        copyto!(d_z, z)
        copyto!(d_in, h_in)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((2*n+n^2)*sizeof(Float64)) dssyax_test(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
    function dmid_test(n::Int, dx::CuDeviceArray{Float64},
                        dl::CuDeviceArray{Float64},
                        du::CuDeviceArray{Float64},
                        d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        xl = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        xu = CuDynamicSharedArray(Float64, n, (2*n)*sizeof(Float64))
        x[tx] = dx[tx]
        xl[tx] = dl[tx]
        xu[tx] = du[tx]
        CUDA.sync_threads()

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

        return
    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 = CuArray{Float64}(undef, n)
        dl = CuArray{Float64}(undef, n)
        du = CuArray{Float64}(undef, n)
        d_out = CuArray{Float64}(undef, n)
        copyto!(dx, x)
        copyto!(dl, xl)
        copyto!(du, xu)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((3*n)*sizeof(Float64)) dmid_test(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
    function dgpstep_test(n,dx::CuDeviceArray{Float64},
                            dl::CuDeviceArray{Float64},
                            du::CuDeviceArray{Float64},
                            alpha::Float64,
                            dw::CuDeviceArray{Float64},
                            d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        xl = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        xu = CuDynamicSharedArray(Float64, n, (2*n)*sizeof(Float64))
        w = CuDynamicSharedArray(Float64, n, (3*n)*sizeof(Float64))
        s = CuDynamicSharedArray(Float64, n, (4*n)*sizeof(Float64))
        x[tx] = dx[tx]
        xl[tx] = dl[tx]
        xu[tx] = du[tx]
        w[tx] = dw[tx]
        CUDA.sync_threads()

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

        return
    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 = CuArray{Float64}(undef, n)
        dl = CuArray{Float64}(undef, n)
        du = CuArray{Float64}(undef, n)
        dw = CuArray{Float64}(undef, n)
        d_out = CuArray{Float64}(undef, n)
        copyto!(dx, x)
        copyto!(dl, xl)
        copyto!(du, xu)
        copyto!(dw, w)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((5*n)*sizeof(Float64)) dgpstep_test(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
    function dbreakpt_test(n,dx::CuDeviceArray{Float64},
                            dl::CuDeviceArray{Float64},
                            du::CuDeviceArray{Float64},
                            dw::CuDeviceArray{Float64},
                            d_nbrpt::CuDeviceArray{Float64},
                            d_brptmin::CuDeviceArray{Float64},
                            d_brptmax::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        xl = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        xu = CuDynamicSharedArray(Float64, n, (2*n)*sizeof(Float64))
        w = CuDynamicSharedArray(Float64, n, (3*n)*sizeof(Float64))
        x[tx] = dx[tx]
        xl[tx] = dl[tx]
        xu[tx] = du[tx]
        w[tx] = dw[tx]
        CUDA.sync_threads()

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

        return
    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 = CuArray{Float64}(undef, n)
        dl = CuArray{Float64}(undef, n)
        du = CuArray{Float64}(undef, n)
        dw = CuArray{Float64}(undef, n)
        d_nbrpt = CuArray{Float64,2}(undef, (n,n))
        d_brptmin = CuArray{Float64,2}(undef, (n,n))
        d_brptmax = CuArray{Float64,2}(undef, (n,n))
        copyto!(dx, x)
        copyto!(dl, xl)
        copyto!(du, xu)
        copyto!(dw, w)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((4*n)*sizeof(Float64)) dbreakpt_test(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
    function dnrm2_test(n::Int, d_in::CuDeviceArray{Float64},
                        d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        x[tx] = d_in[tx]
        CUDA.sync_threads()

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

        return
    end

    @inbounds for i=1:itermax
        h_in = rand(n)
        h_out = zeros(n)
        d_in = CuArray{Float64}(undef, n)
        d_out = CuArray{Float64}(undef, n)
        copyto!(d_in, h_in)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=(n*sizeof(Float64)) dnrm2_test(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
    function nrm2_test(n::Int, d_A::CuDeviceArray{Float64}, d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        wa = CuDynamicSharedArray(Float64, n)
        A = CuDynamicSharedArray(Float64, (n,n), n*sizeof(Float64))

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

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

        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 = CuArray{Float64,2}(undef, (n,n))
        d_out = CuArray{Float64}(undef, n)
        h_wa = zeros(n)
        copyto!(d_A, A)
        CUDA.@sync @cuda threads=(n,n) blocks=nblk shmem=((n^2+n)*sizeof(Float64)) nrm2_test(n,d_A,d_out)
        copyto!(h_wa, d_out)

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

@testset "dcopy" begin
    function dcopy_test(n::Int, d_in::CuDeviceArray{Float64},
                        d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        y = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))

        x[tx] = d_in[tx]
        CUDA.sync_threads()

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

        if bx == 1
            d_out[tx] = y[tx]
        end
        CUDA.sync_threads()

        return
    end

    @inbounds for i=1:itermax
        h_in = rand(n)
        h_out = zeros(n)
        d_in = CuArray{Float64}(undef, n)
        d_out = CuArray{Float64}(undef, n)
        copyto!(d_in, h_in)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((2*n)*sizeof(Float64)) dcopy_test(n,d_in,d_out)
        copyto!(h_out, d_out)

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

@testset "ddot" begin
    function ddot_test(n::Int, d_in::CuDeviceArray{Float64},
                        d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        y = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        x[tx] = d_in[tx]
        y[tx] = d_in[tx]
        CUDA.sync_threads()

        v = ExaTronKernels.ddot(n, x, 1, y, 1)
        if bx == 1
            for i in 1:n
                d_out[i,tx] = v
            end
        end
        CUDA.sync_threads()

        return
    end

    @inbounds for i=1:itermax
        h_in = rand(n)
        h_out = zeros((n,n))
        d_in = CuArray{Float64}(undef, n)
        d_out = CuArray{Float64,2}(undef, (n,n))
        copyto!(d_in, h_in)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((2*n)*sizeof(Float64)) ddot_test(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
    function dscal_test(n::Int, da::Float64,
                        d_in::CuDeviceArray{Float64},
                        d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        x[tx] = d_in[tx]
        CUDA.sync_threads()

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

        return
    end

    for i=1:itermax
        h_in = rand(n)
        h_out = zeros(n)
        da = rand(1)[1]
        d_in = CuArray{Float64}(undef, n)
        d_out = CuArray{Float64}(undef, n)
        copyto!(d_in, h_in)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=(n*sizeof(Float64)) dscal_test(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
    function dtrqsol_test(n::Int, d_x::CuDeviceArray{Float64},
                            d_p::CuDeviceArray{Float64},
                            d_out::CuDeviceArray{Float64},
                            delta::Float64)
        tx = CUDA.threadIdx().x
        ty = CUDA.threadIdx().y

        x = CuDynamicSharedArray(Float64, n)
        p = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))

        x[tx] = d_x[tx]
        p[tx] = d_p[tx]
        CUDA.sync_threads()

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

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

        d_x = CuArray{Float64}(undef, n)
        d_p = CuArray{Float64}(undef, n)
        d_out = CuArray{Float64,2}(undef, (n,n))
        copyto!(d_x, x)
        copyto!(d_p, p)
        CUDA.@sync @cuda threads=n blocks=nblk shmem=((2*n)*sizeof(Float64)) dtrqsol_test(n,d_x,d_p,d_out,delta)

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

@testset "dspcg" begin
    function dspcg_test(n::Int, delta::Float64, rtol::Float64,
                        cg_itermax::Int, dx::CuDeviceArray{Float64},
                        dxl::CuDeviceArray{Float64},
                        dxu::CuDeviceArray{Float64},
                        dA::CuDeviceArray{Float64},
                        dg::CuDeviceArray{Float64},
                        ds::CuDeviceArray{Float64},
                        d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        xl = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        xu = CuDynamicSharedArray(Float64, n, (2*n)*sizeof(Float64))
        g = CuDynamicSharedArray(Float64, n, (3*n)*sizeof(Float64))
        s = CuDynamicSharedArray(Float64, n, (4*n)*sizeof(Float64))
        w = CuDynamicSharedArray(Float64, n, (5*n)*sizeof(Float64))
        wa1 = CuDynamicSharedArray(Float64, n, (6*n)*sizeof(Float64))
        wa2 = CuDynamicSharedArray(Float64, n, (7*n)*sizeof(Float64))
        wa3 = CuDynamicSharedArray(Float64, n, (8*n)*sizeof(Float64))
        wa4 = CuDynamicSharedArray(Float64, n, (9*n)*sizeof(Float64))
        wa5 = CuDynamicSharedArray(Float64, n, (10*n)*sizeof(Float64))
        gfree = CuDynamicSharedArray(Float64, n, (11*n)*sizeof(Float64))
        indfree = CuDynamicSharedArray(Int, n, (12*n)*sizeof(Float64))
        iwa = CuDynamicSharedArray(Int, 2*n, n*sizeof(Int) + (12*n)*sizeof(Float64))

        A = CuDynamicSharedArray(Float64, (n,n), (12*n)*sizeof(Float64)+(3*n)*sizeof(Int))
        B = CuDynamicSharedArray(Float64, (n,n), (12*n+n^2)*sizeof(Float64)+(3*n)*sizeof(Int))
        L = CuDynamicSharedArray(Float64, (n,n), (12*n+2*n^2)*sizeof(Float64)+(3*n)*sizeof(Int))

        @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]
        CUDA.sync_threads()

        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
        CUDA.sync_threads()

        return
    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(Int, n)
        iwa = zeros(Int, 3*n)

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

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

        CUDA.@sync @cuda threads=n blocks=nblk shmem=((3*n)*sizeof(Int)+(12*n+3*(n^2))*sizeof(Float64)) dspcg_test(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
    function dgpnorm_test(n, dx, dxl, dxu, dg, d_out)
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        xl = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        xu = CuDynamicSharedArray(Float64, n, (2*n)*sizeof(Float64))
        g = CuDynamicSharedArray(Float64, n, (3*n)*sizeof(Float64))

        x[tx] = dx[tx]
        xl[tx] = dxl[tx]
        xu[tx] = dxu[tx]
        g[tx] = dg[tx]
        CUDA.sync_threads()

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

        return
    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 = CuArray{Float64}(undef, n)
        dxl = CuArray{Float64}(undef, n)
        dxu = CuArray{Float64}(undef, n)
        dg = CuArray{Float64}(undef, n)
        d_out = CuArray{Float64}(undef, n)

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

        gptime = @timed CUDA.@sync @cuda threads=n blocks=nblk shmem=(4*n*sizeof(Float64)) dgpnorm_test(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
    function dtron_test(n::Int, f::Float64, frtol::Float64, fatol::Float64, fmin::Float64,
                        cgtol::Float64, cg_itermax::Int, delta::Float64, task::Int,
                        disave::CuDeviceArray{Int}, ddsave::CuDeviceArray{Float64},
                        dx::CuDeviceArray{Float64}, dxl::CuDeviceArray{Float64},
                        dxu::CuDeviceArray{Float64}, dA::CuDeviceArray{Float64},
                        dg::CuDeviceArray{Float64}, d_out::CuDeviceArray{Float64})
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        xl = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        xu = CuDynamicSharedArray(Float64, n, (2*n)*sizeof(Float64))
        g = CuDynamicSharedArray(Float64, n, (3*n)*sizeof(Float64))
        xc = CuDynamicSharedArray(Float64, n, (4*n)*sizeof(Float64))
        s = CuDynamicSharedArray(Float64, n, (5*n)*sizeof(Float64))
        wa = CuDynamicSharedArray(Float64, n, (6*n)*sizeof(Float64))
        wa1 = CuDynamicSharedArray(Float64, n, (7*n)*sizeof(Float64))
        wa2 = CuDynamicSharedArray(Float64, n, (8*n)*sizeof(Float64))
        wa3 = CuDynamicSharedArray(Float64, n, (9*n)*sizeof(Float64))
        wa4 = CuDynamicSharedArray(Float64, n, (10*n)*sizeof(Float64))
        wa5 = CuDynamicSharedArray(Float64, n, (11*n)*sizeof(Float64))
        gfree = CuDynamicSharedArray(Float64, n, (12*n)*sizeof(Float64))
        indfree = CuDynamicSharedArray(Int, n, (13*n)*sizeof(Float64))
        iwa = CuDynamicSharedArray(Int, 2*n, n*sizeof(Int) + (13*n)*sizeof(Float64))

        A = CuDynamicSharedArray(Float64, (n,n), (13*n)*sizeof(Float64)+(3*n)*sizeof(Int))
        B = CuDynamicSharedArray(Float64, (n,n), (13*n+n^2)*sizeof(Float64)+(3*n)*sizeof(Int))
        L = CuDynamicSharedArray(Float64, (n,n), (13*n+2*n^2)*sizeof(Float64)+(3*n)*sizeof(Int))

        @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]
        CUDA.sync_threads()

        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
        CUDA.sync_threads()

        return
    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 = CuArray{Float64}(undef, n)
        dxl = CuArray{Float64}(undef, n)
        dxu = CuArray{Float64}(undef, n)
        dA = CuArray{Float64,2}(undef, (n,n))
        dg = CuArray{Float64}(undef, n)
        disave = CuArray{Int}(undef, n)
        ddsave = CuArray{Float64}(undef, n)
        d_out = CuArray{Float64}(undef, n)

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

        @cuda threads=n blocks=nblk shmem=((3*n)*sizeof(Int)+(13*n+3*(n^2))*sizeof(Float64)) dtron_test(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
    function eval_f(n, x, dA, dc)
        f = 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
        CUDA.sync_threads()
        return f
    end

    function eval_g(n, x, g, dA, dc)
        @inbounds for i=1:n
            gval = 0
            @inbounds for j=1:n
                gval += dA[i,j]*x[j]
            end
            g[i] = gval + dc[i]
        end
        CUDA.sync_threads()
        return
    end

    function eval_h(n, scale, x, A, dA)
        tx = CUDA.threadIdx().x

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

    function driver_kernel(n::Int, max_feval::Int, max_minor::Int,
                            x::CuDeviceArray{Float64}, xl::CuDeviceArray{Float64},
                            xu::CuDeviceArray{Float64}, dA::CuDeviceArray{Float64},
                            dc::CuDeviceArray{Float64})
        # We start with a shared memory allocation.
        # The first 3*n*sizeof(Float64) bytes are used for x, xl, and xu.

        g = CuDynamicSharedArray(Float64, n, (3*n)*sizeof(Float64))
        xc = CuDynamicSharedArray(Float64, n, (4*n)*sizeof(Float64))
        s = CuDynamicSharedArray(Float64, n, (5*n)*sizeof(Float64))
        wa = CuDynamicSharedArray(Float64, n, (6*n)*sizeof(Float64))
        wa1 = CuDynamicSharedArray(Float64, n, (7*n)*sizeof(Float64))
        wa2 = CuDynamicSharedArray(Float64, n, (8*n)*sizeof(Float64))
        wa3 = CuDynamicSharedArray(Float64, n, (9*n)*sizeof(Float64))
        wa4 = CuDynamicSharedArray(Float64, n, (10*n)*sizeof(Float64))
        wa5 = CuDynamicSharedArray(Float64, n, (11*n)*sizeof(Float64))
        gfree = CuDynamicSharedArray(Float64, n, (12*n)*sizeof(Float64))
        dsave = CuDynamicSharedArray(Float64, n, (13*n)*sizeof(Float64))
        indfree = CuDynamicSharedArray(Int, n, (14*n)*sizeof(Float64))
        iwa = CuDynamicSharedArray(Int, 2*n, n*sizeof(Int) + (14*n)*sizeof(Float64))
        isave = CuDynamicSharedArray(Int, n, (3*n)*sizeof(Int) + (14*n)*sizeof(Float64))

        A = CuDynamicSharedArray(Float64, (n,n), (14*n)*sizeof(Float64)+(4*n)*sizeof(Int))
        B = CuDynamicSharedArray(Float64, (n,n), (14*n+n^2)*sizeof(Float64)+(4*n)*sizeof(Int))
        L = CuDynamicSharedArray(Float64, (n,n), (14*n+2*n^2)*sizeof(Float64)+(4*n)*sizeof(Int))

        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, max_feval, max_minor,
                                dx, dxl, dxu, dA, dc, d_out)
        tx = CUDA.threadIdx().x
        bx = CUDA.blockIdx().x

        x = CuDynamicSharedArray(Float64, n)
        xl = CuDynamicSharedArray(Float64, n, n*sizeof(Float64))
        xu = CuDynamicSharedArray(Float64, n, (2*n)*sizeof(Float64))

        x[tx] = dx[tx]
        xl[tx] = dxl[tx]
        xu[tx] = dxu[tx]
        CUDA.sync_threads()

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

        if bx == 1
            d_out[tx] = x[tx]
        end
        CUDA.sync_threads()
        return
    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 = CuArray{Float64}(undef, n)
    dxl = CuArray{Float64}(undef, n)
    dxu = CuArray{Float64}(undef, n)
    dA = CuArray{Float64,2}(undef, (n,n))
    dc = CuArray{Float64}(undef, n)
    d_out = CuArray{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)

        CUDA.@sync @cuda threads=n blocks=nblk shmem=((4*n)*sizeof(Int)+(14*n+3*(n^2))*sizeof(Float64)) driver_kernel_test(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:
"""
    @params struct_def

A macro that changes all fields' types to type parameters while respecting the type bounds specificied by the user. For example:
```
@params struct MyType{T}
    f1::T
    f2::AbstractVector{T}
    f3::AbstractVector{<:Real}
    f4
end
```
will define:
```
struct MyType{T, T1 <: T, T2 <: AbstractVector{T}, T3 <: AbstractVector{<:Real}, T4}
    f1::T1
    f2::T1
    f3::T3
    f4::T4
end
```
The default non-parameteric constructor, e.g. `MyType(f1, f2, f3, f4)`, will always work if all the type parameters used in the type header are used in the field types. Using a type parameter in the type header that is not used in the field types such as:
```
struct MyType{T}
    f1
    f2
end
```
is not recommended.
"""
macro params(struct_expr)
    header = struct_expr.args[2]
    fields = @view struct_expr.args[3].args[2:2:end]
    params = []
    for i in 1:length(fields)
        x = fields[i]
        T = gensym()
        if x isa Symbol
            push!(params, T)
            fields[i] = :($x::$T)
        elseif x.head == :(::)
            abstr = x.args[2]
            var = x.args[1]
            push!(params, :($T <: $abstr))
            fields[i] = :($var::$T)
        end
    end
    if header isa Symbol && length(params) > 0
        struct_expr.args[2] = :($header{$(params...)})
    elseif header.head == :curly
        append!(struct_expr.args[2].args, params)
    elseif header.head == :<:
        if struct_expr.args[2].args[1] isa Symbol
            name = struct_expr.args[2].args[1]
            struct_expr.args[2].args[1] = :($name{$(params...)})
        elseif header.head == :<: && struct_expr.args[2].args[1] isa Expr
            append!(struct_expr.args[2].args[1].args, params)
        else
            error("Unidentified type definition.")
        end
    else
        error("Unidentified type definition.")
    end
    esc(struct_expr)
end

struct RaggedArray{TO, TV}
    offsets::TO
    values::TV
end

function RaggedArray(vv::Vector{Vector{T}}) where T
    offsets = [1; 1 .+ accumulate(+, collect(length(v) for v in vv))]
    values = Vector{T}(undef, offsets[end]-1)
    for (i, v) in enumerate(vv)
        r = offsets[i]:offsets[i+1]-1
        values[r] .= v
    end
    RaggedArray(offsets, values)
end

function Base.getindex(ra::RaggedArray, i)
    @assert 1 <= i < length(ra.offsets)
    r = ra.offsets[i]:ra.offsets[i+1]-1
    @assert 1 <= r.start && r.stop <= length(ra.values)
    return @view ra.values[r]
end
function Base.getindex(ra::RaggedArray, i, j)
    @assert 1 <= j < length(ra.offsets)
    r = ra.offsets[j]:ra.offsets[j+1]-1
    @assert 1 <= i <= length(r)
    return ra.values[r[i]]
end
function Base.setindex!(ra::RaggedArray, v, i, j)
    @assert 1 <= j < length(ra.offsets)
    r = ra.offsets[j]:ra.offsets[j+1]-1
    @assert 1 <= i <= length(r)
    ra.values[r[i]] = v
end

function find_varind(black, white, ::Type{TI}=Int) where TI
    nel = length(black)
    nel == length(white) || throw("Black and white vectors should be of the same length")
    varind = zeros(TI, nel)
    k = 1
    for i in 1:nel
        if !black[i] && !white[i]
            varind[i] = k
            k += 1
        end
    end
    return varind
end

function find_black_and_white(dh)
    black = falses(getncells(dh.grid))
    white = falses(getncells(dh.grid))
    if haskey(dh.grid.cellsets, "black")
        for c in grid.cellsets["black"]
            black[c] = true
        end
    end
    if haskey(dh.grid.cellsets, "white")
        for c in grid.cellsets["white"]
            white[c] = true
        end
    end

    return black, white
end

YoungsModulus(p) = getE(p)
PoissonRatio(p) = getν(p)

function compliance(Ke, u, dofs)
    comp = zero(eltype(u))
    for i in 1:length(dofs)
        for j in 1:length(dofs)
            comp += u[dofs[i]]*Ke[i,j]*u[dofs[j]]
        end
    end
    comp
end

function meandiag(K::AbstractMatrix)
    z = zero(eltype(K))
    for i in 1:size(K, 1)
        z += abs(K[i, i])
    end
    return z / size(K, 1)
end

density(var, xmin) = var*(1-xmin) + xmin

macro debug(expr)
    return quote
        if DEBUG[]
            $(esc(expr))
        end
    end
end

@generated function _getproperty(c::T, ::Val{fallback}, ::Val{f}) where {T, fallback, f}
    f ∈ fieldnames(T) && return :(getfield(c, $(QuoteNode(f))))
    return :(getproperty(getfield(c, $(QuoteNode(fallback))), $(QuoteNode(f))))
end
@generated function _setproperty!(c::T, ::Val{fallback}, ::Val{f}, val) where {T, fallback, f}
    f ∈ fieldnames(T) && return :(setfield!(c, $(QuoteNode(f)), val))
    return :(setproperty!(getfield(c, $(QuoteNode(fallback))), $(QuoteNode(f)), val))
end
macro forward_property(T, field)
    quote
        Base.getproperty(c::$(esc(T)), f::Symbol) = _getproperty(c, Val($(QuoteNode(field))), Val(f))
        Base.setproperty!(c::$(esc(T)), f::Symbol, val) = _setproperty!(c, Val($(QuoteNode(field))), Val(f), val)
    end
end

ext: jl


text:
# Set the Attributes of some "builtin" Protected symbols

## NOTE!: Until the attribute setting is cleaned up, in particular in math_functions.jl,
## This file should be loaded last or almost last in order to correct incorrectly set attributes.
## This happens because big lists of functions are handled by macros that set default attributes,
## but the defaults should not apply to a few of them.

## TODO: The use of arrays, vs. single elements and strings vs symbols is outdated
## set_sysattributes takes two args. either can be a string or a symbol or an array of strings or symbols.

set_sysattributes(["Pattern", "SetJ", "SetAttributes", "ClearAttributes", "TimesBy", "AddTo", "Catch"], :HoldFirst)

## Piecewise is also read protected in Mma
set_sysattributes(["Module","LModule", "Clear", "ClearAll", "Condition", "HoldPattern", "HoldForm", "Hold",
                   "DumpHold", "DownValues", "UpValues", "HAge", "Table", "NTable", "For", "If", "While", "Do",
                   "Jxpr", "JuliaExpression", "Protect", "Unprotect", "Function", "Definition", "DirtyQ", "ToSymata",
                   "NIntegrate", "Compile","SymataCall", "Piecewise", "ExportJ", "OwnValues"],
                  :HoldAll)

set_sysattributes(["HoldComplete", "Unevaluated"], [:HoldAllComplete])

set_sysattributes(["Save", "Last"], [:HoldRest] )

set_sysattributes("Attributes",[:HoldAll,:Listable])

set_sysattributes("DivisorSigma",[:NHoldAll,:Listable])

set_sysattributes("Rule",:SequenceHold)

set_sysattributes(["Set","UpSet"],[:HoldFirst, :SequenceHold])

set_sysattributes(["Increment","Decrement"],[:HoldFirst, :ReadProtected])

set_sysattributes(["RuleDelayed","PatternTest"],[:HoldRest, :SequenceHold])

set_sysattributes(["Timing","Allocated","SetDelayed","UpSetDelayed"], [:HoldAll,:SequenceHold])

set_sysattributes(["Pi","E"],[:ReadProtected,:Constant])
set_sysattributes(["EulerGamma", "GoldenRatio", "Catalan"],[:Constant])

set_sysattributes("I", [:ReadProtected,:Locked]) # Careful with this. We mostly use julia symbol :I bound to complex(0,1)

set_sysattributes(["CompoundExpression","Sum","Product"],[:ReadProtected,:HoldAll])

set_sysattributes(["Part","D","LaplaceTransform","InverseLaplaceTransform",
                 "FourierTransform","InverseFourierTransform", "Integrate", "DSolve", "ReadString",
                 "STDOUT", "STDERR", "DevNull", "HoldAllComplete"
                 ],:ReadProtected)

# We kinda need Exp, see the apprules.
set_sysattributes(["Cos", "ArcCos", "Sin", "ArcSin", "Tan", "ArcTan", "CosPi", "SinPi",
                 "Cot", "Cosh", "Conjugate", "Sinh","Minus","Abs","Re","Im", "ReIm", "Exp", "Sqrt",
                 "PolyGamma", "EllipticE", "EllipticF", "EllipticK", "EllipticPi", "LogIntegral", "Mod", "DivRem",
                 "Sign", "SphericalHarmonicY", "SphericalBesselJ", "SphericalBesselY", "Erf", "Gamma", "GammaRegularized", "Divide",
                 "DirichletEta", "PolyLog", "Sqrt", "CubeRoot", "Surd", "BesselJZero"
                 ],
                [:Listable,:NumericFunction])

## These were given Listable attribute in math_functions.jl
clear_attributes(:MeijerG, :HypergeometricPFQ)
set_sysattributes([:MeijerG, :HypergeometricPFQ], [:NumericFunction])

set_sysattributes([:HeavisideTheta],
                [:Listable, :Orderless])

set_sysattributes(["Plus", "Times"],
                [:Flat,:Listable,:NumericFunction,:OneIdentity,:Orderless])

set_sysattributes(["Max", "Min"],
                [:Flat,:NumericFunction,:OneIdentity,:Orderless])

set_sysattributes(["LCM", "GCD"],
                [:Flat,:Listable,:OneIdentity,:Orderless])

set_sysattributes(["Composition"],
            [:Flat,:OneIdentity,])

set_sysattributes("Power",[:Listable,:NumericFunction,:OneIdentity])

# This is not quite what Mma has for these. I don't understand why. Eg. CosIntegral is not a NumericFunction
# Many are not Listable in the docs, although they are in practice.
set_sysattributes(["Pochhammer", "LogIntegral", "LerchPhi", "CosIntegral", "SinIntegral", "FresnelC", "FresnelS", "MittagLefflerE",
                 "HarmonicNumber", "InverseErf", "Log"],
                [:Listable,:NumericFunction,:ReadProtected])

# BellB should be split into BellB and BellY
clear_attributes(:BellB)  ## FIXME: set these attributes only once
set_sysattributes(["BellB"],
                [:NumericFunction,:ReadProtected])

# Note: Simplify is not Listable in Mma
set_sysattributes(["Boole", "EvenQ","OddQ", "PrimeQ", "Range","Limit", "Together", "Apart", "Cyclotomic", "MoebiusMu", "EulerPhi", "Divisors", "DivisorCount",
                 "ToExpression", "Simplify", "FullSimplify", "StringReverse", "ToJExpression", "LetterNumber"]
                ,[:Listable])

set_sysattributes(["DirectedInfinity"], [:Listable, :ReadProtected])

# Some symbols that may be particular to SymPy
# We should use GaussianIntegers instead of Gaussian
set_sysattributes(["PolarLift", "ExpPolar", "ExpandFunc", "Modulus","Force", "Deep", "Gaussian"],
                :Protected)

set_sysattributes(["And", "Or"], [:Flat, :HoldAll, :OneIdentity])

set_sysattributes(["UnicodeForm", "JupyterForm", "InputForm"], :Protected)

## Do we still need these ?
#set_sysattributes(["Plain","Unicode", "IJulia" ], :Protected)

set_sysattributes(["Return","Break","Continue"], :Protected)

set_sysattributes(["HoldFirst","HoldAll","HoldRest","NHoldFirst","NHoldAll","NHoldRest"], :Protected)

# FIXME: Subscript cannot have Attribute 'Protected', because the user must be able to assign, e.g. Subscript(a,1) = 1.
# But, without Protected, we do not get TAB completion.
set_attributes(:Subscript, :NHoldRest)

set_sysattributes(["Int64", "Int", "Float64",  "BigFloat", "BigInt", "AbstractString", "Real", "Float"], :Protected)

set_sysattributes(["O","OO","OOO","OOOO","OOOOO","OOOOOO","OOOOOOO","OOOOOOOO","OOOOOOOOOO"], :Protected)

set_sysattributes(["System", "Main"], :Protected)

# True and False are actually parsed as Julia true and false. But, this is a quick way to get TAB completion for them
set_sysattributes(["True", "False"], :Protected)

set_sysattributes(["Optional","Repeated", "RepeatedNull" ], :Protected)

set_sysattributes([
    "!=","//","<","<=","==",">","===",">=","Age","Algebraics","All","Alternatives","Apply","Array","ArrayDepth",
    "AtomQ","BF","BI","Big","Binary","Blank",
    "BlankNullSequence","BlankSequence","Booleans","BuiltIns","ByteCount","Cancel",
    "Chop", "Collect","Comparison","CompiledFunction","Complex", "ComplexInfinity",
    "Complexes","ConditionalExpression",
    "Conditions","ConstantArray","ConstantQ","Depth","DirtyQ","Dump","Evaluate","Example","Except",
    "ExpToTrig", "Expand", "ExpandA", "ExpandCos", "ExpandCosRule",
    "ExpandSin","ExpandSinCos","ExpandSinRule","Factor","FactorInteger",
    "FixedQ","Flat","FullForm","Head", "Help",
    "Indeterminate","Infinity","Integer","IntegerDigits","Integers",
    "JVar","Keys","LeafCount","Length","List","ListQ","Listable",
    "Map","MatchQ","Missing","N","None","Not",
    "Null","Numerator","NumericFunction","NumericQ","OneIdentity","Order","Orderless",
    "Pack","Permutations", "PossibleClosedForm","Precision","Primes","Print",
    "PrintTests","Println", "Protected",
    "RatSimp","Rational","Rationals","RealRoots","Reals", "Replace", "ReplaceAll",
    "ReplaceRepeated","Reverse","Roots","Sequence","SequenceHold",
    "Series","Solve","Span","String","StringJoin","StringLength","StringQ","Subdivide","Symbol",
    "Syms","Temporary","TensorRank","Throw",
    "Time","ToLowerCase","ToString","ToSymPy","ToUpperCase", "Trace","TraceDownValues",
    "TraceUpValues","TrigSimp", "Undefined", "Unfix", "Unpack", "UserSyms", "Values", "Normal",
    "ans", "nothing" ], :Protected)


text:
""" Checks if `l` is a typle, list or vector type """
isa_list(l) = l isa Vector || l isa Tuple || l isa Array


text:
import Libdl: dlext

DIR = @__DIR__

if Sys.isapple()
    const DEP_DIR = "/usr/local/opt/librdkafka/lib"
elseif DIR != nothing
    const DEP_DIR = joinpath(DIR, "..", "deps", "usr")
else
    const DEP_DIR = Pkg.dir("RDKafka", "deps", "usr")
end

function find_librdkafka(dep_dir)
    base = "librdkafka"
    for fname in [
        "$base.$dlext",
        "$base.$dlext.1",
        "$base.1.$dlext",
    ]
        path = abspath(joinpath(dep_dir, fname))
        if isfile(path)
            return path
        end
    end
    return nothing
end

const LIBRDKAFKA = find_librdkafka(DEP_DIR)


text:
#  Copyright 2022 California Institute of Technology
#
#  Licensed under the Apache License, Version 2.0 (the "License");
#  you may not use this file except in compliance with the License.
#  You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
#  Unless required by applicable law or agreed to in writing, software
#  distributed under the License is distributed on an "AS IS" BASIS,
#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#  See the License for the specific language governing permissions and
#  limitations under the License.
#
# Author: Philip G. Brodrick, philip.g.brodrick@jpl.nasa.gov

using DataFrames
using CSV
using Interpolations
using Logging
using Statistics
using Plots
using ModularIndices

function nanargmax(input::Array)
    x = copy(input)
    x[isnan.(x)] .= -Inf;
    return argmax(x);
end

function nanargmin(input::Array)
    x = copy(input)
    x[isnan.(x)] .= Inf;
    return argmin(x);
end

function read_envi_wavelengths(filename::String)
    header_name = splitext(filename)[1] * ".hdr"
    header = readlines(header_name)
    found = false
    for line in header
        if occursin("wavelength = {", line) || occursin("wavelength= {", line) ||  occursin("wavelength={", line)
            header = line
            found = true
            break
        end
    end

    if !found
        @error "No wavelength found in " * header_name
        return nothing
    end

    wavelengths = [parse(Float64, strip(x)) for x in split(split(split(header, "{")[2], "}")[1],",")]
    return wavelengths
end

function get_good_bands_mask(wavelengths::Array{Float64}, wavelength_pairs)
    good_bands = ones(Bool, length(wavelengths))

    for wvp in wavelength_pairs
        wavelength_diff = wavelengths .- wvp[1]
        wavelength_diff[wavelength_diff .< 0] .= maximum(filter(!isnan, wavelength_diff))
        lower_index = nanargmin(wavelength_diff)

        wavelength_diff = wvp[2] .- wavelengths
        wavelength_diff[wavelength_diff .< 0] .= maximum(filter(!isnan, wavelength_diff))
        upper_index = nanargmin(wavelength_diff)
        good_bands[lower_index:upper_index] .= false
    end

    return good_bands
end

mutable struct SpectralLibrary
    file_name::String
    class_header_name::String
    spectral_starting_column::Int64
    truncate_end_columns::Int64
    class_valid_keys
    scale_factor::Float64
    wavelength_regions_ignore
    SpectralLibrary(file_name::String, class_header_name::String, spectral_starting_column::Int64, truncate_end_columns::Int64, class_valid_keys = nothing,
                    scale_factor = 1.0, wavelength_regions_ignore= [[0,440],[1310,1490],[1770,2050],[2440,2880]]) =
                    new(file_name, class_header_name, spectral_starting_column, truncate_end_columns, class_valid_keys, scale_factor, wavelength_regions_ignore)

    spectra
    classes
    good_bands
    wavelengths
end

function load_data!(library::SpectralLibrary)

    df = DataFrame(CSV.File(library.file_name))

    try
        library.spectra = Matrix(df[:,library.spectral_starting_column:end-library.truncate_end_columns])
    catch e
        throw(ArgumentError("Could not read spectra from endmember library.  Try adjusting spectral_starting_column or truncate_end_columns, or reworking the library."))
    end

    try
        library.classes = convert(Array{AbstractString}, df[!,library.class_header_name])
    catch e
        throw(ArgumentError("Could not read classes from endmember library.  Try adjusting class_header_name, or reworking the library."))
    end

    try
        library.wavelengths = parse.(Float64, names(df[:,library.spectral_starting_column:end-library.truncate_end_columns]))
    catch e
        throw(ArgumentError("Could not read wavelengths from endmember library.  Try adjusting class_header_name, or reworking the library."))
    end

    wl_order = sortperm(library.wavelengths)
    library.spectra = library.spectra[:,wl_order]
    library.wavelengths = library.wavelengths[wl_order]

    if isnothing(library.class_valid_keys)
        library.class_valid_keys = unique(library.classes)
    end

    good_bands = get_good_bands_mask(library.wavelengths, library.wavelength_regions_ignore)
    library.good_bands = good_bands

    return library
end

function filter_by_class!(library::SpectralLibrary)
    if isnothing(library.class_valid_keys)
        @info "No class valid keys provided, no filtering occuring"
        return
    end

    valid_classes = zeros(Bool, size(library.spectra)[1])
    for cla in library.class_valid_keys
        valid_classes[library.classes .== cla] .= true
    end

    library.spectra = library.spectra[valid_classes,:]
    library.classes = library.classes[valid_classes]
end

function remove_wavelength_region_inplace!(library::SpectralLibrary, set_as_nans::Bool=false)
    if set_as_nans
        library.spectra[:,.!library.good_bands] .= NaN
        library.wavelengths[.!library.good_bands] .= NaN
    else
        library.spectra = library.spectra[:, library.good_bands]
        library.wavelengths = library.wavelengths[library.good_bands]
    end
end

function interpolate_library_to_new_wavelengths!(library::SpectralLibrary, new_wavelengths::Array{Float64})
    old_spectra = copy(library.spectra)

    library.spectra = zeros((size(library.spectra)[1], length(new_wavelengths)))
    for _s in 1:size(old_spectra)[1]
        fit = LinearInterpolation(library.wavelengths, old_spectra[_s,:], extrapolation_bc=Flat());
        library.spectra[_s,:] = fit(new_wavelengths)
    end
    library.wavelengths = new_wavelengths

    good_bands = get_good_bands_mask(library.wavelengths, library.wavelength_regions_ignore)
    library.good_bands = good_bands
end

function scale_library!(library::SpectralLibrary, scaling_factor=nothing)
    if isnothing(scaling_factor)
        library.spectra =  library.spectra ./ library.scale_factor
    else
        library.spectra /= scaling_factor
    end
end

function brightness_normalize!(library::SpectralLibrary)
    library.spectra = library.spectra ./ sqrt.(mean(library.spectra[:,library.good_bands].^2, dims=2))
end

function plot_mean_endmembers(endmember_library::SpectralLibrary, output_name::String)
    plots = []
    for (_u, u) in enumerate(endmember_library.class_valid_keys)
        mean_spectra = mean(endmember_library.spectra[endmember_library.classes .== u,:],dims=1)[:]
        plot!(endmember_library.wavelengths, mean_spectra, label=u, xlim=[300,3200])
    end
    xlabel!("Wavelength [nm]")
    ylabel!("Reflectance")
    xticks!([500, 1000, 1500, 2000, 2500, 3000])
    plot!(dpi=400)
    savefig(output_name)
end

function plot_endmembers(endmember_library::SpectralLibrary, output_name::String)

    for (_u, u) in enumerate(endmember_library.class_valid_keys)
        if _u == 1
            plot(endmember_library.wavelengths, endmember_library.spectra[endmember_library.classes .== u,:]', lab="", xlim=[300,3200], color=palette(:tab10)[Mod(_u)],dpi=400)
        else
            plot!(endmember_library.wavelengths, endmember_library.spectra[endmember_library.classes .== u,:]', lab="",xlim=[300,3200], color=palette(:tab10)[Mod(_u)])
        end
    end
    xlabel!("Wavelenth [nm]")
    ylabel!("Reflectance")
    xticks!([500, 1000, 1500, 2000, 2500, 3000])
    for (_u, u) in enumerate(endmember_library.class_valid_keys)
        plot!([1:2],[0,0.3], color=palette(:tab10)[Mod(_u)], labels=u)
    end
    savefig(output_name)
end

function plot_endmembers_individually(endmember_library::SpectralLibrary, output_name::String)
    plots = []
    spectra = endmember_library.spectra
    classes = endmember_library.classes
    for (_u, u) in enumerate(endmember_library.class_valid_keys)
        sp = spectra[classes .== u,:]
        sp[broadcast(isnan,sp)] .= 0
        brightness = sum(sp, dims=2)
        toplot = spectra[classes .== u,:] ./ brightness
        #push!(plots, plot(endmember_library.wavelengths, toplot', title=u, color=palette(:tab10)[_u], xlabel="Wavelength [nm]", ylabel="Reflectance"))
        push!(plots, plot(endmember_library.wavelengths, toplot', title=u, xlabel="Wavelength [nm]", ylabel="Reflectance"))
        xticks!([500, 1000, 1500, 2000, 2500])
    end
    plot(plots...,size=(1000,600),dpi=400)
    savefig(output_name)
end


text:
using GLVisualize, GeometryTypes, Reactive, GLAbstraction, Colors

if !isdefined(:runtests)
    window = glscreen()
    timesignal = loop(linspace(0f0,1f0, 360))
end

description = """
Animating unicode particles.
"""

const S = -5f0
const W = 10f0
const N = 1000
const ps = (rand(Point3f0, N)*W)+S
const velocity = rand(Vec3f0, N)*0.01f0
const gravity  = Vec3f0(0,0,-0.04)

upper_bound(x) = x>S+W
lower_bound(x) = x<S
function let_it_snow(position, t)
    @inbounds for i=1:length(ps)
        pos = Vec(position[i])
        p = Point3f0(pos+gravity+velocity[i])
        if any(upper_bound, p) || any(lower_bound, p)
            position[i] = Point3f0(rand(linspace(S,S+W, 1000)),rand(linspace(S,S+W, 1000)), S+W)
            velocity[i] = Vec3f0(0)
        else
            position[i] = p
        end
    end
    position
end
particles       = foldp(let_it_snow, ps, timesignal)
rotation 		= map(rotationmatrix_z, const_lift(*, timesignal, 2f0*pi))
color_ramp      = colormap("Blues", 50)
colors          = RGBA{Float32}[color_ramp[rand(1:50)] for i=1:N]

snowflakes = visualize(
    ('❄', particles),
    color = colors,
    scale = Vec2f0(0.6),
    billboard = true, model = rotation
)

_view(snowflakes, window, camera=:perspective)

if !isdefined(:runtests)
    renderloop(window)
end


text:
# A brief test script for the JCM

include("jcm.jl")

function testplanting()
    tr1 = jcm.Tree(jcm.Species(1), 0, 0)
    tr2 = jcm.Tree(jcm.Species(2), 5, 3)
    tr3 = jcm.Tree(jcm.Species(3), 4, -1)
    jcm.planttree!(tr1)
    jcm.planttree!(tr2)
    jcm.planttree!(tr3)
    println(jcm.forest)
end

function testcompetition()
    # Expected outcome:
    # Update 1: tr2 & tr3 killed (overlapped by tr1)
    # Update 2: tr1 killed (overlapped by tr4)
    # Update 3: tr4 & tr5 survive
    tr1 = jcm.Tree(jcm.Species(1), 0, 0)
    tr2 = jcm.Tree(jcm.Species(2), -2, 3)
    tr3 = jcm.Tree(jcm.Species(3), 0, -2)
    tr4 = jcm.Tree(jcm.Species(4), 4, 7)
    tr5 = jcm.Tree(jcm.Species(5), 12, 27)
    tr1.size = 5
    tr2.size = 3
    tr3.size = 1
    tr4.size = 11
    tr5.size = 9
    jcm.planttree!(tr1)
    jcm.planttree!(tr2)
    jcm.planttree!(tr3)
    jcm.planttree!(tr4)
    jcm.planttree!(tr5)
    jcm.recordindividual(tr1)
    jcm.recordindividual(tr2)
    jcm.recordindividual(tr3)
    jcm.recordindividual(tr4)
    jcm.recordindividual(tr5)
    jcm.settings["verbosity"] = "Debug"
    jcm.settings["datafile"] = "jcm_test_data.csv"
    jcm.settings["datafreq"] = 1
    jcm.run(3, true)
    println(jcm.forest)
end

max_stars_repo_name: veddox/jcm

max_issues_repo_name: veddox/jcm

max_forks_repo_name: veddox/jcm


text:
using Random, Test
using NDTools
using PaddedViews
Random.seed!(42)

# include("utils.jl")
include("offset_types.jl")
include("MutablePaddedViews.jl")
include("datatype_tools.jl")
include("size_tools.jl")
include("iteration_tools.jl")
include("selection_tools.jl")
include("generation_tools.jl")

return


text:
## DSP.jl ##

mutable struct DFTSetup{T}
    setup::Ptr{Cvoid}
    direction::Int

    function DFTSetup(T::DataType, setup::Ptr{Cvoid}, direction::Int)
        dftsetup = new{T}(setup, direction)
        finalizer(plan_destroy, dftsetup)
        dftsetup
    end
end

mutable struct Biquad{T}
    setup::Ptr{Cvoid}
    sections::Int

    function Biquad(T::DataType, setup::Ptr{Cvoid}, sections::Int)
        biquadsetup = new{T}(setup, sections)
        finalizer(biquaddestroy, biquadsetup)
        biquadsetup
    end
end

## === FUNCTIONS == ##

for (T, suff) in ((Float64, "D"), (Float32, ""))

    """
    Convolution between an input Vector{T} 'X', and a kernel/filter Vector{T} 'K'.
    Result vector should have at least length(X) + length(K) - 1 elements

    Returns: 'result'. Computation result is also stored in 'result' argument.
    """
    @eval begin
        function conv!(result::Vector{$T}, X::Vector{$T}, K::Vector{$T})
            ksize = length(K)
            xsize = length(X)
            rsize = length(result)
            if (rsize < xsize + ksize - 1)
                error("'result' must have at least length(X) + length(K) - 1 elements")
            end
            xpadded::Vector{$T} = [zeros($T, ksize - 1); X; zeros($T, ksize)]
            fptr = get_fptr($(string("vDSP_conv", suff)))
            ccall(fptr, Cvoid,
                (Ptr{$T}, Int64, Ptr{$T}, Int64, Ptr{$T}, Int64, UInt64, UInt64),
                xpadded, 1, pointer(K, ksize), -1, result, 1, rsize, ksize)
            return result
        end
    end

    """
    Convolution between an input Vector{T} 'X', and a kernel/filter Vector{T} 'K'.

    Returns: Vector{T} with length = length(X) + length(K) - 1
    """
    @eval begin
        function conv(X::Vector{$T}, K::Vector{$T})
            result = Array{$T}(undef, length(X) + length(K) - 1)
            conv!(result, X, K)
        end
    end

    """
    In-place convolution between an input Vector{T} 'X', and a kernel/filter Vector{T} 'K'.

    Returns: 'X'. 'X' is overwritten with computation result.
    """
    @eval begin
        function conv!(X::Vector{$T}, K::Vector{$T})
            conv!(X, X, K)
        end
    end

    """
    Cross-correlation of two Vector{T}'s 'X' and 'Y'.
    Result vector should have at least length(X) + length(Y) - 1 elements

    Returns: 'result'. The result of the computation is also stored in 'result'
    """
    @eval begin
        function xcorr!(result::Vector{$T}, X::Vector{$T}, Y::Vector{$T})
            ysize = length(Y)
            xsize = length(X)
            rsize = length(result)
            if (rsize < xsize + ysize - 1)
                error("'result' must have at least length(X) + length(Y) - 1 elements")
            end
            xpadded::Vector{$T} = [zeros($T, ysize - 1); X; zeros($T, ysize)]
            fptr = get_fptr($(string("vDSP_conv", suff)))
            ccall(fptr, Cvoid,
                (Ptr{$T}, Int64, Ptr{$T}, Int64, Ptr{$T}, Int64, UInt64, UInt64),
                xpadded, 1, Y, 1, result, 1, rsize, ysize)
            return result
        end
    end

    """
    Cross-correlation of two Vector{T}'s 'X' and 'Y'.

    Returns: Vector{T} with length(X) + length(Y) - 1
    """
    @eval begin
        function xcorr(X::Vector{$T}, Y::Vector{$T})
            result = Array{$T}(undef, length(X) + length(Y) - 1)
            xcorr!(result, X, Y)
        end
    end

    """
    In-place cross-correlation of two Vector{T}'s 'X' and 'Y'.

    Returns: 'X'. 'X' is overwritten with the result of the cross-correlation.
    """
    @eval begin
        function xcorr!(X::Vector{$T}, Y::Vector{$T})
            xcorr!(X, X, Y)
        end
    end

    """
    Performs auto-correlation of a 1-D signal with itself.

    Returns: Vector{T} with length = 2*length(X) - 1
    """
    @eval begin
        function xcorr(X::Vector{$T})
            xcorr(X, X)
        end
    end
end

## == Biquadratic/IIR filtering
for (T, suff) in ((Float64, "D"),)

    """
    Initializes a vDSP_biquad_setup for use with vDSP_biquad. A multi-section filter
    can be initialized with a single call to biquad_create_setup. coefficients must
    contain 5 coefficients for each section. The three feed-forward coefficients are
    specified first, followed by the two feedback coefficients. Returns a Biquad object.

    Returns: Biquad{T}
    """
    @eval begin
        function biquadcreate(coefficients::Vector{$T}, sections::Int)
            if length(coefficients) < 5 * sections
                error("Incomplete biquad specification provided - coefficients must
                            contain 5 elements for each filter section")
            end
            fptr = get_fptr($(string("vDSP_biquad_CreateSetup", suff)))
            setup = ccall(fptr, Ptr{Cvoid},
                (Ptr{$T}, UInt64),
                coefficients, sections)
            return Biquad($T, setup, sections)
        end
    end

    """
    Filters an input array X with the specified Biquad filter and filter delay values provided
    in delays; only numelem elements are filtered. After execution, delays contains the final
    state data of the filter.

    Returns: Vector{T}
    """
    @eval begin
        function biquad(X::Vector{$T}, delays::Vector{$T}, numelem::Int, biquad::Biquad{$T})
            if length(delays) < (2 * (biquad.sections) + 2)
                error("Incomplete delay specification provided - delays must contain 2*M + 2
                                values where M is the number of sections in the biquad")
            end
            result::Vector{$T} = similar(X)
            fptr = get_fptr($(string("vDSP_biquad", suff)))
            ccall(fptr, Cvoid,
                (Ptr{Cvoid}, Ptr{$T}, Ptr{$T}, Int64, Ptr{$T}, Int64, UInt64),
                biquad.setup, delays, X, 1, result, 1, numelem)
            return result
        end
    end

    """
    Frees all resources associated with a particular Biquad previously
    created through a call to biquad_create_setup. This is called automatically
    when the setup object is no longer visible to the garbage collector.

    Returns: Cvoid
    """
    @eval begin
        function biquaddestroy(biquad::Biquad{$T})
            fptr = get_fptr($(string("vDSP_biquad_DestroySetup", suff)))
            ccall(fptr, Cvoid,
                (Ptr{Cvoid},),
                biquad.setup)
        end
    end
end

## == WINDOW GENERATION == ##

"""
    blackman(length, [rtype=Float64])

Generates a Blackman window of length 'length'. Default return type
is Vector{Float64}, but if rtype=Float32, Vector{Float32}
will be returned.
"""
function blackman(length::Int, rtype::DataType = Float64)
    result::Vector{rtype} = Array{rtype}(undef, length)
    blackman!(result, length, 0)
end

"""
    hamming(length, [rtype=Float64])

Generates a Hamming window of length 'length'. Default return type
is Vector{Float64}, but if rtype=Float32, Vector{Float32}
will be returned.
"""
function hamming(length::Int, rtype::DataType = Float64)
    result::Vector{rtype} = Array{rtype}(undef, length)
    hamming!(result, length, 0)
end

"""
    hanning(length, [rtype=Float64])

Generates a denormalized Hanning window of length 'length'. Default
return type is Vector{Float64}, but if rtype=Float32, Vector{Float32}
will be returned.
"""
function hanning(length::Int, rtype::DataType = Float64)
    result::Vector{rtype} = Array{rtype}(undef, length)
    hanning!(result, length, 0)
end

"""
Alias function for `hanning`
"""
function hann(length::Int, rtype::DataType = Float64)
    hanning(length, rtype)
end

for (T, suff) in ((Float32, ""), (Float64, "D"))

    """
    Generates a Blackman window of length 'length' and stores it in `result'. By
    default, a full window is created, but if flag=1, only the first (n+1)/2 points
    will be calculated.

    Returns: Vector{$T}
    """
    @eval begin
        function blackman!(result::Vector{$T}, length::Int, flag::Int = 0)
            fptr = get_fptr($(string("vDSP_blkman_window", suff)))
            ccall(fptr, Cvoid,
                (Ptr{$T}, UInt64, Int64),
                result, length, flag)
            return result
        end
    end

    """
    Generates a Hamming window of length 'length' and stores it in `result'. By
    default, a full window is created, but if flag=1, only the first (n+1)/2 points
    will be calculated.

    Returns: Vector{$T}
    """
    @eval begin
        function hamming!(result::Vector{$T}, length::Int, flag::Int = 0)
            fptr = get_fptr($(string("vDSP_hamm_window", suff)))
            ccall(fptr, Cvoid,
                (Ptr{$T}, UInt64, Int64),
                result, length, flag)
            return result
        end
    end

    """
    Generates a Hanning window of length 'length' and stores it in `result'.
    Different window options can be set using the flag argument; default value
    is zero (denormalized Hanning window)

    flag = 0: Returns a denormalized Hanning window
    flag = 1: Returns a denormalized Hanning window with only (n+1)/2 points
    flag = 2: Returns a normalized Hanning window
    flag = 3: Returns a normalized Hanning window with only (n+1)/2 points

    Returns: Vector{$T}
    """
    @eval begin
        function hanning!(result::Vector{$T}, length::Int, flag::Int = 0)
            fptr = get_fptr($(string("vDSP_hann_window", suff)))
            ccall(fptr, Cvoid,
                (Ptr{$T}, UInt64, Int64),
                result, length, flag)
            return result
        end
    end

    """
    Alias function for hanning!

    Returns: Vector{$T}
    """
    @eval begin
        function hann!(result::Vector{$T}, length::Int, flag::Int = 0)
            hanning!(result, length, flag)
        end
    end

end

## == Discrete Cosine Transform (DCT) == ##
"""
Initializes a new DCT setup object. 'dct_type' must be 2, 3, 4 corresponding to Type II, III and IV.
DCT 'length' must be equal to f*(2^n) where f = 1,3,5,15 and n >= 4. If you have a previous DCT setup
object, that can be passed in as 'previous'. The returned DCT setup will share the underlying data
storage of the previous setup object.

Returns: DFTSetup
"""
function plan_dct(length::Int, dct_type::Int, previous = C_NULL)
    n = trailing_zeros(length)
    f = length >> n
    if dct_type < 2 && dct_type > 4
        error("DCT type ", dct_type, " is not supported. Only DCT types 2, 3 and 4 are supported")
    elseif !(n >= 4 && f in (1, 3, 5, 15))
        error("Invalid DCT length. Length must be equal to f*(2^n) where f = 1,3,5,15 and n >= 4")
    end
    fptr = get_fptr("vDSP_DCT_CreateSetup")
    setup::Ptr{Cvoid} = ccall(fptr, Ptr{Cvoid},
        (Ptr{Cvoid}, UInt64, UInt64),
        previous, length, dct_type)
    return DFTSetup(Float32, setup, 0)
end

"""
Computes the DCT of a given input vector X using the parameters setup in
the DFTSetup object.

Returns: Vector{Float32}
"""
function dct(X::Vector{Float32}, setup::DFTSetup)
    result = similar(X)
    fptr = get_fptr("vDSP_DCT_Execute")

    ccall(fptr, Cvoid,
        (Ptr{Cvoid}, Ptr{Float32}, Ptr{Float32}),
        setup.setup, X, result)
    return result
end

"""
Computes the DCT of a given input vector X using a DCT Type 'dct_type' (defaults to type II).
This function does not require a separate call to dct_setup.

Returns: Vector{Float32}
"""
function dct(X::Vector{Float32}, dct_type::Int = 2)
    setup = plan_dct(length(X), dct_type)
    return dct(X, setup)
end

"""
Deinitializes a DFTSetup object created by plan_dct
"""
function plan_destroy(setup::DFTSetup)
    fptr = get_fptr("vDSP_DFT_DestroySetup")
    ccall(fptr, Cvoid,
        (Ptr{Cvoid},),
        setup.setup)
end

# const     FFT_FORWARD         = +1
# const     FFT_INVERSE         = -1
# const     SIGNAL_STRIDE       = 1
#
# typealias vDSP_Stride Clong
# typealias vDSP_Length Culong
# typealias FFTDirection Cint
#
# immutable DSPDoubleSplitComplex
#     realp::Ptr{Float64}
#     imagp::Ptr{Float64}
# end
#
# DSPDoubleSplitComplex(realp::Vector{Float64},imagp::Vector{Float64})=DSPDoubleSplitComplex(pointer(realp),pointer(imagp))
#
#
# function plan_fft(n,radix::Integer=2)
#     @assert isinteger(log2(n))
#     logn=round(Int,log2(n))
#     ccall(("vDSP_create_fftsetupD",libacc),Ptr{Cvoid},(Cuint,Cint),logn,radix)
# end
#
# function newfft(r::Vector{Complex128},plan)
#     n=length(r)
#     @assert isinteger(log2(n))
#     logn=round(Int,log2(n))
#
#     realp,imagp=real(r),imag(r)         # keep references to avoid garbage collection
#     vals=DSPDoubleSplitComplex(realp,imagp)
#     retr,reti=Array{Float64}(undef,n),Array{Float64}(undef,n)        # keep references to avoid garbage collection
#     ret=DSPDoubleSplitComplex(retr,reti)
#
#
#     ccall(("vDSP_fft_zopD",libacc),Cvoid,
#           (Ptr{Cvoid},DSPDoubleSplitComplex,vDSP_Stride,     DSPDoubleSplitComplex,vDSP_Stride,  vDSP_Length,FFTDirection),
#            plan,     vals,                 SIGNAL_STRIDE,   ret,                  SIGNAL_STRIDE,logn,       FFT_FORWARD)
#
#     retr+im*reti
# End


text:
using JuMP, KNITRO, MathProgComplex

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "Matpower case9 KNITRO" begin
    instancepath = joinpath(Pkg.dir("MathProgComplex"), "test", "instances")
    pb_real, ~ = import_from_dat(joinpath(instancepath,"case9real.dat"))
    mysolver = KnitroSolver(KTR_PARAM_OUTLEV=3,
                          KTR_PARAM_MAXIT=600,
                          KTR_PARAM_SCALE=0,
                          KTR_PARAM_FEASTOL=1.0,
                          KTR_PARAM_OPTTOL=1.0,
                          KTR_PARAM_FEASTOLABS=1e-6,
                          KTR_PARAM_OPTTOLABS=1e-3,
                          KTR_PARAM_BAR_INITPT=2,
                          KTR_PARAM_PRESOLVE=0,
                          KTR_PARAM_HONORBNDS=0)

    m, variables_jump = get_JuMP_cartesian_model(pb_real, mysolver)
    solve(m)
    @test getobjectivevalue(m) ≈ 1.45883471040144e+003 atol=1e-6
end


text:
using ArgParse,FileIO,JLD,KUparser,Knet
macro msg(_x) :(if logging>0; join(STDERR,[Dates.format(now(),"HH:MM:SS"), $_x,'\n'],' '); end) end
macro log(_x) :(@msg($(string(_x))); $(esc(_x))) end

include("conll17.jl")

# for i in `cut -f2 languages.tsv | head -2`; do echo $i; julia gparser-mlp-static.jl --train $i-ud-train.conllu $i-ud-dev.conllu --vectors $i.vectors.xz --parse $i-ud-dev.txt --udpipe $i.udpipe > foo; ./eval.py $i-ud-dev.conllu foo; done

# Some example data
grctrn = "/mnt/ai/data/nlp/conll17/ud-treebanks-conll2017/UD_Ancient_Greek/grc-ud-train.conllu"
grcdev = "/mnt/ai/data/nlp/conll17/ud-treebanks-conll2017/UD_Ancient_Greek/grc-ud-dev.conllu"
grctxt = "/mnt/ai/data/nlp/conll17/ud-treebanks-conll2017/UD_Ancient_Greek/grc-ud-dev.txt"
grcvec = "/mnt/ai/data/nlp/conll17/word-embeddings-conll17/Ancient_Greek/grc.vectors.xz"
grcudp = "/mnt/ai/data/nlp/conll17/UDPipe/udpipe-ud-2.0-conll17-170315/models/ancient_greek-ud-2.0-conll17-170315.udpipe"

function main(args="")
    isa(args, AbstractString) && (args=split(args))
    s = ArgParseSettings()
    s.description="Using greedy parser with oracle (static) trained MLP model"
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        # Parsing options
        ("--parse"; help="file in plain text format to be parsed")
        ("--udpipe"; help="UDpipe model file")
        ("--load"; help="load model from jld file")
        ("--ptype"; default="ArcEagerR1"; help="Parser type")
        ("--flist"; default="Flist.zn11pv"; help="Feature list to use")

        # Training options
        ("--train"; nargs='*'; help="training file(s) in conllu format")
        ("--vectors"; help="word vectors")
        ("--save"; help="save model to jld file")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"); help="array type: Array for cpu, KnetArray for gpu")
        ("--otype"; default="Adam()"; help="Optimization algorithm and parameters.")
        ("--hidden"; nargs='*'; arg_type=Int; help="sizes of hidden layers, e.g. --hidden 128 64 for a net with two hidden layers")
        ("--epochs"; arg_type=Int; default=10; help="number of epochs for training")
        ("--batchsize"; arg_type=Int; default=100; help="minibatch size")
        ("--dropout"; arg_type=Float64; default=0.0; help="dropout probability")
        ("--l1"; arg_type=Float64; default=0.0; help="L1 regularization")
        ("--l2"; arg_type=Float64; default=0.0; help="L2 regularization")
        ("--fast"; action=:store_true; help="skip loss calculation for faster run")
        #TODO ("--gcheck"; arg_type=Int; default=0; help="check N random gradients per parameter")
        ("--seed"; arg_type=Int; default=-1; help="random number seed: use a nonnegative int for repeatable results")
        ("--logging"; arg_type=Int; default=1; help="Logging level: 0 to turn off")
    end

    # Process args
    o = parse_args(args, s; as_symbols=true)
    logging = o[:logging]
    @msg(s.description)
    @msg(string("opts=",[(k,v) for (k,v) in o]...))
    if o[:seed] > 0; setseed(o[:seed]); end
    atype = eval(parse(o[:atype]))
    ptype = eval(parse(o[:ptype]))
    flist = eval(parse(o[:flist]))

    # Load model
    if o[:load] != nothing
        @log a = load(o[:load])
        model = a["model"]
        optim = a["optim"]
        vocab = a["vocab"]
    else
        model = optim = nothing
        vocab = Vocab(); vocab.postags=UPOSTAG; vocab.deprels=UDEPREL
    end

    # Load raw text
    if o[:parse] != nothing
        # we need to test with predicted toks and tags from txt data.
        parse_conllu = tempname()
        @log run(pipeline(`udpipe --tokenize --tag $(o[:udpipe]) $(o[:parse])`, parse_conllu))
        @log parse_sentences = readconllu(parse_conllu, vocab)
    end

    # Load training data
    if !isempty(o[:train])
        readc(f)=readconllu(f,vocab)
        @log train_corpora = map(readc, o[:train])
    end

    # Load vectors: this has to be done after all conllu files have been loaded
    if o[:vectors] != nothing
        @log readvectors(o[:vectors], vocab)
    end

    # Fix words without vectors, this may happen if the model loaded does not have all the words in the corpora read
    if size(vocab.wvecs,2) < length(vocab.words)
        z = fill!(similar(vocab.wvecs, size(vocab.wvecs,1), length(vocab.words) - size(vocab.wvecs,2)), 0)
        @log vocab.wvecs = hcat(vocab.wvecs, z)
    end

    # Train model
    if !isempty(o[:train])
        data = map(train_corpora) do corpus
            @log (p,x,y) = oparse(ptype, corpus, flist)
            @log minibatch(x, y, o[:batchsize]; atype=atype)
        end
        if model == nothing
            xtrn,ytrn = data[1][1]
            model = initmodel(size(xtrn,1), o[:hidden]..., size(ytrn,1); atype=atype)
            optim = initoptim(model, o[:otype])
        end
        if o[:fast]
            report = (ep->nothing)
        else
            report = (ep->(@msg((ep,map(d->accuracy(model,d,mlp),data)...))))
            @msg((:epoch,:trnacc,:devacc))
        end
        report(0)
        for epoch=1:o[:epochs]
            for (x,y) in data[1]
                grads = softgrad(model,x,y,mlp;l1=o[:l1],l2=o[:l2],pdrop=o[:dropout])
                update!(model,grads,optim)
            end
            report(epoch)
        end
    end

    # Save model
    if o[:save] != nothing
        @log save(File(format"JLD",o[:save]), "model", model, "optim", optim, "vocab", vocab)
    end

    # Parse
    if o[:parse] != nothing
        pred(x,y)=copy!(y, mlp(model, atype(x)))
        @log parses = gparse(ptype, parse_sentences, flist, pred; nbatch=o[:batchsize])
        @log writeconllu(parse_sentences, parses, parse_conllu)
        rm(parse_conllu)
    end
end

# function savedata()
#     @date c = readconllu(corpus, wordvecs)
#     @date (p,x,y) = oparse(pt, c, ft)
#     @date save("foo.en.jld", "x", x, "y", y, "c", c, "p", p)
# end

function minibatch(x, y, batch; atype=Array{Float32})
    nx = size(x,2)
    data = Any[]
    for i=1:batch:nx
        j=min(nx,i+batch-1)
        push!(data, (atype(view(x,:,i:j)),atype(view(y,:,i:j))))
    end
    return data
end

function mlp(w,x;pdrop=0)
    for i=1:2:length(w)-2
        x = relu(w[i]*x .+ w[i+1])
        x = dropout(x, pdrop)
    end
    return w[end-1]*x .+ w[end]
end

# # For use with gparser until we change the interface:
# function predict(w,x,y)
#     x = KnetArray(x)
#     for i=1:2:length(w)-2
#         x = relu(w[i]*x .+ w[i+1])
#     end
#     # ccall(("cudaDeviceSynchronize","libcudart"),UInt32,()) # uncomment for profiling
#     copy!(y, w[end-1]*x .+ w[end])
# end

using Knet: @cuda, Cptr

# target may be a subarray, TODO: handle in Knet
function Base.copy!{T}(dest::SubArray{T}, src::KnetArray{T})
    if length(dest) < length(src); throw(BoundsError()); end
    @cuda(cudart,cudaMemcpy,(Cptr,Cptr,Csize_t,UInt32),
          pointer(dest), pointer(src), length(src)*sizeof(T), 2)
    return dest
end

# Calculate cross entropy loss of a model with weights w for one minibatch (x,p)
# Use non-zero l1 or l2 for regularization
function softloss(w,x,p,model;l1=0,l2=0,o...)
    y = model(w,x;o...)
    logphat = logp(y)
    J = -sum(p .* logphat) / size(x,2)  # dividing by number of instances for per-instance average
    if l1 != 0; J += l1 * sum(sumabs(wi)  for wi in w[1:2:end]); end
    if l2 != 0; J += l2 * sum(sumabs2(wi) for wi in w[1:2:end]); end
    return J
end

softgrad = grad(softloss)

function avgloss(w,data,model) # average loss for the whole dataset
    sum = cnt = 0
    for (x,y) in data
        sum += softloss(w,x,y,model)
        cnt += 1
    end
    return sum/cnt
end

function accuracy(w,data,model)
    corr = cnt = 0
    for (x,y) in data
        ypred = mlp(w,x)
        corr += sum(y .* (ypred .== maximum(ypred,1)))
        cnt += size(y,2)
    end
    return corr/cnt
end

function initmodel(d...; atype=Array{Float32})
    init(d...)=atype(xavier(Float32,d...))
    bias(d...)=atype(zeros(Float32,d...))
    w = Any[]
    for i=1:length(d)-1
        push!(w, init(d[i+1],d[i]))
        push!(w, bias(d[i+1],1))
    end
    return w
end

# This should work for any combination of tuple/array/dict
initoptim{T<:Number}(::KnetArray{T},otype)=eval(parse(otype))
initoptim{T<:Number}(::Array{T},otype)=eval(parse(otype))
initoptim(a::Associative,otype)=Dict(k=>initoptim(v,otype) for (k,v) in a)
initoptim(a,otype)=map(x->initoptim(x,otype), a)

# function train(dtrn,dtst;
#                hidden = [ 128 ],
#                optimizer = Adam,
#                epochs = 10,
#                l1 = 0, l2 = 0, pdrop = 0
#                )
#     (x,y)=dtrn[1]
#     w = initmodel(size(x,1), hidden..., size(y,1))
#     p = oparams(w, optimizer)
#     report(ep)=println((ep,accuracy(w,dtrn,mlp),accuracy(w,dtst,mlp)))
#     report(0)
#     for epoch=1:epochs
#         for (x,y) in dtrn
#             g = softgrad(w,x,y,mlp;l1=l1,l2=l2,pdrop=pdrop)
#             update!(w,g,p)
#         end
#         report(epoch)
#     end
#     return w
# end

# julia> dtrn,dtst = loaddata()
# julia> train(dtrn,dtst)
# (0,0.008304813f0,0.01f0)
# (1,0.9219519f0,0.9071f0)
# (2,0.9373984f0,0.9188f0)
# (3,0.9448262f0,0.9239f0)
# (4,0.9502219f0,0.923f0)
# (5,0.95345455f0,0.9229f0)
# (6,0.95614976f0,0.923f0)
# (7,0.95832354f0,0.9228f0)
# (8,0.9598797f0,0.9223f0)
# (9,0.962f0,0.9258f0)
# (10,0.9629759f0,0.9233f0)

# julia> train(dtrn,dtst;hidden=1000)
# (0,0.0073877005f0,0.0074f0)
# (1,0.9364385f0,0.9173f0)
# (2,0.9535187f0,0.9284f0)
# (3,0.9634011f0,0.9316f0)
# (4,0.97001874f0,0.9315f0)
# (5,0.97462296f0,0.9295f0)
# (6,0.9792781f0,0.9307f0)
# (7,0.982254f0,0.9316f0)
# (8,0.98386633f0,0.9325f0)
# (9,0.98568714f0,0.9291f0)
# (10,0.98813367f0,0.9312f0)

# julia> train(dtrn,dtst;hidden=10000)
# (0,0.0077700536f0,0.0077f0)
# (1,0.9415294f0,0.9223f0)
# (2,0.95814705f0,0.9313f0)
# (3,0.9689037f0,0.9358f0)
# (4,0.97555614f0,0.932f0)
# (5,0.98095185f0,0.9337f0)
# (6,0.9839626f0,0.934f0)
# (7,0.98727006f0,0.9358f0)
# (8,0.9889813f0,0.9356f0)
# (9,0.9904198f0,0.9358f0)
# (10,0.99142516f0,0.9338f0)

if basename(PROGRAM_FILE)==basename(@__FILE__); main(ARGS); end


text:
# NOTE: s₃ = image(+, s₁, s₂) satisfies indices(s₃) == indices(s₁) ∪ indices(s₂)

function +̄ end
function -̄ end
function *̄ end
function ^̄ end

# fallback methods

image(::typeof(-), s₁::VectorSpace, s₂::VectorSpace) = image(+, s₁, s₂)
image(::typeof(-̄), s₁::VectorSpace, s₂::VectorSpace) = image(+̄, s₁, s₂)

# Parameter space

image(::typeof(+), ::ParameterSpace, ::ParameterSpace) = ParameterSpace()
image(::typeof(+̄), ::ParameterSpace, ::ParameterSpace) = ParameterSpace()

# Sequence spaces

for f ∈ (:+, :+̄, :*, :*̄)
    @eval image(::typeof($f), s₁::TensorSpace{<:NTuple{N,BaseSpace}}, s₂::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
        TensorSpace(map((s₁ᵢ, s₂ᵢ) -> image($f, s₁ᵢ, s₂ᵢ), spaces(s₁), spaces(s₂)))
end

function image(::typeof(^), s::SequenceSpace, n::Int)
    n < 0 && return throw(DomainError(n, "^ is only defined for positive integers"))
    n == 0 && return s
    n == 1 && return s
    n == 2 && return image(*, s, s)
    # power by squaring
    t = trailing_zeros(n) + 1
    n >>= t
    while (t -= 1) > 0
        s = image(*, s, s)
    end
    new_s = s
    while n > 0
        t = trailing_zeros(n) + 1
        n >>= t
        while (t -= 1) ≥ 0
            s = image(*, s, s)
        end
        new_s = image(*, new_s, s)
    end
    return new_s
end

function image(::typeof(^̄), s::SequenceSpace, n::Int)
    n < 0 && return throw(DomainError(n, "^̄ is only defined for positive integers"))
    n == 0 && return s
    n == 1 && return s
    n == 2 && return image(*̄, s, s)
    # power by squaring
    t = trailing_zeros(n) + 1
    n >>= t
    while (t -= 1) > 0
        s = image(*̄, s, s)
    end
    new_s = s
    while n > 0
        t = trailing_zeros(n) + 1
        n >>= t
        while (t -= 1) ≥ 0
            s = image(*̄, s, s)
        end
        new_s = image(*̄, new_s, s)
    end
    return new_s
end

# Taylor

image(::typeof(+), s₁::Taylor, s₂::Taylor) = union(s₁, s₂)
image(::typeof(+̄), s₁::Taylor, s₂::Taylor) = intersect(s₁, s₂)

image(::typeof(*), s₁::Taylor, s₂::Taylor) = Taylor(order(s₁) + order(s₂))
image(::typeof(*̄), s₁::Taylor, s₂::Taylor) = intersect(s₁, s₂)

# Fourier

image(::typeof(+), s₁::Fourier, s₂::Fourier) = union(s₁, s₂)
image(::typeof(+̄), s₁::Fourier, s₂::Fourier) = intersect(s₁, s₂)

function image(::typeof(*), s₁::Fourier{T}, s₂::Fourier{S}) where {T,S}
    frequency(s₁) == frequency(s₂) || return throw(DomainError)
    R = promote_type(T, S)
    return Fourier(order(s₁) + order(s₂), convert(R, frequency(s₁)))
end
image(::typeof(*̄), s₁::Fourier, s₂::Fourier) = intersect(s₁, s₂)

# Chebyshev

image(::typeof(+), s₁::Chebyshev, s₂::Chebyshev) = union(s₁, s₂)
image(::typeof(+̄), s₁::Chebyshev, s₂::Chebyshev) = intersect(s₁, s₂)

image(::typeof(*), s₁::Chebyshev, s₂::Chebyshev) = Chebyshev(order(s₁) + order(s₂))
image(::typeof(*̄), s₁::Chebyshev, s₂::Chebyshev) = intersect(s₁, s₂)

# Cartesian spaces

for f ∈ (:+, :+̄)
    @eval begin
        function image(::typeof($f), s₁::CartesianPower, s₂::CartesianPower)
            nb_cartesian_product(s₁) == nb_cartesian_product(s₂) || return throw(DimensionMismatch)
            return CartesianPower(image($f, space(s₁), space(s₂)), nb_cartesian_product(s₁))
        end
        function image(::typeof($f), s₁::CartesianProduct, s₂::CartesianProduct)
            nb_cartesian_product(s₁) == nb_cartesian_product(s₂) || return throw(DimensionMismatch)
            return CartesianProduct(map((s₁ᵢ, s₂ᵢ) -> image($f, s₁ᵢ, s₂ᵢ), spaces(s₁), spaces(s₂)))
        end
        function image(::typeof($f), s₁::CartesianPower, s₂::CartesianProduct)
            nb_cartesian_product(s₁) == nb_cartesian_product(s₂) || return throw(DimensionMismatch)
            return CartesianProduct(map(s₂ᵢ -> image($f, space(s₁), s₂ᵢ), spaces(s₂)))
        end
        function image(::typeof($f), s₁::CartesianProduct, s₂::CartesianPower)
            nb_cartesian_product(s₁) == nb_cartesian_product(s₂) || return throw(DimensionMismatch)
            return CartesianProduct(map(s₁ᵢ -> image($f, s₁ᵢ, space(s₂)), spaces(s₁)))
        end
    end
end


text:
using VariantCallFormat
using Documenter

makedocs(;
    modules=[VariantCallFormat],
    authors="Kenta Sato, Ben J. Ward, Rasmus Henningsson, The BioJulia Organisation and other contributors.",
    repo="https://github.com/rasmushenningsson/VariantCallFormat.jl/blob/{commit}{path}#L{line}",
    sitename="VariantCallFormat.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://rasmushenningsson.github.io/VariantCallFormat.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "VCF and BCF formatted files" => "vcf-bcf.md",
        "Reference" => "reference.md",
    ],
)

deploydocs(;
    repo="github.com/rasmushenningsson/VariantCallFormat.jl",
    devbranch="main",
)