text:
using TransportBasedInference: evaluate

@testset "Test tools for expanded function Nx = 1 " begin

    Nx = 1
    Ne = 500
    X = zeros(Nx, Ne)

    X .= randn(Nx, Ne)

    Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8); CstPhyHermiteBasis(8); CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]
    for b in Blist
        B = MultiBasis(b, Nx)

        idx = reshape([0], (1,1))
        truncidx = idx[1:2:end,:]
        Nψ = 1

        coeff = randn(Nψ)
        f = ExpandedFunction(B, idx, coeff)

        ## Test evaluate_offdiagbasis
        ψoff = evaluate_offdiagbasis(f, X)
        ψofftrunc = evaluate_offdiagbasis(f, X, truncidx)

        @test norm(ψoff[:,1:2:end] - ψofftrunc) < 1e-8

        ψofft = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = MultiFunction(MultiBasis(f.B.B,Nx-1), f.idx[j,1:Nx-1])
            for i=1:Ne
                ψofft[i,j] += fj(X[:,i][1:end-1])
            end
        end

        @test norm(ψofft - ψoff)<1e-8

        ## Test evaluate_diagbasis
        ψdiag = evaluate_diagbasis(f, X)
        ψdiagtrunc = evaluate_diagbasis(f, X, truncidx)

        @test norm(ψdiag[:,1:2:end] - ψdiagtrunc) < 1e-8

        ψdiagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = MultiFunction(MultiBasis(f.B.B,1), f.idx[j,Nx:Nx])
            for i=1:Ne
                ψdiagt[i,j] += fj([X[:,i][end]])
            end
        end

        @test norm(ψdiagt - ψdiag)<1e-8

        ## Test grad_xd_diagbasis
        Gψ_xd_diag = grad_xd_diagbasis(f, X)
        Gψ_xd_diagtrunc = grad_xd_diagbasis(f, X, truncidx)
        @test Gψ_xd_diag[:,1:2:end] == Gψ_xd_diagtrunc

        Gψ_xd_diagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                Gψ_xd_diagt[i,j] += ForwardDiff.derivative(fj, X[:,i][end])
            end
        end

        @test norm(Gψ_xd_diagt - Gψ_xd_diag)<1e-10

        ## Test hess_xd_diagbasis

        Hψ_xd_diag = hess_xd_diagbasis(f, X)
        Hψ_xd_diagtrunc = hess_xd_diagbasis(f, X, truncidx)
        @test Hψ_xd_diag[:,1:2:end] == Hψ_xd_diagtrunc

        Hψ_xd_diagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                Hψ_xd_diagt[i,j] += ForwardDiff.derivative(x->ForwardDiff.derivative(fj, x), X[:,i][end])
            end
        end

        @test norm(Hψ_xd_diagt - Hψ_xd_diag)<1e-10

        ## Test function evaluation

        ψ = evaluate(f, X)
        ψt = zeros(Ne)

        for i=1:Ne
            ψt[i] = f(X[:,i])
        end
        @test norm(ψ - ψt)<1e-10

        ## Test grad_xd

        dψxd = grad_xd(f, X)
        dψxdt = zeros(Ne)

        for i=1:Ne
            dψxdt[i] = ForwardDiff.gradient(f, X[:,i])[end]
        end
        @test norm(dψxd - dψxdt)<1e-10

        ## Test hess_xd

        d2ψxd = hess_xd(f, X)
        d2ψxdt = zeros(Ne)

        for i=1:Ne
            d2ψxdt[i] = ForwardDiff.hessian(f, X[:,i])[end, end]
        end

        @test norm(d2ψxd - d2ψxdt)<1e-10

        ## Test grad_coeff

        dψ_coefft = evaluate_basis(f, X)
        dψ_coeff = grad_coeff(f, X)
        dψ_coeffexpanded = grad_coeff(f, X)

        @test norm(dψ_coeff - dψ_coefft) < 1e-10
        @test norm(dψ_coeffexpanded - dψ_coefft) < 1e-10

        ## Test grad_coeff_grad_xd

        dψ_coeff_xd = grad_coeff_grad_xd(f, X)
        dψ_coeff_xdt = zeros(Ne, Nψ)

        for j=1:Nψ
            foffj = MultiFunction(MultiBasis(f.B.B,Nx-1), f.idx[j,1:end-1])
            fdiagj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                dψ_coeff_xdt[i,j] = foffj(X[:,i][1:end-1])*ForwardDiff.derivative(fdiagj, X[:,i][end])
            end
        end

        @test norm(dψ_coeff_xd - dψ_coeff_xdt)<1e-10
    end
end

@testset "Test tools for expanded function Nx = 1 II " begin

    Nx = 1
    Ne = 500
    X = zeros(Nx, Ne)

    X .= randn(Nx, Ne)

    Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8); CstPhyHermiteBasis(8); CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]
    for b in Blist
        B = MultiBasis(b, Nx)

        idx = reshape([0; 1; 2; 3; 4], (5,1))
        truncidx = idx[1:2:end,:]
        Nψ = 5

        coeff = randn(Nψ)
        f = ExpandedFunction(B, idx, coeff)

        ## Test evaluate_offdiagbasis
        ψoff = evaluate_offdiagbasis(f, X)
        ψofftrunc = evaluate_offdiagbasis(f, X, truncidx)

        @test ψoff[:,1:2:end] == ψofftrunc

        ψofft = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = MultiFunction(MultiBasis(f.B.B,Nx-1), f.idx[j,1:Nx-1])
            for i=1:Ne
                ψofft[i,j] += fj(X[:,i][1:end-1])
            end
        end

        @test norm(ψofft - ψoff)<1e-8

        ## Test evaluate_diagbasis
        ψdiag = evaluate_diagbasis(f, X)
        ψdiagtrunc = evaluate_diagbasis(f, X, truncidx)

        @test ψdiag[:,1:2:end] == ψdiagtrunc

        ψdiagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = MultiFunction(MultiBasis(f.B.B,1), f.idx[j,Nx:Nx])
            for i=1:Ne
                ψdiagt[i,j] += fj([X[:,i][end]])
            end
        end

        @test norm(ψdiagt - ψdiag)<1e-8

        ## Test grad_xd_diagbasis
        Gψ_xd_diag = grad_xd_diagbasis(f, X)
        Gψ_xd_diagtrunc = grad_xd_diagbasis(f, X, truncidx)
        @test Gψ_xd_diag[:,1:2:end] == Gψ_xd_diagtrunc

        Gψ_xd_diagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                Gψ_xd_diagt[i,j] += ForwardDiff.derivative(fj, X[:,i][end])
            end
        end

        @test norm(Gψ_xd_diagt - Gψ_xd_diag)<1e-10

        ## Test hess_xd_diagbasis

        Hψ_xd_diag = hess_xd_diagbasis(f, X)
        Hψ_xd_diagtrunc = hess_xd_diagbasis(f, X, truncidx)
        @test Hψ_xd_diag[:,1:2:end] == Hψ_xd_diagtrunc

        Hψ_xd_diagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                Hψ_xd_diagt[i,j] += ForwardDiff.derivative(x->ForwardDiff.derivative(fj, x), X[:,i][end])
            end
        end

        @test norm(Hψ_xd_diagt - Hψ_xd_diag)<1e-10

        ## Test function evaluation

        ψ = evaluate(f, X)
        ψt = zeros(Ne)

        for i=1:Ne
            ψt[i] = f(X[:,i])
        end
        @test norm(ψ - ψt)<1e-10

        ## Test grad_xd

        dψxd = grad_xd(f, X)
        dψxdt = zeros(Ne)

        for i=1:Ne
            dψxdt[i] = ForwardDiff.gradient(f, X[:,i])[end]
        end
        @test norm(dψxd - dψxdt)<1e-10

        ## Test hess_xd

        d2ψxd = hess_xd(f, X)
        d2ψxdt = zeros(Ne)

        for i=1:Ne
            d2ψxdt[i] = ForwardDiff.hessian(f, X[:,i])[end, end]
        end

        @test norm(d2ψxd - d2ψxdt)<1e-10

        ## Test grad_coeff

        dψ_coefft = evaluate_basis(f, X)
        dψ_coeff = grad_coeff(f, X)
        dψ_coeffexpanded = grad_coeff(f, X)

        @test norm(dψ_coeff - dψ_coefft) < 1e-10
        @test norm(dψ_coeffexpanded - dψ_coefft) < 1e-10

        ## Test grad_coeff_grad_xd

        dψ_coeff_xd = grad_coeff_grad_xd(f, X)
        dψ_coeff_xdt = zeros(Ne, Nψ)

        for j=1:Nψ
            foffj = MultiFunction(MultiBasis(f.B.B,Nx-1), f.idx[j,1:end-1])
            fdiagj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                dψ_coeff_xdt[i,j] = foffj(X[:,i][1:end-1])*ForwardDiff.derivative(fdiagj, X[:,i][end])
            end
        end

        @test norm(dψ_coeff_xd - dψ_coeff_xdt)<1e-10
    end

end

@testset "Test tools for expanded function Nx = 2" begin

    Nx = 2
    Ne = 500
    X = zeros(Nx, Ne)

    X .= randn(Nx, Ne)

    Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8); CstPhyHermiteBasis(8); CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]
    for b in Blist
        B = MultiBasis(b, Nx)

        idx = [0 0; 0 1; 1 0; 1 1; 1 2; 2 1]
        truncidx = idx[1:2:end,:]
        Nψ = 6

        coeff = randn(Nψ)
        f = ExpandedFunction(B, idx, coeff)

        ## Test evaluate_offdiagbasis
        ψoff = evaluate_offdiagbasis(f, X)
        ψofftrunc = evaluate_offdiagbasis(f, X, truncidx)

        @test norm(ψoff[:,1:2:end] - ψofftrunc) < 1e-8

        ψofft = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = MultiFunction(MultiBasis(f.B.B,Nx-1), f.idx[j,1:Nx-1])
            for i=1:Ne
                ψofft[i,j] += fj(X[:,i][1:end-1])
            end
        end

        @test norm(ψofft - ψoff)<1e-8

        ## Test evaluate_diagbasis
        ψdiag = evaluate_diagbasis(f, X)
        ψdiagtrunc = evaluate_diagbasis(f, X, truncidx)

        @test norm(ψdiag[:,1:2:end] - ψdiagtrunc) < 1e-8

        ψdiagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = MultiFunction(MultiBasis(f.B.B,1), f.idx[j,Nx:Nx])
            for i=1:Ne
                ψdiagt[i,j] += fj([X[:,i][end]])
            end
        end

        @test norm(ψdiagt - ψdiag)<1e-8

        ## Test grad_xd_diagbasis
        Gψ_xd_diag = grad_xd_diagbasis(f, X)
        Gψ_xd_diagtrunc = grad_xd_diagbasis(f, X, truncidx)
        @test Gψ_xd_diag[:,1:2:end] == Gψ_xd_diagtrunc

        Gψ_xd_diagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                Gψ_xd_diagt[i,j] += ForwardDiff.derivative(fj, X[:,i][end])
            end
        end

        @test norm(Gψ_xd_diagt - Gψ_xd_diag)<1e-10

        ## Test hess_xd_diagbasis

        Hψ_xd_diag = hess_xd_diagbasis(f, X)
        Hψ_xd_diagtrunc = hess_xd_diagbasis(f, X, truncidx)
        @test Hψ_xd_diag[:,1:2:end] == Hψ_xd_diagtrunc

        Hψ_xd_diagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                Hψ_xd_diagt[i,j] += ForwardDiff.derivative(x->ForwardDiff.derivative(fj, x), X[:,i][end])
            end
        end

        @test norm(Hψ_xd_diagt - Hψ_xd_diag)<1e-10

        ## Test function evaluation

        ψ = evaluate(f, X)
        ψt = zeros(Ne)

        for i=1:Ne
            ψt[i] = f(X[:,i])
        end
        @test norm(ψ - ψt)<1e-10

        ## Test grad_xd

        dψxd = grad_xd(f, X)
        dψxdt = zeros(Ne)

        for i=1:Ne
            dψxdt[i] = ForwardDiff.gradient(f, X[:,i])[end]
        end
        @test norm(dψxd - dψxdt)<1e-10

        ## Test hess_xd

        d2ψxd = hess_xd(f, X)
        d2ψxdt = zeros(Ne)

        for i=1:Ne
            d2ψxdt[i] = ForwardDiff.hessian(f, X[:,i])[end, end]
        end

        @test norm(d2ψxd - d2ψxdt)<1e-10

        ## Test grad_coeff

        dψ_coefft = evaluate_basis(f, X)
        dψ_coeff = grad_coeff(f, X)
        dψ_coeffexpanded = grad_coeff(f, X)

        @test norm(dψ_coeff - dψ_coefft) < 1e-10
        @test norm(dψ_coeffexpanded - dψ_coefft) < 1e-10

        ## Test grad_coeff_grad_xd

        dψ_coeff_xd = grad_coeff_grad_xd(f, X)
        dψ_coeff_xdt = zeros(Ne, Nψ)

        for j=1:Nψ
            foffj = MultiFunction(MultiBasis(f.B.B,Nx-1), f.idx[j,1:end-1])
            fdiagj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                dψ_coeff_xdt[i,j] = foffj(X[:,i][1:end-1])*ForwardDiff.derivative(fdiagj, X[:,i][end])
            end
        end

        @test norm(dψ_coeff_xd - dψ_coeff_xdt)<1e-10
    end

end

@testset "Test tools for expanded function Nx = 3" begin

    Nx = 3
    Ne = 500
    X = zeros(Nx, Ne)

    X .= randn(Nx, Ne)

    Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8); CstPhyHermiteBasis(8); CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]
    for b in Blist
        B = MultiBasis(b, Nx)

        idx = [0 0 0 ;0  0 1; 0 1 0; 1 0 0;1 1 0; 0 1 1; 1 0 1; 1 1 1; 1 2 0; 2 1 0]
        truncidx = idx[1:2:end,:]
        Nψ = 10

        coeff = randn(Nψ)
        f = ExpandedFunction(B, idx, coeff)

        ## Test evaluate_offdiagbasis
        ψoff = evaluate_offdiagbasis(f, X)
        ψofftrunc = evaluate_offdiagbasis(f, X, truncidx)

        @test norm(ψoff[:,1:2:end] - ψofftrunc) < 1e-8

        ψofft = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = MultiFunction(MultiBasis(f.B.B,Nx-1), f.idx[j,1:Nx-1])
            for i=1:Ne
                ψofft[i,j] += fj(X[:,i][1:end-1])
            end
        end

        @test norm(ψofft - ψoff)<1e-8

        ## Test evaluate_diagbasis
        ψdiag = evaluate_diagbasis(f, X)
        ψdiagtrunc = evaluate_diagbasis(f, X, truncidx)

        @test norm(ψdiag[:,1:2:end] - ψdiagtrunc) < 1e-8

        ψdiagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = MultiFunction(MultiBasis(f.B.B,1), f.idx[j,Nx:Nx])
            for i=1:Ne
                ψdiagt[i,j] += fj([X[:,i][end]])
            end
        end

        @test norm(ψdiagt - ψdiag)<1e-8

        ## Test grad_xd_diagbasis
        Gψ_xd_diag = grad_xd_diagbasis(f, X)
        Gψ_xd_diagtrunc = grad_xd_diagbasis(f, X, truncidx)
        @test Gψ_xd_diag[:,1:2:end] == Gψ_xd_diagtrunc

        Gψ_xd_diagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                Gψ_xd_diagt[i,j] += ForwardDiff.derivative(fj, X[:,i][end])
            end
        end

        @test norm(Gψ_xd_diagt - Gψ_xd_diag)<1e-10

        ## Test hess_xd_diagbasis

        Hψ_xd_diag = hess_xd_diagbasis(f, X)
        Hψ_xd_diagtrunc = hess_xd_diagbasis(f, X, truncidx)
        @test Hψ_xd_diag[:,1:2:end] == Hψ_xd_diagtrunc

        Hψ_xd_diagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                Hψ_xd_diagt[i,j] += ForwardDiff.derivative(x->ForwardDiff.derivative(fj, x), X[:,i][end])
            end
        end

        @test norm(Hψ_xd_diagt - Hψ_xd_diag)<1e-10

        ## Test function evaluation

        ψ = evaluate(f, X)
        ψt = zeros(Ne)

        for i=1:Ne
            ψt[i] = f(X[:,i])
        end
        @test norm(ψ - ψt)<1e-10

        ## Test grad_xd

        dψxd = grad_xd(f, X)
        dψxdt = zeros(Ne)

        for i=1:Ne
            dψxdt[i] = ForwardDiff.gradient(f, X[:,i])[end]
        end
        @test norm(dψxd - dψxdt)<1e-10

        ## Test hess_xd

        d2ψxd = hess_xd(f, X)
        d2ψxdt = zeros(Ne)

        for i=1:Ne
            d2ψxdt[i] = ForwardDiff.hessian(f, X[:,i])[end, end]
        end

        @test norm(d2ψxd - d2ψxdt)<1e-10

        ## Test grad_coeff

        dψ_coefft = evaluate_basis(f, X)
        dψ_coeff = grad_coeff(f, X)
        dψ_coeffexpanded = grad_coeff(f, X)

        @test norm(dψ_coeff - dψ_coefft) < 1e-10
        @test norm(dψ_coeffexpanded - dψ_coefft) < 1e-10

        ## Test grad_coeff_grad_xd

        dψ_coeff_xd = grad_coeff_grad_xd(f, X)
        dψ_coeff_xdt = zeros(Ne, Nψ)

        for j=1:Nψ
            foffj = MultiFunction(MultiBasis(f.B.B,Nx-1), f.idx[j,1:end-1])
            fdiagj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                dψ_coeff_xdt[i,j] = foffj(X[:,i][1:end-1])*ForwardDiff.derivative(fdiagj, X[:,i][end])
            end
        end

        @test norm(dψ_coeff_xd - dψ_coeff_xdt)<1e-10
    end
end

@testset "Test tools for expanded function Nx = 3" begin

    Nx = 3
    Ne = 500

    X = randn(Nx, Ne)

    Blist = [ProHermiteBasis(8); PhyHermiteBasis(8); CstProHermiteBasis(8); CstPhyHermiteBasis(8); CstLinProHermiteBasis(8); CstLinPhyHermiteBasis(8)]
    for b in Blist
        B = MultiBasis(b, Nx)

        idx = [0 0 0]
        truncidx = idx[1:2:end,:]
        Nψ = 1

        coeff = randn(Nψ)
        f = ExpandedFunction(B, idx, coeff)

        ## Test evaluate_offdiagbasis
        ψoff = evaluate_offdiagbasis(f, X)
        ψofftrunc = evaluate_offdiagbasis(f, X, truncidx)

        @test norm(ψoff[:,1:2:end] - ψofftrunc) < 1e-8

        ψofft = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = MultiFunction(MultiBasis(f.B.B,Nx-1), f.idx[j,1:Nx-1])
            for i=1:Ne
                ψofft[i,j] += fj(X[:,i][1:end-1])
            end
        end

        @test norm(ψofft - ψoff)<1e-8

        ## Test evaluate_diagbasis
        ψdiag = evaluate_diagbasis(f, X)
        ψdiagtrunc = evaluate_diagbasis(f, X, truncidx)

        @test norm(ψdiag[:,1:2:end] - ψdiagtrunc) < 1e-8

        ψdiagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = MultiFunction(MultiBasis(f.B.B,1), f.idx[j,Nx:Nx])
            for i=1:Ne
                ψdiagt[i,j] += fj([X[:,i][end]])
            end
        end

        @test norm(ψdiagt - ψdiag)<1e-8

        ## Test grad_xd_diagbasis
        Gψ_xd_diag = grad_xd_diagbasis(f, X)
        Gψ_xd_diagtrunc = grad_xd_diagbasis(f, X, truncidx)
        @test Gψ_xd_diag[:,1:2:end] == Gψ_xd_diagtrunc

        Gψ_xd_diagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                Gψ_xd_diagt[i,j] += ForwardDiff.derivative(fj, X[:,i][end])
            end
        end

        @test norm(Gψ_xd_diagt - Gψ_xd_diag)<1e-10

        ## Test hess_xd_diagbasis

        Hψ_xd_diag = hess_xd_diagbasis(f, X)
        Hψ_xd_diagtrunc = hess_xd_diagbasis(f, X, truncidx)
        @test Hψ_xd_diag[:,1:2:end] == Hψ_xd_diagtrunc

        Hψ_xd_diagt = zeros(Ne, Nψ)

        for j=1:Nψ
            fj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                Hψ_xd_diagt[i,j] += ForwardDiff.derivative(x->ForwardDiff.derivative(fj, x), X[:,i][end])
            end
        end

        @test norm(Hψ_xd_diagt - Hψ_xd_diag)<1e-10

        ## Test function evaluation

        ψ = evaluate(f, X)
        ψt = zeros(Ne)

        for i=1:Ne
            ψt[i] = f(X[:,i])
        end
        @test norm(ψ - ψt)<1e-10

        ## Test grad_xd

        dψxd = grad_xd(f, X)
        dψxdt = zeros(Ne)

        for i=1:Ne
            dψxdt[i] = ForwardDiff.gradient(f, X[:,i])[end]
        end
        @test norm(dψxd - dψxdt)<1e-10

        ## Test hess_xd

        d2ψxd = hess_xd(f, X)
        d2ψxdt = zeros(Ne)

        for i=1:Ne
            d2ψxdt[i] = ForwardDiff.hessian(f, X[:,i])[end, end]
        end

        @test norm(d2ψxd - d2ψxdt)<1e-10

        ## Test grad_coeff

        dψ_coefft = evaluate_basis(f, X)
        dψ_coeff = grad_coeff(f, X)
        dψ_coeffexpanded = grad_coeff(f, X)

        @test norm(dψ_coeff - dψ_coefft) < 1e-10
        @test norm(dψ_coeffexpanded - dψ_coefft) < 1e-10

        ## Test grad_coeff_grad_xd

        dψ_coeff_xd = grad_coeff_grad_xd(f, X)
        dψ_coeff_xdt = zeros(Ne, Nψ)

        for j=1:Nψ
            foffj = MultiFunction(MultiBasis(f.B.B,Nx-1), f.idx[j,1:end-1])
            fdiagj = f.B.B[f.idx[j,end]+1]
            for i=1:Ne
                dψ_coeff_xdt[i,j] = foffj(X[:,i][1:end-1])*ForwardDiff.derivative(fdiagj, X[:,i][end])
            end
        end

        @test norm(dψ_coeff_xd - dψ_coeff_xdt)<1e-10
    end
end


text:
# ---
# title: 914. X of a Kind in a Deck of Cards
# id: problem914
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Math
# link: <https://leetcode.com/problems/x-of-a-kind-in-a-deck-of-cards/description/>
# hidden: true
# ---
#
# In a deck of cards, each card has an integer written on it.
#
# Return `true` if and only if you can choose `X >= 2` such that it is possible
# to split the entire deck into 1 or more groups of cards, where:
#
#   * Each group has exactly `X` cards.
#   * All the cards in each group have the same integer.
#
#
#
# **Example 1:**
#
#
#
#     Input: deck = [1,2,3,4,4,3,2,1]
#     Output: true
#     **Explanation** : Possible partition [1,1],[2,2],[3,3],[4,4].
#
#
# **Example 2:**
#
#
#
#     Input: deck = [1,1,1,2,2,2,3,3]
#     Output: false´
#     **Explanation** : No possible partition.
#
#
# **Example 3:**
#
#
#
#     Input: deck = [1]
#     Output: false
#     **Explanation** : No possible partition.
#
#
# **Example 4:**
#
#
#
#     Input: deck = [1,1]
#     Output: true
#     **Explanation** : Possible partition [1,1].
#
#
# **Example 5:**
#
#
#
#     Input: deck = [1,1,2,2,2,2]
#     Output: true
#     **Explanation** : Possible partition [1,1],[2,2],[2,2].
#
#
#
#
# **Constraints:**
#
#   * `1 <= deck.length <= 10^4`
#   * `0 <= deck[i] < 10^4`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
"""
    function addParameters_!(p::Program_,platform::String)

Generate the variables of the model and declare them.
"""
function addParameters_!(p::Program_,platform::String)

    #Parameter declare###########################################################################

    append!(p.declareVar.args,
        (quote
            simulationBox = Array(com.simulationBox)
        end).args
    )

    if length(p.agent.declaredSymbols["Local"])>0
        append!(p.declareVar.args,
            (quote
                localV = Array($FLOAT.([com.local_;Base.zeros(nMax-N,$(length(p.agent.declaredSymbols["Local"])))]))
            end).args
        )

        push!(p.args,:localV)

        if !isempty(p.update["Local"])
            append!(p.declareVar.args,
                (quote
                    localVCopy = zeros($FLOAT,size(localV)[1],$(length(keys(p.update["Local"]))))
                end).args
            )

            push!(p.args,:localVCopy)
        end
    end

    if length(p.agent.declaredSymbols["Identity"])>0
        append!(p.declareVar.args,
            (quote
                identityV = Array($INT.([com.identity_;Base.zeros(Int,nMax-N,$(length(p.agent.declaredSymbols["Identity"])))]))
            end).args
        )

        push!(p.args,:identityV)

        if !isempty(p.update["Identity"])
            append!(p.declareVar.args,
                (quote
                    identityVCopy = zeros($INT,size(identityV)[1],$(length(keys(p.update["Identity"]))))
                end).args
            )

            push!(p.args,:identityVCopy)
        end
    end

    if length(p.agent.declaredSymbols["Global"])>0
        append!(p.declareVar.args,
            (quote
                globalV = Array($FLOAT.(com.global_))
            end).args
        )

        push!(p.args,:globalV)
        if !isempty(p.update["Global"]) && !isempty([ v for v in keys(p.update["Global"]) if v in p.agent.declaredSymbols["Global"] ])
            append!(p.declareVar.args,
                (quote
                    globalVCopy = zeros($FLOAT,$(length(p.update["Global"])))
                end).args
            )

            push!(p.args,:globalVCopy)
        end
    end

    if length(p.agent.declaredSymbols["GlobalArray"])>0
        for (j,i) in enumerate(p.agent.declaredSymbols["GlobalArray"])
            append!(p.declareVar.args,
            (quote
                $(Meta.parse(string(i))) = Array(copy(com.globalArray_[$j]))
            end).args
            )

            push!(p.args,Meta.parse(string(i)))
        end
        if !isempty(p.update["GlobalArray"])
            for (j,i) in enumerate(keys(p.update["GlobalArray"]))
                append!(p.declareVar.args,
                (quote
                    $(Meta.parse(string(i,GLOBALARRAYCOPY))) = Array(copy(com.globalArray_[$j]))
                end).args
                )

                push!(p.args,Meta.parse(string(i,GLOBALARRAYCOPY)))
            end
        end
    end

    if length(p.agent.declaredSymbols["Medium"])>0
        if p.agent.dims >= 1
            push!(p.declareVar.args, :(Nx_ = com.mediumN[1]+2))
            push!(p.declareVar.args, :(dxₘ_ = (com.simulationBox[1,2]-com.simulationBox[1,1])/(Nx_-2)))

            push!(p.args,:Nx_)
            push!(p.args,:dxₘ_)
        end
        if p.agent.dims >= 2
            push!(p.declareVar.args, :(Ny_ = com.mediumN[2]+2))
            push!(p.declareVar.args, :(dyₘ_ = (com.simulationBox[2,2]-com.simulationBox[2,1])/(Ny_-2)))

            push!(p.args,:Ny_)
            push!(p.args,:dyₘ_)
        end
        if p.agent.dims >= 3
            push!(p.declareVar.args, :(Nz_ = com.mediumN[3]+2))
            push!(p.declareVar.args, :(dzₘ_ = (com.simulationBox[3,2]-com.simulationBox[3,1])/(Nz_-2)))

            push!(p.args,:Nz_)
            push!(p.args,:dzₘ_)
        end
        if p.agent.dims == 1
            push!(p.declareVar.args, :(mediumV = Base.zeros($FLOAT,size(com.medium_).+(2,0))))
            push!(p.declareVar.args, :(mediumV[2:end-1,:] = com.medium_))
        elseif p.agent.dims == 2
            push!(p.declareVar.args, :(mediumV = Base.zeros($FLOAT,size(com.medium_).+(2,2,0))))
            push!(p.declareVar.args, :(mediumV[2:end-1,2:end-1,:] = com.medium_))
        elseif p.agent.dims == 3
            push!(p.declareVar.args, :(mediumV = Base.zeros($FLOAT,size(com.medium_).+(2,2,2,0))))
            push!(p.declareVar.args, :(mediumV[2:end-1,2:end-1,2:end-1,:] = com.medium_))
        end
        push!(p.declareVar.args, :(mediumV = Array(mediumV)))

        push!(p.args,:mediumV)

        if "UpdateMedium" in keys(p.agent.declaredUpdates)
            push!(p.declareVar.args, :(mediumVCopy = copy(mediumV)))
            push!(p.args,:mediumVCopy)
        end

    end

    return nothing
end

"""
    function addUpdate_!(p::Program_,platform::String)

Generate the functions to update all the modified values.
"""
function addUpdate_!(p::Program_,platform::String)

    create = false

    #Check is there is something to update
    for i in keys(p.update)
        if !isempty(p.update[i])
            create = true
        end
    end

    if create
        gen = quote end #General update
        up = quote end
        for i in keys(p.agent.declaredSymbols)
            if  i == "GlobalArray"
                for j in keys(p.update[i])
                    n = Meta.parse(string(j,GLOBALARRAYCOPY))
                    push!(gen.args,:($j.=$n))
                end
            elseif i in ["Local","Identity"]
                var = Meta.parse(string(lowercase(i),"V"))
                varCopy = Meta.parse(string(lowercase(i),"VCopy"))
                for j in keys(p.update[i])
                    pos = findfirst(p.agent.declaredSymbols[i] .== j)
                    posCopy = p.update[i][j]
                    push!(up.args,:($var[ic1_,$pos]=$varCopy[ic1_,$posCopy]))
                end
            elseif i == "Global"
                var = Meta.parse(string(lowercase(i),"V"))
                varCopy = Meta.parse(string(lowercase(i),"VCopy"))

                upb = quote end #Assign updates of the global to the first thread
                for j in keys(p.update[i])
                    pos = findfirst(p.agent.declaredSymbols[i] .== j)
                    posCopy = p.update[i][j]
                    push!(upb.args,:($var[$pos]=$varCopy[$posCopy]))
                end

                push!(up.args,
                    :(begin
                        if ic1_ == 1
                            $(upb)
                        end
                    end))
            elseif i == "Medium"
                if "UpdateMedium" in keys(p.agent.declaredUpdates)
                    push!(gen.args,:(mediumV.=mediumVCopy))
                end
            end

        end

        #Construct global and local update function
        f = simpleFirstLoopWrapInFunction_(platform,:updateLocGlob_!,up)
        push!(p.declareF.args,f)
        #Construct general update function
        push!(p.declareF.args,
            :(begin
                function update_!(ARGS_)
                    $gen
                    @platformAdapt updateLocGlob_!(ARGS_)
                    return
                end
            end)
            )

        #Add to execution cue
        push!(p.execInloop.args,:(update_!(ARGS_)))
    end

    return nothing
end

"""
    function addCopyInitialisation_!(p::Program_,platform::String)

Generate the functions to update all the modified values.
"""
function addCopyInitialisation_!(p::Program_,platform::String)

    create = false

    for i in keys(p.agent.declaredUpdates)
        if !(i in ["UpdateLocalInteraction","UpdateInteraction","EventDivision"]) && !emptyquote_(p.agent.declaredUpdates[i])
            create = true
        end
    end

    if create
        gen = quote end #General update
        up = quote end
        for i in keys(p.agent.declaredSymbols)
            if  i == "GlobalArray"
                for j in keys(p.update[i])
                    n = Meta.parse(string(j,GLOBALARRAYCOPY))
                    push!(gen.args,:($n.=$j))
                end
            elseif i in ["Local","Identity"]
                var = Meta.parse(string(lowercase(i),"V"))
                varCopy = Meta.parse(string(lowercase(i),"VCopy"))
                for j in keys(p.update[i])
                    pos = findfirst(p.agent.declaredSymbols[i] .== j)
                    posCopy = p.update[i][j]
                    push!(up.args,:($varCopy[ic1_,$posCopy]=$var[ic1_,$pos]))
                end
            elseif i == "Global"
                var = Meta.parse(string(lowercase(i),"V"))
                varCopy = Meta.parse(string(lowercase(i),"VCopy"))

                upb = quote end #Assign updates of the global to the first thread
                for j in keys(p.update[i])
                    pos = findfirst(p.agent.declaredSymbols[i] .== j)
                    posCopy = p.update[i][j]
                    push!(upb.args,:($varCopy[$posCopy]=$var[$pos]))
                end

                push!(up.args,
                    :(begin
                        if ic1_ == 1
                            $(upb)
                        end
                    end))
            end

        end

        #Construct global and local update function
        f = simpleFirstLoopWrapInFunction_(platform,:initialiseLocGlobCopy_!,up)
        push!(p.declareF.args,f)
        #Construct general update function
        push!(p.declareF.args,
            :(begin
                function initialiseCopy_!(ARGS_)
                    $gen
                    @platformAdapt initialiseLocGlobCopy_!(ARGS_)
                    return
                end
            end)
            )

        #Add to execution cue
        push!(p.execInit.args,:(initialiseCopy_!(ARGS_)))
    end

    return nothing
end


text:
using AbstractGPs
using ConjugateComputationVI
using Distributions
using FiniteDifferences
using KernelFunctions
using LinearAlgebra
using Random
using StatsFuns
# using TemporalGPs
using Test
using Zygote

using ConjugateComputationVI:
    approx_posterior,
    batch_quadrature,
    canonical_from_expectation,
    canonical_from_natural,
    expectation_from_canonical,
    gaussian_reconstruction_term,
    natural_from_canonical,
    optimise_approx_posterior,
    update_approx_posterior

function generate_synthetic_problem(rng::AbstractRNG)
    f = GP(Matern52Kernel())
    x = collect(range(-2.0, 2.0; length=7))
    Σ = Diagonal(rand(rng, 7) .+ 1e-2)
    y = rand(rng, f(x, Σ))
    return f, x, Σ, y
end

@testset "ConjugateComputationVI.jl" begin
    include("quadrature.jl")
    include("parametrisation_bijections.jl")
    include("approximate_inference.jl")
end


text:
Graph = Dict{Int,Dict{Int,Int}}

neighbors(g::Graph, u::Int) = keys(g[u])

nodes(g::Graph) = keys(g)

hasnode(g::Graph, u::Int) = haskey(g, u)

hasedge(g::Graph, u::Int, v::Int) = hasnode(g, u) && haskey(g[u], v)

function ne(g::Graph)
    m = 0
    for u in keys(g)
        m += length(g[u])
    end
    return div(m, 2)
end

function addnode!(g::Graph, u::Int)
    if !haskey(g, u)
        g[u] = valtype(g)()
    end
    nothing
end

function delnode!(g::Graph, u::Int)
    if haskey(g, u)
        for v in keys(g[u])
            deledge!(g, u, v)
        end
        delete!(g, u)
    end
    nothing
end

function addedge!(g::Graph, u::Int, v::Int)
    if !hasedge(g, u, v)
        addnode!(g, u)
        addnode!(g, v)
        cnt = 1
        for i in keys(g[u])
            if haskey(g[v], i)
                cnt += 1
                g[u][i] += 1
                g[i][u] += 1
                g[v][i] += 1
                g[i][v] += 1
            end
        end
        g[u][v] = cnt
        g[v][u] = cnt
    end
    nothing
end

function deledge!(g::Graph, u::Int, v::Int)
    if haskey(g, u) && haskey(g, v)
        delete!(g[u], v)
        delete!(g[v], u)
        for i in keys(g[u])
            if haskey(g[v], i)
                g[u][i] -= 1
                g[i][u] -= 1
                g[v][i] -= 1
                g[i][v] -= 1
            end
        end
    end
    nothing
end


text:
include("intcode.jl")

function main()
    # Part 1: run the program interactively
    prog = open(f -> parse.(Int64, split(readline(f), ",")), "day25_input.txt")
    input = Channel(1)
    output = Channel(1)
    @sync begin
        # Screen printing loop
        @async begin
            try
                while true
                    print(Char(take!(output)))
                end
            catch InvalidStateException
            end
        end
        # Input loop
        @async begin
            try
                while true
                    command = readline()
                    if command == "q"
                        close(input)
                        close(output)
                        break
                    end
                    [put!(input, c) for c in Int32.(collect("$command\n"))]
                end
            catch InvalidStateException
            end
        end

        # Start the program
        @async begin
            try
                run_intcode!(copy(prog), input, output)
            catch InvalidStateException
            end
        end
    end
    # (answer was 529920)
end
main()

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
module BitmapGraphs

using LightGraphs, SparseArrays

export AbstractBMGraph, BMGraph, invert!, set_diag!, BitRow, SRow, nchunks, mk_msk, SBMGraph
import Base: setindex
#utils for iterating over bit-chunks
@inline _blsr(x) = x & (x-1)

struct LBitRow{chunk_T}
    len::Int
    chunks::chunk_T
end

#future feature
struct BitRow{chunk_T}
    chunks::chunk_T
end
struct SRow{N}
    chunks::NTuple{N, UInt64}
end
const _BitRow = Union{LBitRow, BitRow, SRow}# BitRow

#opt out of simd. Fixme once github.com/JuliaLang/julia/pull/31113 has landed.
Base.SimdLoop.simd_index(v::_BitRow, j::Int64, i) = j
Base.SimdLoop.simd_inner_length(v::_BitRow, j::Int64) = 1
Base.SimdLoop.simd_outer_range(v::_BitRow) = v

nchunks(br::_BitRow) = length(br.chunks)

Base.@propagate_inbounds function Base.in(br::_BitRow, idx)
    i1,i2 = Base.get_chunks_id(idx)
    @boundscheck ((0<i1<len(br)) || return false)
    @inbounds u = br.chunks[i1]
    return iszero(u & (1<< (i2 & 63)))
end

function Base.iterate(br::_BitRow)
    nchunks(br)==0 && return nothing
    return iterate(br, (1, @inbounds br.chunks[1]))
end

#@inline Base.length(br::BitRow) = sum(count_ones, br.chunks)
@inline Base.length(br::LBitRow) = br.len

@inline function Base.iterate(br::_BitRow, s)
    chunks = br.chunks
    i1, c = s
    while c==0
        i1 % UInt >= length(chunks) % UInt && return nothing
        i1 += 1
        @inbounds c = chunks[i1]
    end
    tz = trailing_zeros(c) + 1
    c = _blsr(c)
    return ((i1-1)<<6 + tz, (i1, c))
end

struct BMEdgeIter{chunkT, T}
    ne::Int
    chunks::chunkT
end
@inline Base.eltype(::BMEdgeIter{chunkT, T}) where {chunkT, T} = T
@inline Base.length(ei::BMEdgeIter) = ei.ne

@inline function Base.iterate(br::BMEdgeIter)
    length(br.chunks)>0 || return nothing
    return iterate(br, (1, 1, @inbounds br.chunks[1]))
end

@inline function Base.iterate(br::BMEdgeIter{chunkT, T}, s) where {chunkT, T}
    chunks = br.chunks
    nchunks = size(chunks,1)
    row, col, c = s
    while c==0
        if ((row % UInt) >= nchunks)
            col >= size(chunks, 2) && return nothing
            col += 1
            row,i2 = Base.get_chunks_id(col)
            @inbounds c = chunks[row, col] & (~UInt(0) << (i2 & 63) )
            continue
        end
        row += 1
        @inbounds c = chunks[row, col]
    end
    tz = trailing_zeros(c) + 1
    c = _blsr(c)
    return T((row-1)<<6 + tz, col) ,  (row, col, c)
end

#TODO utils for bitrows:
#intersection, union, xor, setminus, complement, count.
#need to come in broadcast and broadcast! variants.

abstract type AbstractBMGraph <:LightGraphs.AbstractGraph{Int} end

mutable struct BMGraph <: AbstractBMGraph
    adj_chunks::Matrix{UInt64}
    adj_mat::BitMatrix
    degrees::Vector{Int}
    ne::Int
end

function BMGraph(n::Int)
    nc = Base.num_bit_chunks(n)
    adj_mat = falses(nc*64, n)
    adj_chunks = reshape(adj_mat.chunks, nc, n)
    degs = zeros(Int, n)
    return BMGraph(adj_chunks, adj_mat, degs, 0)
end

function BMGraph(g::AbstractGraph)
    vertices(g) == Base.OneTo(nv(g))||throw(ArgumentError("BMGraphs only support `vertices isa Base.OneTo`"))
    res = BMGraph(nv(g))
    for ed in edges(g)
        add_edge!(res, src(ed), dst(ed))
    end
    return res
end
nchunks(g::BMGraph) = size(g.adj_chunks, 1)

function Base.copy(bm::BMGraph)
    adj_mat = copy(bm.adj_mat)
    adj_chunks = reshape(adj_mat.chunks, size(bm.adj_chunks))
    return BMGraph(adj_chunks, adj_mat, copy(bm.degs), bm.ne)
end

#interface:
LightGraphs.edgetype(::AbstractBMGraph) = LightGraphs.SimpleGraphs.SimpleEdge{Int}
Base.eltype(::AbstractBMGraph) = Int
LightGraphs.nv(g::AbstractBMGraph) = size(g.adj_chunks, 2)
LightGraphs.ne(g::BMGraph) = g.ne
LightGraphs.vertices(g::AbstractBMGraph) = Base.OneTo(nv(g))
LightGraphs.edges(g::AbstractBMGraph) = BMEdgeIter{typeof(g.adj_chunks), LightGraphs.SimpleGraphs.SimpleEdge{Int}}(ne(g), g.adj_chunks)
LightGraphs.is_directed(g::AbstractBMGraph) = false
LightGraphs.is_directed(::Type{<:AbstractBMGraph}) = false

LightGraphs.has_vertex(g::BMGraph, v) =  g in vertices(g)

Base.@propagate_inbounds LightGraphs.has_edge(g::BMGraph, e::AbstractEdge) = has_edge(g, src(e), dst(e))
Base.@propagate_inbounds LightGraphs.has_edge(g::BMGraph, e::Tuple{<:Integer, <:Integer}) = has_edge(g, e...)
Base.@propagate_inbounds LightGraphs.has_edge(g::BMGraph, s, d) = g.adj_mat[s, d]
Base.@propagate_inbounds function LightGraphs.outneighbors(g::BMGraph, v)
    vvv = view(g.adj_chunks, :, v)
    return LBitRow(g.degrees[v], vvv)
end
Base.@propagate_inbounds LightGraphs.inneighbors(g::AbstractBMGraph, v...) = outneighbors(g, v...)
Base.@propagate_inbounds LightGraphs.indegree(g::BMGraph, v::Integer) = g.degrees[v]
Base.@propagate_inbounds LightGraphs.outdegree(g::BMGraph, v) = g.degrees[v]
#Base.@propagate_inbounds LightGraphs.degree(g::BMGraph, v) = g.degrees[v] #causes method ambiguity errors

LightGraphs.add_edge!(g::BMGraph, e::AbstractEdge) = add_edge!(g, src(e), dst(e))
LightGraphs.add_edge!(g::BMGraph, e::Tuple{<:Integer, <:Integer}) = add_edge!(g, e...)
Base.@propagate_inbounds function LightGraphs.add_edge!(g::BMGraph, s, d)
    @boundscheck begin
        0 < d <= nv(g) || throw(BoundsError(g, d))
        0 < s <= nv(g) || throw(BoundsError(g, s))
    end
    @inbounds oldv = g.adj_mat[s,d]
    @inbounds g.adj_mat[s,d] = true
    @inbounds g.adj_mat[d,s] = true
    if !oldv
        g.ne += 1
        @inbounds g.degrees[s] += 1
        s != d && @inbounds g.degrees[d] += 1
    end
    return !oldv
end

Base.@propagate_inbounds function LightGraphs.rem_edge!(g::BMGraph, s, d)
    @boundscheck begin
        0 < d <= nv(g) || throw(BoundsError(g, d))
        0 < s <= nv(g) || throw(BoundsError(g, s))
    end
    @inbounds oldv = g.adj_mat[s,d]
    @inbounds g.adj_mat[s,d] = false
    @inbounds g.adj_mat[d,s] = false
    if oldv
        g.ne -= 1
        @inbounds g.degrees[s] -= 1
        s != d && @inbounds g.degrees[d] -= 1
    end
    return oldv
end

"""
invert!(bm::BMGraph,  no_selfloops = true)

Inverts the graph: All formerly present edges are now missing, and vice versa. This includes self-loops.
If `no_selfloops` is set, then all self-loops are removed afterwards.
"""
function invert!(bm::BMGraph, no_selfloops = true)
    nr, nc = size(bm)
    @inbounds for i=1:nc
        bm.degrees[i] = nc - bm.degrees[i]
    end
    bm.ne = sum(bm.degrees)
    @simd for i=1:length(bm.adj_chunks)
        @inbounds bm.adj_chunks[i] = ~bm.adj_chunks[i]
    end
    msk_e = Base._msk_end(nc)
    for j = 1:nr
        @inbounds bm.adj_chunks[j, nr] &= msk_e
    end
    if no_selfloops
        set_diag!(bm, false)
    end

    return bm
end

"""
set_diag!(bm::BMGraph, b::Bool)

Sets the diagonal to b, i.e. adds or removes all self-loops.
"""
function set_diag!(bm::BMGraph, b::Bool)
    if b
        for i=1:nv(bm)
            add_edge!(bm, i, i)
        end
    else
        for i=1:nv(bm)
            rem_edge!(bm, i, i)
        end
    end
    bm
end

"""
induced_subgraph(BMGraph, g::Lightgraphs.AbstractGraph, vertices)

Returns the induced subgraph as a BMGraph. Assumes that input vertices and `neighbors` are sorted.
"""
function LightGraphs.induced_subgraph(::Type{BMGraph}, g::LightGraphs.AbstractGraph, vertices)
    n = length(vertices)
    bm = BMGraph(n)
    for (i,v) in enumerate(vertices)
        vi = 1
        #vv = vertices[1]
        @inbounds for u in neighbors(g, v)
            u <  vertices[vi] && continue
            while vertices[vi] < u
                vi < n || @goto done
                vi += 1
            end
            u < vertices[vi] && continue
            add_edge!(bm, i, vi)
            vi += 1
        end
        @label done
    end
    return bm
end

_true(args...) = true
_false(args...) = false

"""
induced(::Type{BMGraph}, mat::SparseMatrixCSC, vertices, edgepred = _true, bailout = _false)

Returns the induced subgraph, of the graph with adjacency matrix `mat` for which `edgepred(i, j, mat[i,j])` evaluates to `true`.
The matrix is assumed symmetric and `vertices` is assumed sorted.

The predicate `bailout(bm, i)` is called after each constructed row `i`. If it returns `true`, then the partial result is returned.
"""
function LightGraphs.induced_subgraph(::Type{BMGraph}, mat::SparseMatrixCSC, vertices, edgepred = _true, bailout = _false)
    n = length(vertices)
    bm = BMGraph(n)
    n > 1 || return bm
    for (i,v) in enumerate(vertices)
        vi = 1
        ran =  mat.colptr[v]:mat.colptr[v+1]-1
        for ui in ran #neighbors(g, v)
            @inbounds u = mat.rowval[ui]
            @inbounds u <  vertices[vi] && continue

            while @inbounds vertices[vi] < u
                vi < n || @goto done
                vi += 1
            end
            @inbounds u <  vertices[vi] && continue
            edgepred(v, u, @inbounds mat.nzval[ui]) && (@inbounds bm.adj_mat[vi, i] = true)
            vi += 1
            vi < n || @goto done
        end
        @label done
        deg = sum(count_ones, view(bm.adj_chunks, :, i))
        bm.degrees[i] = deg
        bm.ne += deg
        bailout(bm, i) && return bm
    end
    return bm
end

#example algorithm
function LightGraphs.gdistances(g::BMGraph, s)
    res = fill(-1, nv(g))
    nchunks = size(g.adj_chunks, 1)
    visited = zeros(UInt, nchunks)
    todo = zeros(UInt, nchunks)
    nxt = zeros(UInt, nchunks)
    done = false
    dist = 1
    todo .= outneighbors(g, s).chunks
    visited .= outneighbors(g, s).chunks
    while !done
        done = true
        for i in BitRow(todo)
            @simd for j = 1:nchunks
                @inbounds nxt[j] |= g.adj_chunks[j, i] & ~visited[j]
            end
            res[i] = dist
            done = false
        end
        done && break
        visited .|= nxt
        todo .= nxt
        fill!(nxt, 0)
        dist += 1
    end
    res[s] = 0
    res
end

###############
#Static graphs

#Base.in, Base.iterate, are defined above
@inline Base.length(a::SRow) = count_ones(a)
Base.eltype(::SRow) = Int

(Base.:&)(a::SRow{N}, b::SRow{N}) where N = SRow(ntuple(i->(@inbounds a.chunks[i] & b.chunks[i]), N))
(Base.:|)(a::SRow{N}, b::SRow{N}) where N = SRow(ntuple(i->(@inbounds a.chunks[i] | b.chunks[i]), N))
(Base.xor)(a::SRow{N}, b::SRow{N}) where N = SRow(ntuple(i->(@inbounds xor(a.chunks[i], b.chunks[i])), N))
(Base.:~)(a::SRow{N}) where N = SRow(ntuple(i->(@inbounds ~a.chunks[i]), N))
Base.zero(::Type{SRow{N}}) where N = SRow(ntuple(i->UInt64(0), N))
Base.iszero(a::SRow{N}) where N = all(iszero, a.chunks)
Base.isempty(a::SRow{N}) where N = all(iszero, a.chunks)

Base.@propagate_inbounds function Base.setindex!(a::Base.RefValue{SRow{N}}, b, idx) where N
    i1,i2 = Base.get_chunks_id(idx)
    @boundscheck ((0<i1<=N) || throw(BoundsError(a, idx)))
    @inbounds u = a[].chunks[i1]
    nu = ifelse(b, u | (1<< (i2 & 63)), u & ~(1<< (i2 & 63)) )
    ptr = convert(Ptr{UInt64}, pointer_from_objref(a))
    GC.@preserve a unsafe_store!(ptr, nu, i1)
    b
end

Base.@propagate_inbounds function mk_msk(::Type{SRow{N}}, idx) where N
    i1,i2 = Base.get_chunks_id(idx)
    @boundscheck ((0<i1<=N) || throw(BoundsError(SRow{N}, idx)))
    r = Ref(zero(SRow{N}))
    ptr = convert(Ptr{UInt64}, pointer_from_objref(r))
    GC.@preserve r begin
    for i=1:i1-1
        unsafe_store!(ptr, -1%UInt, i)
    end
    unsafe_store!(ptr, Base._msk_end(idx), i1)
    end
    return r[]
end

#=
function pop_(a::SRow{N}) where N
    i = 1
    @inbounds while i <= N && iszero(a.chunks[i])
        i += 1
    end
    i == N+1 && return -1, a
    @inbounds idx = trailing_zeros(a.chunks[i]) + 1 + (i-1)<<6
    r = Ref(a)
    ptr = convert(Ptr{UInt64}, pointer_from_objref(a))
    unsafe_store!(ptr, _blsr(a.chunks[i]), i)
    return idx, a[]
end
=#

Base.@propagate_inbounds Base.setindex(a::SRow{N}, b, idx) where N = begin r=Ref(a); setindex!(r, b, idx); r[] end

@inline function Base.count_ones(a::SRow{N}) where N
    res = 0
    for i=1:N
        @inbounds res+= count_ones(a.chunks[i])
    end
    res
end
Base.@propagate_inbounds function Base.count_ones(a::SRow{N}, k) where N
    res = 0
    i1, i2 = Base.get_chunks_id(k)
    for i = 1:i1-1
        res += count_ones(a.chunks[i])
    end
    res += count_ones(a.chunks[i1] & Base.msk_end(i1))
    res
end

##SBMGraph
struct SBMGraph{N} <:AbstractBMGraph
    adj_chunks::Matrix{UInt64}
end
@inline SBMGraph(g::BMGraph) = SBMGraph(g, Val(nchunks(g)))
SBMGraph(g::BMGraph, ::Val{N}) where N = SBMGraph{N}(g.adj_chunks)

nchunks(g::SBMGraph{N}) where N = N
Base.@propagate_inbounds function Base.getindex(g::SBMGraph{N}, idx) where N
    @boundscheck checkbounds(g.adj_chunks, 1, idx)
    off = N*(idx-1)
    return SRow(ntuple(i -> (@inbounds g.adj_chunks[i + off]), N))
end
Base.@propagate_inbounds LightGraphs.outneighbors(g::SBMGraph, idx) = g[idx]
function LightGraphs.ne(g::SBMGraph)
    s = sum(count_ones, g.adj_chunks)
    for i=1:nv(g)
        i in g[i] && (s+=1)
    end
    return s>>1
end
Base.@propagate_inbounds LightGraphs.has_edge(g::SBMGraph, s, d) = d in g[s]
Base.@propagate_inbounds LightGraphs.indegree(g::SBMGraph, v::Integer) = count_ones(g[v])
Base.@propagate_inbounds LightGraphs.outdegree(g::SBMGraph, v::Integer) = count_ones(g[v])

LightGraphs.gdistances(g::SBMGraph{N}, s) where N = gdistances!(g, s, fill(0, nv(g)))
function LightGraphs.gdistances!(g::SBMGraph{N}, s, res) where N
    resize!(res, nv(g))
    fill!(res, 0)
    T = SRow{N}
    @inbounds todo = visited = g[s]
    nxt = zero(T)
    dist = 1
    done = false
    while !done
        done = true
        for i in todo
            @inbounds nxt |= g[i] & ~visited
            @inbounds res[i] = dist
            done = false
        end
        done && break
        visited |= nxt
        todo = nxt
        nxt = zero(T)
        dist += 1
    end
    res[s] = 0
    res
end

max_stars_repo_name: wookay/BitmapGraphs

max_issues_repo_name: wookay/BitmapGraphs

max_forks_repo_name: wookay/BitmapGraphs


text:
struct Verlet{TGV, TF<:AbstractFloat, TM} <: SecondOrderNonMetropolisSampler
    ∇V!::TGV
    M::TM
    Δt::TF
end

"""
    Veret(∇V!, M, Δt)

Set up the Verlet integrator.

### Fields

* ∇V!   - In place gradient of the potential
* M     - Mass (either scalar or vector)
* Δt    - Time step
"""
function Veret(∇V!::TGV, M::TM, Δt::TF) where {TGV, TF<:AbstractFloat,TM}
    return Verlet(∇V!,M, Δt);
end

mutable struct VerletState{Tq, Tx<:AbstractVector{Tq}, TF<:AbstractFloat} <:SecondOrderNonMetropolisSamplerState
    x::Tx
    ∇V::Tq
    p_mid::Tq
end

function InitState!(x₀, sampler::Verlet)
    ∇V = similar(x₀[1])
    sampler.∇V!(∇V , x₀);
    return VerletState(x₀, copy(∇V) , similar(x₀[1]));
end

function InitState(x₀, sampler::Verlet)
    ∇V = similar(x₀[1])
    sampler.∇V!(∇V , x₀);
    return VerletState(copy(x₀), copy(∇V) , similar(x₀[1]));
end

function UpdateState!(state::VerletState, sampler::Verlet)
    @. state.p_mid = state.p - 0.5 * sampler.Δt * state.∇V;
    @. state.x = state.x + sampler.Δt * state.p_mid/sampler.M;
    sampler.∇V!(state.∇V, state.x);
    @. state.p = state.p_mid - 0.5 * sampler.Δt * state.∇V;
    state
end

max_stars_repo_name: cortner/JuBasicMD

max_issues_repo_name: cortner/JuBasicMD

max_forks_repo_name: cortner/JuBasicMD


text:
const ENABLE_THREADS = Ref(true)

function disable_threads()
    ENABLE_THREADS[] = false
end

function enable_threads()
    ENABLE_THREADS[] = true
end

macro threads(ex)
    return esc(threads_m(ex))
end

macro threads(min_threads_len::Int, ex)
    return esc(threads_m(ex, min_threads_len))
end

function threads_m(ex::Expr, min_threads_len::Int=1024)
    range = ex.args[1].args[2]
    enable = GlobalRef(BQCESubroutine, :ENABLE_THREADS)
    return quote
        if $enable[] && length($range) > $min_threads_len
            Threads.@threads $ex
        else
            $ex
        end
    end
end


text:
"""
Affine form with center `c`, affine components `γ` and error `Δ`. This
is the AF1 form described in Messine 2002.

Variant where Δ is an interval
"""
struct BasicAffine1{N,T<:Real,S<:AbstractLinearization}
    c::T               # mid-point
    γ::SVector{N,T}    # affine terms
    Δ::Interval{T}     # error term
end

function Base.show(io::IO, C::BasicAffine1{N,T,S}) where {N,T,S}
    print(io, "⟨", C.c, "; ", C.γ, "; ", C.Δ, "⟩")
end

# ==(C::BasicAffine1ine, D::BasicAffine1ine) = C.c == D.c && C.γ == D.γ

"""
Make a `BasicAffine1` based on an interval, which is number `i` of `n` total variables
specifying a linearization mode of type `S`.
"""
function BasicAffine1(X::Interval{T}, ::Val{N}, i::Int, ::S) where {N,T,S}
    c = mid(X)
    r = radius(X)
    γ = SVector{N,T}(ntuple(j->i==j ? r : zero(r), N))
    return BasicAffine1{N,T,S}(c, γ, Interval{T}(0))
end

"""
Make a `BasicAffine1` based on an interval, which is number `i` of `n` total variables
using the default linearization mode.
"""
function BasicAffine1(X::Interval{T}, ::Val{N}, i::Int) where {N,T}
    return BasicAffine1(X, Val(N), i, DEFAULT_LIN())
end

## Promotion rules
promote_rule(::Type{BasicAffine1{N,T,S}}, ::Type{BasicAffine1{N,Q,S}}) where {N, T<:Real, Q<:Real, S} =
    BasicAffine1{N, promote_type(T,Q), S}

promote_rule(::Type{BasicAffine1{N,T,S}}, ::Type{Q}) where {N, T<:Real, Q<:Real, S} =
    BasicAffine1{N, promote_type(T,Q), S}

promote_rule(::Type{BigFloat}, ::Type{BasicAffine1{N,T,S}}) where {N, T<:Real,S} =
    BasicAffine1{N, promote_type(T,BigFloat), S}

# convert methods:
convert(::Type{BasicAffine1{N,T,S}}, x::Bool) where {N,T,S} = convert(BasicAffine1{N,T,S}, Int(x))
convert(::Type{BasicAffine1{N,T,S}}, x::Real) where {N,T,S} = BasicAffine1{N,T,S}(atomic(Interval{T}, x))
convert(::Type{BasicAffine1{N,T,S}}, x::T) where {N,T<:Real,S} = BasicAffine1{N,T,S}(x)
convert(::Type{BasicAffine1{N,T,S}}, x::Interval{T}) where {N,T,S} = x
convert(::Type{BasicAffine1{N,T,S}}, x::Interval) where {N,T,S} = BasicAffine1{N,T,S}(atomic(Interval{T}, x))

for f in (:+,:-)
    @eval function ($f)(x::BasicAffine1{N,T,S}, y::BasicAffine1{N,T,S}) where {N,T<:Real,S<:AbstractLinearization}
        BasicAffine1{N,T,S}(($f).(x.c, y.c), ($f).(x.γ,y.γ), ($f).(x.Δ,y.Δ))
    end
end

interval(C::BasicAffine1{N,T,S}) where {N,T,S} = C.c + sum(abs.(C.γ))*(-1..1) + C.Δ

function *(x::BasicAffine1{N,T,S}, y::BasicAffine1{N,T,S}) where {N,T,S}
    c = x.c * y.c

    γ = x.c .* y.γ + y.c .* x.γ

    Δ = (x.γ ⋅ y.γ) * (0..1)  # ϵ_i^2

    if N > 1
        Δ += sum(x.γ[i] * y.γ[j] for i in 1:N, j in 1:N if i ≠ j) * (-1..1)  # ϵ_i * ϵ_j
    end

    Δ += (x.c + sum(abs.(x.γ))*(-1..1)) * y.Δ
    Δ += (y.c + sum(abs.(y.γ))*(-1..1)) * x.Δ

    Δ += x.Δ * y.Δ

    return BasicAffine1{N,T,S}(c, γ, Δ)
end

# Scalar Operations
*(x::BasicAffine1{N,T,S}, α::T) where {N,T<:Real,S<:AbstractLinearization} = BasicAffine1{N,T,S}(α*x.c, α.*x.γ, α*x.Δ)
*(α::T, x::BasicAffine1{N,T,S}) where {N,T<:Real,S<:AbstractLinearization} = x * α

+(x::BasicAffine1{N,T,S}, α::T) where {N,T<:Real,S<:AbstractLinearization} = BasicAffine1{N,T,S}(α+x.c, x.γ, x.Δ)
+(α::T, x::BasicAffine1{N,T,S}) where {N,T<:Real,S<:AbstractLinearization} = x + α

-(x::BasicAffine1{N,T,S}) where {N,T<:Real,S<:AbstractLinearization} = BasicAffine1{N,T,S}(-x.c, .-(x.γ), -x.Δ)
-(x::BasicAffine1{N,T,S}, α::T) where {N,T<:Real,S<:AbstractLinearization} = BasicAffine1{N,T,S}(x.c - α, x.γ, x.Δ)
-(α::T, x::BasicAffine1{N,T,S}) where {N,T<:Real,S<:AbstractLinearization} = α + (-x)

/(x::BasicAffine1{N,T,S}, α::T) where {N,T<:Real,S<:AbstractLinearization} = BasicAffine1{N,T,S}(x.c/α, x.γ/α, x.Δ/α)
/(α::T, x::BasicAffine1{N,T,S}) where {N,T<:Real,S<:AbstractLinearization} = α*inv(x)

function power_by_squaring(x, p::Integer)
    if p == 1
        return copy(x)
    elseif p == 0
        return one(x)
    elseif p == 2
        return x*x
    elseif p < 0
        isone(x) && return copy(x)
        isone(-x) && return iseven(p) ? one(x) : copy(x)
        Base.throw_domerr_powbysq(x, p)
    end
    t = trailing_zeros(p) + 1
    p >>= t
    x = x
    while (t -= 1) > 0
        x *= x
    end
    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x *= x
        end
        y *= x
    end
    return y
end

# Ran into the following error so re-defined power by squaring:
#=
ERROR: LoadError: MethodError: Cannot `convert` an object of type
  BasicAffine1{2, Float64, MinRange} to an object of type
  Union{}
Closest candidates are:
  convert(::Type{Union{}}, ::Any) at essentials.jl:203
  convert(::Type{var"#s53"} where var"#s53"<:Polynomials.AbstractPolynomial, ::Any) at C:\Users\wilhe\.julia\packages\Polynomials\1aa8e\src\common.jl:309
  convert(::Type{T}, ::Any) where T<:VecElement at baseext.jl:8
  ...
Stacktrace:
  [1] convert(#unused#::Core.TypeofBottom, x::BasicAffine1{2, Float64, MinRange})
    @ Base .\essentials.jl:203
  [2] to_power_type(x::BasicAffine1{2, Float64, MinRange})
    @ Base .\intfuncs.jl:240
  [3] power_by_squaring(x_::BasicAffine1{2, Float64, MinRange}, p::Int64)
    @ Base .\intfuncs.jl:255
  [4] ^(x::BasicAffine1{2, Float64, MinRange}, n::Int64)
    @ AffineArithmetic ~\Desktop\Package Development\AffineArithmetic.jl\src\basic_affine1.jl:14
=#
function ^(x::BasicAffine1{N,T,S}, n::Integer) where {N,T<:Real,S<:AbstractLinearization}

    invert = false

    if n < 0
        invert = true
        n = -n
    end

    result = power_by_squaring(x, n)

    if invert
        result = inv(result)
    end

    return result
end

Base.literal_pow(::typeof(^), x::BasicAffine1, ::Val{p}) where {T,p} = x^p

"General formula for affine approximation of nonlinear functions"
function affine_approx(x::BasicAffine1{N,T,S}, α, ζ, δ) where {N,T<:Real,S<:AbstractLinearization}

    c = α * x.c + ζ
    γ = α .* x.γ
    δ += α * x.Δ  # interval

    return BasicAffine1{MinRange,N,T}(c, γ, δ)
end
# to_power_type(x) = convert(Base._return_type(*, Tuple{typeof(x), typeof(x)}), x)
# to_power_type(x) = convert(Base._return_type(*, Tuple{typeof(x), typeof(x)}), x)
#=
function power_by_squaring(x_, p::Integer)
    x = to_power_type(x_)
    if p == 1
        return copy(x)
    elseif p == 0
        return one(x)
    elseif p == 2
        return x*x
    elseif p < 0
        isone(x) && return copy(x)
        isone(-x) && return iseven(p) ? one(x) : copy(x)
        throw_domerr_powbysq(x, p)
    end
    t = trailing_zeros(p) + 1
    p >>= t
    while (t -= 1) > 0
        x *= x
    end
    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x *= x
        end
        y *= x
    end
    return y
end
=#
function Base.sqrt(x::BasicAffine1{N,T,MinRange}, X=interval(x)) where {N,T<:Real,S<:AbstractLinearization}

    a, b = X.lo, X.hi

    # min-range:  de Figuereido book, pg. 64
    α = 1/(2*√b)
    ζ = (√b)/2
    δ = ((√b - √a)^2/(2*√b))*(-1..1)

    return affine_approx(x, α, ζ, δ)
end

function Base.inv(x::BasicAffine1{N,T,MinRange}, X=interval(x)) where {N,T<:Real,S<:AbstractLinearization}

    a, b = X.lo, X.hi

    # min-range:  de Figuereido book, pg. 70
    α = -1/(b^2)
    d = interval(1/b - α*b, 1/a - α*a)
    ζ = mid(d)
    δ = radius(d) * (-1..1)

    if a < 0
        ζ = -ζ
    end

    return affine_approx(x, α, ζ, δ)
end

function Base.inv(x::BasicAffine1{N,T,Chebyshev}, X=interval(x)) where {N,T<:Real,S<:AbstractLinearization}

    a, b = X.lo, X.hi

    # min-range:  de Figuereido book, pg. 70
    α = -1/(b^2)
    d = interval(1/b - α*b, 1/a - α*a)
    ζ = mid(d)
    δ = radius(d) * (-1..1)

    if a < 0
        ζ = -ζ
    end

    return affine_approx(x, α, ζ, δ)
end

# Directly applies Lemma 1 of Rump 2015 to convex and concave functions of
# fixed monotonicity
f_inverse = Dict{Symbol,Any}()

f_inverse[:exp]   = :log
f_inverse[:exp2]  = :log2
f_inverse[:exp10] = :log10
f_inverse[:expm1] = :log1p

f_inverse[:log]   = :log
f_inverse[:log2]  = :exp2
f_inverse[:log10] = :exp10
f_inverse[:log1p] = :expm1

#f_inverse[:sqrt] = # TODO
f_inverse[:acosh] = :cosh

# Define montonicity
f_curv = Dict{Symbol,Any}()

f_curv[:exp]   = :convex
f_curv[:exp2]  = :convex
f_curv[:exp10] = :convex
f_curv[:expm1] = :convex

f_curv[:log]   = :concave
f_curv[:log2]  = :concave
f_curv[:log10] = :concave
f_curv[:log1p] = :concave

#f_curv[:sqrt] = :concave
f_curv[:acosh] = :concave

for f in (:exp, :exp2, :exp10, :expm1, :log, :log2, :log10, :log1p, :acosh)
    pval = f_curv[f] == :convex ? diffrule(:Base, f, :a) : diffrule(:Base, f, :b)
    finv = f_inverse[f]
    @eval function ($f)(x::BasicAffine1{N,T,MinRange}, X=interval(x)) where {N,T<:Real,S<:AbstractLinearization}

        a, b = X.lo, X.hi

        p = $pval
        q = (($f)(a) + ($f)(b) - p*(a + b))/2
        Δ = abs(($f)(b) - ($f)(a) - p*(b - a))/2

        return affine_approx(x, α, ζ, δ)
    end
    @eval function ($f)(x::BasicAffine1{N,T,Chebyshev}, X=interval(x)) where {N,T<:Real,S<:AbstractLinearization}

        a, b = X.lo, X.hi

        p = (($f)(b) - ($f)(a))/(b - a)
        ξ = ($finv)(p)
        q = (($f)(a) + ($f)(ξ) - p*(a + ξ))/2
        Δ = abs(($f)(ξ) - ($f)(a) - p*(ξ - a))/2

        return affine_approx(x, α, ζ, δ)
    end
end


text:
using Test, GoogleSheetsCSVExporter, CSV

url = "https://docs.google.com/spreadsheets/d/10lDb5TlSZwlsxUSyzLcWDMI9LFRK1S_MRuz8pI_FFxA/edit#gid=0"
table = GoogleSheetsCSVExporter.fromURI(url) |> CSV.File
rows = map(row -> (row.a, row.b, row.c), table)
@test rows[1] == (1, 2, 3)
@test rows[2] === (4, missing, 5)


text:
import TimeZones: astimezone

"""
    AnchoredInterval{P, T}(anchor::T, [inclusivity::Inclusivity]) where {P, T} -> AnchoredInterval{P, T}
    AnchoredInterval{P, T}(anchor::T, [closed_left::Bool, closed_right::Bool]) where {P, T} -> AnchoredInterval{P, T}

`AnchoredInterval` is a subtype of `AbstractInterval` that represents a non-iterable range
or span of values defined not by two endpoints but instead by a single `anchor` point and
the value type `P` which represents the size of the range. When `P` is positive, the
`anchor` represents the lesser endpoint (the beginning of the range); when `P` is negative,
the `anchor` represents the greater endpoint (the end of the range).

The interval represented by an `AnchoredInterval` value may be closed (both endpoints are
included in the interval), open (neither endpoint is included), or half-open. This openness
is defined by an `Inclusivity` value, which defaults to half-open (with the lesser endpoint
included for positive values of `P` and the greater endpoint included for negative values).

### Why?

`AnchoredIntervals` are most useful in cases where a single value is used to stand in for a
range of values. This happens most often with dates and times, where "HE15" is often used as
shorthand for (14:00..15:00].

To this end, `HourEnding` is a type alias for `AnchoredInterval{Hour(-1)}`. Similarly,
`HourBeginning` is a type alias for `AnchoredInterval{Hour(1)}`.

### Rounding

While the user may expect an `HourEnding` or `HourBeginning` value to be anchored to a
specific hour, the constructor makes no guarantees that the anchor provided is rounded:

```julia
julia> HourEnding(DateTime(2016, 8, 11, 2, 30))
HourEnding{DateTime}(2016-08-11T02:30:00, Inclusivity(false, true))
```

The `HE` and `HB` pseudoconstructors round the input up or down to the nearest hour, as
appropriate:

```julia
julia> HE(DateTime(2016, 8, 11, 2, 30))
HourEnding{DateTime}(2016-08-11T03:00:00, Inclusivity(false, true))

julia> HB(DateTime(2016, 8, 11, 2, 30))
HourBeginning{DateTime}(2016-08-11T02:00:00, Inclusivity(true, false))
```

### Example

```julia
julia> AnchoredInterval{Hour(-1)}(DateTime(2016, 8, 11, 12))
HourEnding{DateTime}(2016-08-11T12:00:00, Inclusivity(false, true))

julia> AnchoredInterval{Day(1)}(DateTime(2016, 8, 11))
AnchoredInterval{1 day, DateTime}(2016-08-11T00:00:00, Inclusivity(true, false))

julia> AnchoredInterval{Minute(5)}(DateTime(2016, 8, 11, 12, 30), true, true)
AnchoredInterval{5 minutes, DateTime}(2016-08-11T12:30:00, Inclusivity(true, true))
```

See also: [`Interval`](@ref), [`Inclusivity`](@ref), [`HE`](@ref), [`HB`](@ref)
"""
struct AnchoredInterval{P, T} <: AbstractInterval{T}
    anchor::T
    inclusivity::Inclusivity
end

function AnchoredInterval{P, T}(i::T, x::Bool, y::Bool) where {P, T}
    return AnchoredInterval{P, T}(i, Inclusivity(x, y))
end

# When an interval is anchored to the lesser endpoint, default to Inclusivity(false, true)
# When an interval is anchored to the greater endpoint, default to Inclusivity(true, false)
function AnchoredInterval{P, T}(i::T) where {P, T}
    return AnchoredInterval{P, T}(i::T, Inclusivity(P ≥ zero(P), P ≤ zero(P)))
end

function AnchoredInterval{P, T}(i, inc...) where {P, T}
    return AnchoredInterval{P, T}(convert(T, i), inc...)
end

AnchoredInterval{P}(i::T, inc...) where {P, T} = AnchoredInterval{P, T}(i, inc...)

"""
    HourEnding{T<:TimeType} <: AbstractInterval{T}

A type alias for `AnchoredInterval{Hour(-1), T}` which is used to denote a 1-hour period of
time which ends at a time instant (of type `T`).
"""
const HourEnding{T} = AnchoredInterval{Hour(-1), T} where T <: TimeType
HourEnding(a::T, args...) where T = HourEnding{T}(a, args...)

"""
    HourBeginning{T<:TimeType} <: AbstractInterval{T}

A type alias for `AnchoredInterval{Hour(1), T}` which is used to denote a 1-hour period of
time which begins at a time instant (of type `T`).
"""
const HourBeginning{T} = AnchoredInterval{Hour(1), T} where T <: TimeType
HourBeginning(a::T, args...) where T = HourBeginning{T}(a, args...)

"""
    HE(anchor, args...) -> HourEnding

`HE` is a pseudoconstructor for [`HourEnding`](@ref) that rounds the anchor provided up to the
nearest hour.
"""
HE(a, args...) = HourEnding(ceil(a, Hour), args...)

"""
    HB(anchor, args...) -> HourBeginning

`HB` is a pseudoconstructor for [`HourBeginning`](@ref) that rounds the anchor provided down to the
nearest hour.
"""
HB(a, args...) = HourBeginning(floor(a, Hour), args...)

function Base.copy(x::AnchoredInterval{P, T}) where {P, T}
    return AnchoredInterval{P, T}(anchor(x), inclusivity(x))
end

##### ACCESSORS #####

# We would typically compute `first` and `last` using `min` and `max` respectively, but we
# can get unexpected behaviour if adding the span to the anchor endpoint produces a value
# that is no longer comparable (e.g., `NaN`).

function Base.first(interval::AnchoredInterval{P}) where P
    P < zero(P) ? (interval.anchor + P) : (interval.anchor)
end

function Base.last(interval::AnchoredInterval{P}) where P
    P < zero(P) ? (interval.anchor) : (interval.anchor + P)
end

anchor(interval::AnchoredInterval) = interval.anchor
span(interval::AnchoredInterval{P}) where P = abs(P)

##### CONVERSION #####

function Base.convert(::Type{Interval}, interval::AnchoredInterval{P, T}) where {P, T}
    return Interval{T}(first(interval), last(interval), inclusivity(interval))
end

function Base.convert(::Type{Interval{T}}, interval::AnchoredInterval{P, T}) where {P, T}
    return Interval{T}(first(interval), last(interval), inclusivity(interval))
end

# Conversion methods which currently aren't needed but could prove useful. Commented out
# since these are untested.

#=
function Base.convert(::Type{AnchoredInterval{P, T}}, interval::Interval{T}) where {P, T}
    @assert abs(P) == span(interval)
    anchor = P < zero(P) ? last(interval) : first(interval)
    AnchoredInterval{P, T}(last(interval), inclusivity(interval))
end

function Base.convert(::Type{AnchoredInterval{P}}, interval::Interval{T}) where {P, T}
    @assert abs(P) == span(interval)
    anchor = P < zero(P) ? last(interval) : first(interval)
    AnchoredInterval{P, T}(anchor, inclusivity(interval))
end
=#

function Base.convert(::Type{AnchoredInterval{Ending}}, interval::Interval{T}) where {T}
    AnchoredInterval{-span(interval), T}(last(interval), inclusivity(interval))
end

function Base.convert(::Type{AnchoredInterval{Beginning}}, interval::Interval{T}) where {T}
    AnchoredInterval{span(interval), T}(first(interval), inclusivity(interval))
end

Base.convert(::Type{T}, interval::AnchoredInterval{P, T}) where {P, T} = anchor(interval)

# Date/DateTime attempt to convert to Int64 instead of falling back to convert(T, ...)
Dates.Date(interval::AnchoredInterval{P, Date}) where P = convert(Date, interval)
Dates.DateTime(interval::AnchoredInterval{P, DateTime}) where P = convert(DateTime, interval)

##### DISPLAY #####

function Base.show(io::IO, interval::T) where T <: AnchoredInterval
    if get(io, :compact, false)
        print(io, interval)
    else
        print(io, "$T(")
        show(io, anchor(interval))
        print(io, ", ")
        show(io, inclusivity(interval))
        print(io, ")")
    end
end

function Base.print(io::IO, interval::AnchoredInterval{P, T}) where {P, T <: Union{Date, AbstractDateTime}}
    # Print to io in order to keep properties like :limit and :compact
    if get(io, :compact, false)
        io = IOContext(io, :limit=>true)
    end

    print(io, description(interval))
end

##### ARITHMETIC #####

Base.:+(a::T, b) where {T <: AnchoredInterval} = T(anchor(a) + b, inclusivity(a))

Base.:+(a, b::AnchoredInterval) = b + a
Base.:-(a::AnchoredInterval, b) = a + -b

# Required for StepRange{<:AnchoredInterval}
Base.:-(a::AnchoredInterval, b::AnchoredInterval) = anchor(a) - anchor(b)

Base.:-(a::T, b::AnchoredInterval{P, T}) where {P, T <: Number} = a + -b

function Base.:-(a::AnchoredInterval{P, T}) where {P, T <: Number}
    inc = inclusivity(a)
    AnchoredInterval{-P, T}(-anchor(a), Inclusivity(last(inc), first(inc)))
end

##### EQUALITY #####

function Base.:(==)(a::AnchoredInterval{P, T}, b::AnchoredInterval{P, T}) where {P, T}
    return anchor(a) == anchor(b) && inclusivity(a) == inclusivity(b)
end

# Required for min/max of AnchoredInterval{LaxZonedDateTime} when the anchor is AMB or DNE
function Base.isless(a::AnchoredInterval{P, T}, b::AnchoredInterval{P, T}) where {P, T}
    return (
        anchor(a) < anchor(b) ||
        (anchor(a) == anchor(b) && first(inclusivity(a)) && !first(inclusivity(b)))
    )
end

##### RANGE #####

function Base.steprange_last(start::T, step, stop::T) where {T <: AnchoredInterval}
    T(Base.steprange_last(anchor(start), step, anchor(stop)), inclusivity(start))
end

# Infer step for two-argument StepRange{<:AnchoredInterval}
if VERSION < v"0.7.0-DEV.4003"
    function Base.colon(start::AnchoredInterval{P, T}, stop::AnchoredInterval{P, T}) where {P,T}
        return colon(start, abs(P), stop)
    end
else
    # https://github.com/JuliaLang/julia/pull/26074
    function Base.:(:)(start::AnchoredInterval{P, T}, stop::AnchoredInterval{P, T}) where {P,T}
        return (:)(start, abs(P), stop)
    end
end

function Base.length(r::StepRange{<:AnchoredInterval})
    return length(anchor(r.start):r.step:anchor(r.stop))
end

##### SET OPERATIONS #####

function Base.isempty(interval::AnchoredInterval{P, T}) where {P, T}
    return P == zero(P) && !isclosed(interval)
end

function Base.intersect(a::AnchoredInterval{P, T}, b::AnchoredInterval{Q, T}) where {P,Q,T}
    interval = invoke(intersect, Tuple{AbstractInterval{T}, AbstractInterval{T}}, a, b)

    sp = isa(P, Period) ? canonicalize(typeof(P), span(interval)) : span(interval)
    if P ≤ zero(P)
        anchor = last(interval)
        new_P = -sp
    else
        anchor = first(interval)
        new_P = sp
    end

    return AnchoredInterval{new_P, T}(anchor, inclusivity(interval))
end

##### UTILITIES #####

function canonicalize(target_type::Type{<:Period}, p::P) where P <: Period
    Q, max_val = coarserperiod(P)
    if (value(p) % max_val == 0) && (max_val > 1)
        p = canonicalize(target_type, convert(Q, p))
    end
    return p
end

canonicalize(target_type::Type{P}, p::P) where P <: Period = p

##### TIME ZONES #####

function astimezone(i::AnchoredInterval{P, ZonedDateTime}, tz::TimeZone) where P
    return AnchoredInterval{P, ZonedDateTime}(astimezone(anchor(i), tz), inclusivity(i))
end


text:
using DifferentialEquations, ParameterizedFunctions
van! = @ode_def VanDerPol begin
  dy = μ*((1-x^2)*y - x)
  dx = 1*y
end μ

prob = ODEProblem(van!,[0.0,2.0],(0.0,6.3),1e6)

sol = solve(prob,Tsit5())

sol = solve(prob,alg_hints = [:stiff])

sol = solve(prob)

using Plots; gr()
sol = solve(prob,alg_hints = [:stiff],reltol=1e-6)
plot(sol,denseplot=false)

plot(sol,ylims = (-10.0,10.0))

function lorenz!(du,u,p,t)
    σ,ρ,β = p
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
end
u0 = [1.0,0.0,0.0]
p = (10,28,8/3)
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan,p)

using BenchmarkTools
@btime solve(prob);

@btime solve(prob,alg_hints = [:stiff]);

using SciMLTutorials
SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is 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.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Performs a quick test using ACAS X DASC data set and compares
the results to a reference set.
This is primarily used to see if refactoring has changed results
inadvertantly.
"""
module QuickTest

export quicktest

using GrammarExpts, ACASX_MC, ACASX_GE, ACASX_MCTS
using RLESUtils, FileUtils
using Base.Test

const DIR = dirname(@__FILE__)
const SRCROOT = joinpath(DIR, "..", "..", "..", "data", "quicktest_ref")
const DSTROOT = joinpath(DIR, "..", "..", "..", "results")
const FILES = ["ACASX_MC1/acasx_mc_log_current_best.csv.gz",
            "ACASX_MC1/acasx_mc_log_derivtreevis.json",
            "ACASX_MC1/acasx_mc_log_derivtreevis.tex",
            "ACASX_MC1/acasx_mc_log_members.csv.gz",
            "ACASX_MC1/acasx_mc_log_parameters.csv.gz",
            "ACASX_MC1/acasx_mc_log_result.csv.gz",
            "ACASX_GE/acasx_ge_log_code.csv.gz",
            "ACASX_GE/acasx_ge_log_current_best.csv.gz",
            "ACASX_GE/acasx_ge_log_derivtreevis.json",
            "ACASX_GE/acasx_ge_log_derivtreevis.tex",
            "ACASX_GE/acasx_ge_log_parameters.csv.gz",
            "ACASX_GE/acasx_ge_log_result.csv.gz",
            "ACASX_MCTS/acasx_mcts_log_current_best.csv.gz",
            "ACASX_MCTS/acasx_mcts_log_derivtreevis.json",
            "ACASX_MCTS/acasx_mcts_log_derivtreevis.tex",
            "ACASX_MCTS/acasx_mcts_log_expression.csv.gz",
            "ACASX_MCTS/acasx_mcts_log_members.csv.gz",
            "ACASX_MCTS/acasx_mcts_log_parameters.csv.gz",
            "ACASX_MCTS/acasx_mcts_log_result.csv.gz"]

function quicktest{T<:AbstractString}(b_run=true;
    srcroot::AbstractString=SRCROOT,
    dstroot::AbstractString=DSTROOT,
    files::AbstractVector{T}=FILES)

    if b_run
        acasx_mc1()
        acasx_ge()
        acasx_mcts()
    end

    nfailed = 0
    for f in files
        src = joinpath(srcroot, f)
        dst = joinpath(dstroot, f)
        b = isidentical(src, dst)
        s = b ? "[pass]" : ">FAIL<"
        if !b
            nfailed += 1
        end
        println("$s $f")
    end
    println("$nfailed files failed")
    nfailed
end

end #module


text:
# TODO: Document this
function columns_tuple(
    src::DataFrames.DataFrame,
    names::NTuple{N, Symbol},
) where N
    tuple((src[name] for name in names)...)
end

# TODO: Document this
function tuple_iterator(
    src::DataFrames.DataFrame,
    names::NTuple{N, Symbol},
) where N
    zip(columns_tuple(src, names)...)
end

# TODO: Document this
function gen_column(src, f, column_names)
    map(f, tuple_iterator(src, column_names))
end

# TODO: Document this
function new_column(src, constant)
    [copy(constant) for _ in 1:size(src, 1)]
end

# TODO: Document this
function copy_column(src, column_name)
    copy(src[column_name])
end

# TODO: One wrapper fpr creating a column.


text:
using Tinker
using Gtk.ShortNames, GtkReactive, Graphics, Colors, Images, IntervalSets

# Opens a second window displaying the current value of rectview for the
# selection. Makes it easier to check if the view is correct.

ctx = value(Tinker.active_context)

win = Window("rectview", size(ctx.image,2), size(ctx.image,1));
c = canvas(UserUnit);
push!(win, c);

redraw = draw(c, ctx.rectview, ctx.zr) do cnvs, rv, zr
    #fill!(cnvs, colorant"white")
    set_coordinates(cnvs, zr)
    bg = fill(colorant"white", (size(ctx.image)))
    bview = view(bg,1:size(rv,1),1:size(rv,2))
    copy!(bview,rv)
    copy!(bg, bview)
    copy!(cnvs,bg) # make this not fill the whole canvas
end

showall(win);

nothing


text:
"""
This module defines extensions to the `Base.IO` interface to support:
 - an `unread!` function for pushing excess bytes back into a stream,
 - `startwrite`, `closewrite`, `startread` and `closeread` for streams
    with transactional semantics.
"""
module IOExtras

using ..Sockets
using MbedTLS: MbedException

export bytes, ByteView, CodeUnits, IOError, isioerror,
       unread!,
       startwrite, closewrite, startread, closeread,
       tcpsocket, localport, peerport

"""
    bytes(s::String)

Get a `Vector{UInt8}`, a vector of bytes of a string.
"""
function bytes end
bytes(s::SubArray{UInt8}) = unsafe_wrap(Array, pointer(s), length(s))

const CodeUnits = Union{Vector{UInt8}, Base.CodeUnits}
bytes(s::Base.CodeUnits) = bytes(String(s))
bytes(s::String) = codeunits(s)
bytes(s::SubString{String}) = codeunits(s)

bytes(s::Vector{UInt8}) = s

"""
    isioerror(exception)

Is `exception` caused by a possibly recoverable IO error.
"""
isioerror(e) = false
isioerror(::Base.EOFError) = true
isioerror(::Base.IOError) = true
isioerror(e::ArgumentError) = e.msg == "stream is closed or unusable"
isioerror(::MbedException) = true

"""
The request terminated with due to an IO-related error.

Fields:
 - `e`, the error.
"""
struct IOError <: Exception
    e
    message
end

Base.show(io::IO, e::IOError) = print(io, "IOError(", e.e, " ", e.message, ")\n")

"""
    unread!(::IO, bytes)

Push bytes back into a connection (to be returned by the next read).
"""
function unread!(io::IOBuffer, bytes)
    l = length(bytes)
    if l == 0
        return
    end

    @assert bytes == io.data[io.ptr - l:io.ptr-1]

    if io.seekable
        seek(io, io.ptr - (l + 1))
        return
    end

    println("WARNING: Can't unread! non-seekable IOBuffer")
    println("         Discarding $(length(bytes)) bytes!")
    @assert false
    return
end

function unread!(io, bytes)
    if length(bytes) == 0
        return
    end
    println("WARNING: No unread! method for $(typeof(io))!")
    println("         Discarding $(length(bytes)) bytes!")
    return
end

_doc = """
    startwrite(::IO)
    closewrite(::IO)
    startread(::IO)
    closeread(::IO)

Signal start/end of write or read operations.
"""
"$_doc"
startwrite(io) = nothing
"$_doc"
closewrite(io) = nothing
"$_doc"
startread(io) = nothing
"$_doc"
closeread(io) = nothing

using MbedTLS: SSLContext
tcpsocket(io::SSLContext)::TCPSocket = io.bio
tcpsocket(io::TCPSocket)::TCPSocket = io

localport(io) = try !isopen(tcpsocket(io)) ? 0 :
                    Sockets.getsockname(tcpsocket(io))[2]
                catch
                    0
                end

peerport(io) = try !isopen(tcpsocket(io)) ? 0 :
                  Sockets.getpeername(tcpsocket(io))[2]
               catch
                   0
               end

end

const ByteView = typeof(view(UInt8[], 1:0))

"""
Read from an `IO` stream until `find_delimiter(bytes)` returns non-zero.
Return view of bytes up to the delimiter.
"""
function Base.readuntil(io::IO,
                        find_delimiter::Function #= Vector{UInt8} -> Int =#
                       )::ByteView

    # Fast path, buffer already contains delimiter...
    if !eof(io)
        bytes = readavailable(io)
        if (l = find_delimiter(bytes)) > 0
            if l < length(bytes)
                unread!(io, view(bytes, l+1:length(bytes)))
            end
            return view(bytes, 1:l)
        end

        # Otherwise, wait for delimiter...
        buf = Vector{UInt8}(bytes)
        while !eof(io)
            bytes = readavailable(io)
            append!(buf, bytes)
            if (l = find_delimiter(buf)) > 0
                if l < length(buf)
                    n = length(buf) - l
                    bl = length(bytes)
                    unread!(io, view(bytes, 1+bl-n:bl))
                end
                return view(buf, 1:l)
            end
        end
    end

    throw(EOFError())
end


text:
ENV["GKSwstype"] = "nul"
using Test
using IntervalSets

import Plots

import ClimaCore
import ClimaCorePlots

OUTPUT_DIR = mkpath(get(ENV, "CI_OUTPUT_DIR", tempname()))

@testset "spectral element 2D cubed-sphere" begin
    R = 6.37122e6

    domain = ClimaCore.Domains.SphereDomain(R)
    mesh = ClimaCore.Meshes.EquiangularCubedSphere(domain, 6)
    grid_topology = ClimaCore.Topologies.Topology2D(mesh)
    quad = ClimaCore.Spaces.Quadratures.GLL{5}()
    space = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)
    coords = ClimaCore.Fields.coordinate_field(space)

    u = map(coords) do coord
        u0 = 20.0
        α0 = 45.0
        ϕ = coord.lat
        λ = coord.long

        uu = u0 * (cosd(α0) * cosd(ϕ) + sind(α0) * cosd(λ) * sind(ϕ))
        uv = -u0 * sind(α0) * sind(λ)
        ClimaCore.Geometry.UVVector(uu, uv)
    end

    field_fig = Plots.plot(u.components.data.:1)
    @test field_fig !== nothing

    fig_png = joinpath(OUTPUT_DIR, "2D_cubed_sphere_field.png")
    Plots.png(field_fig, fig_png)
    @test isfile(fig_png)
end

@testset "spectral element rectangle 2D" begin
    domain = ClimaCore.Domains.RectangleDomain(
        ClimaCore.Geometry.XPoint(0) .. ClimaCore.Geometry.XPoint(2π),
        ClimaCore.Geometry.YPoint(0) .. ClimaCore.Geometry.YPoint(2π),
        x1periodic = true,
        x2periodic = true,
    )

    n1, n2 = 2, 2
    Nq = 4
    mesh = ClimaCore.Meshes.RectilinearMesh(domain, n1, n2)
    grid_topology = ClimaCore.Topologies.Topology2D(mesh)
    #quad = ClimaCore.Spaces.Quadratures.GLL{Nq}()
    quad = ClimaCore.Spaces.Quadratures.ClosedUniform{Nq + 1}()
    space = ClimaCore.Spaces.SpectralElementSpace2D(grid_topology, quad)
    coords = ClimaCore.Fields.coordinate_field(space)

    space_fig = Plots.plot(space)
    @test space_fig !== nothing

    sinxy = map(coords) do coord
        cos(coord.x + coord.y)
    end

    field_fig = Plots.plot(sinxy)
    @test field_fig !== nothing

    space_png = joinpath(OUTPUT_DIR, "2D_rectangle_space.png")
    field_png = joinpath(OUTPUT_DIR, "2D_rectangle_field.png")
    Plots.png(space_fig, space_png)
    Plots.png(field_fig, field_png)
    @test isfile(space_png)
    @test isfile(field_png)
end

@testset "hybrid finite difference / spectral element 2D" begin
    FT = Float64
    helem = 10
    velem = 40
    npoly = 4

    vertdomain = ClimaCore.Domains.IntervalDomain(
        ClimaCore.Geometry.ZPoint{FT}(0),
        ClimaCore.Geometry.ZPoint{FT}(1000);
        boundary_tags = (:bottom, :top),
    )
    vertmesh = ClimaCore.Meshes.IntervalMesh(vertdomain, nelems = velem)
    vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(vertmesh)

    horzdomain = ClimaCore.Domains.IntervalDomain(
        ClimaCore.Geometry.XPoint{FT}(-500) ..
        ClimaCore.Geometry.XPoint{FT}(500),
        periodic = true,
    )
    horzmesh = ClimaCore.Meshes.IntervalMesh(horzdomain; nelems = helem)
    horztopology = ClimaCore.Topologies.IntervalTopology(horzmesh)

    quad = ClimaCore.Spaces.Quadratures.GLL{npoly + 1}()
    horzspace = ClimaCore.Spaces.SpectralElementSpace1D(horztopology, quad)

    hv_center_space = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(
        horzspace,
        vert_center_space,
    )
    hv_face_space =
        ClimaCore.Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)

    coords = ClimaCore.Fields.coordinate_field(hv_center_space)

    xcoords_fig = Plots.plot(coords.x)
    @test xcoords_fig !== nothing

    zcoords_fig = Plots.plot(coords.z)
    @test zcoords_fig !== nothing

    xcoords_png = joinpath(OUTPUT_DIR, "hybrid_xcoords_center_field.png")
    zcoords_png = joinpath(OUTPUT_DIR, "hybrid_zcoords_center_field.png")
    Plots.png(xcoords_fig, xcoords_png)
    Plots.png(zcoords_fig, zcoords_png)
    @test isfile(xcoords_png)
    @test isfile(zcoords_png)
end

@testset "hybrid finite difference / spectral element 3D" begin
    FT = Float64
    xelem = 10
    yelem = 5
    velem = 40
    npoly = 4

    vertdomain = ClimaCore.Domains.IntervalDomain(
        ClimaCore.Geometry.ZPoint{FT}(0),
        ClimaCore.Geometry.ZPoint{FT}(1000);
        boundary_tags = (:bottom, :top),
    )
    vertmesh = ClimaCore.Meshes.IntervalMesh(vertdomain, nelems = velem)
    vert_center_space = ClimaCore.Spaces.CenterFiniteDifferenceSpace(vertmesh)

    xdomain = ClimaCore.Domains.IntervalDomain(
        ClimaCore.Geometry.XPoint{FT}(-500) ..
        ClimaCore.Geometry.XPoint{FT}(500),
        periodic = true,
    )
    ydomain = ClimaCore.Domains.IntervalDomain(
        ClimaCore.Geometry.YPoint{FT}(-100) ..
        ClimaCore.Geometry.YPoint{FT}(100),
        periodic = true,
    )

    horzdomain = ClimaCore.Domains.RectangleDomain(xdomain, ydomain)
    horzmesh = ClimaCore.Meshes.RectilinearMesh(horzdomain, xelem, yelem)
    horztopology = ClimaCore.Topologies.Topology2D(horzmesh)

    quad = ClimaCore.Spaces.Quadratures.GLL{npoly + 1}()
    horzspace = ClimaCore.Spaces.SpectralElementSpace2D(horztopology, quad)

    hv_center_space = ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace(
        horzspace,
        vert_center_space,
    )
    hv_face_space =
        ClimaCore.Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)

    coords = ClimaCore.Fields.coordinate_field(hv_center_space)

    xcoords_fig = Plots.plot(coords.x, slice = (:, 0.0, :))
    @test xcoords_fig !== nothing

    ycoords_fig = Plots.plot(coords.y, slice = (0.0, :, :))
    @test ycoords_fig !== nothing

    xzcoords_fig = Plots.plot(coords.z, slice = (:, 0.0, :))
    @test xzcoords_fig !== nothing

    yzcoords_fig = Plots.plot(coords.z, slice = (0.0, :, :))
    @test yzcoords_fig !== nothing

    xcoords_png = joinpath(OUTPUT_DIR, "hybrid_xcoords_center_field.png")
    ycoords_png = joinpath(OUTPUT_DIR, "hybrid_ycoords_center_field.png")
    xzcoords_png = joinpath(OUTPUT_DIR, "hybrid_xzcoords_center_field.png")
    yzcoords_png = joinpath(OUTPUT_DIR, "hybrid_yzcoords_center_field.png")

    Plots.png(xcoords_fig, xcoords_png)
    Plots.png(ycoords_fig, ycoords_png)
    Plots.png(xzcoords_fig, xzcoords_png)
    Plots.png(yzcoords_fig, yzcoords_png)

    @test isfile(xcoords_png)
    @test isfile(ycoords_png)
    @test isfile(xzcoords_png)
    @test isfile(yzcoords_png)
end


text:
#FEAST for general nonlinear eigenvalue problems, using Beyn's algorithm for Rayleigh Ritz.
function nlfeast_beyn(Tf,x0,nc,emid,ra,rb,eps,maxit;log=false)
    #Tf is single value nonlinear residual T(z)

    #integrand = (I-T^-1(z)T(l))x(zI-l)^-1
    function integrand(z,x,lest,data,resvecs)
            (n,m0)=size(x)
            TinvTxl=zeros(ComplexF64,n,m0)
            Tfz=Tf(z)

            TinvTxl=\(Tfz,resvecs)

            gamma=spdiagm(0 => 1 ./(z.-lest))
            Qk=x*gamma-TinvTxl*gamma
            return Qk
    end

    #solve projected problem Q'*T(z)Q=0 using Beyn's algorithm
    function rrsolve(Q)
        (n,m0)=size(Q)

        v0=rand(ComplexF64,m0,m0)
        Tfr(z)=Q'*Tf(z)*Q

        (lest,xq)=beyn(Tfr,v0,5*nc,emid,ra,rb)

        return (lest,Q*xq)
    end

    #Block method for applying residual
    function blockTf(l,x)
        (n,m0)=size(x)
        out=zero(x)
        for i in 1:m0
            out[:,i]=Tf(l[i])*x[:,i]
        end
        return out
    end

    return feast_core(blockTf,integrand,rrsolve,x0,nc,emid,ra,rb,eps,maxit;log=log)
end

#iterative FEAST for general nonlinear eigenvalue problems, using Beyn's algorithm for Rayleigh Ritz.
function inlfeast_beyn(Tf,x0,alpha,isMaxit,nc,emid,ra,rb,eps,maxit;log=false)
    #Tf is single value nonlinear residual T(z)

    #integrand = (I-T^-1(z)T(l))x(zI-l)^-1
    function integrand(z,x,lest,data,resvecs)
            (n,m0)=size(x)
            TinvTxl=zeros(ComplexF64,n,m0)
            Tfz=Tf(z)

            maxits=0

            #TinvTxl=\(Tfz,resvecs)
            int=zeros(ComplexF64,n,m0)
            for i in 1:m0
                #(int[:,i],history)=bicgstabl(M,resvecs[:,i],1,max_mv_products=isMaxit,tol=alpha,initial_zero=true,log=true)
                #(int[:,i],history)=idrs(M,rhs[:,i];maxiter=isMaxit,tol=alpha,log=true)
                (int[:,i],history)=gmres(Tfz,resvecs[:,i],restart=isMaxit,tol=alpha,initially_zero=true,maxiter=isMaxit,log=true)
                #(int[:,i],history)=minres(M,resvecs[:,i],tol=alpha,initially_zero=true,maxiter=isMaxit,log=true)
                nlinits=size(history[:resnorm],1)
                data[:linIts][i,nc,data[:iterations]]=nlinits
                data[:linResiduals][i,nc,data[:iterations]]=history[:resnorm][nlinits]
                #data[:linResiduals][i,nc,data[:iterations]]=norm(resvecs[:,i]-M*int[:,i])/norm(resvecs[:,i])

                if(nlinits>maxits)
                    maxits=nlinits
                end
            end

            return (x-int)*spdiagm(0 => 1 ./(z.-lest))
    end

    #solve projected problem Q'*T(z)Q=0 using Beyn's algorithm
    function rrsolve(Q)
        (n,m0)=size(Q)

        v0=rand(ComplexF64,m0,m0)
        Tfr(z)=Q'*Tf(z)*Q

        (lest,xq)=beyn(Tfr,v0,5*nc,emid,ra,rb)

        return (lest,Q*xq)
    end

    #Block method for applying residual
    function blockTf(l,x)
        (n,m0)=size(x)
        out=zero(x)
        for i in 1:m0
            out[:,i]=Tf(l[i])*x[:,i]
        end
        return out
    end

    return feast_core(blockTf,integrand,rrsolve,x0,nc,emid,ra,rb,eps,maxit;log=log)
end

#FEAST for quadratic eigenvalue problem
function nlfeast_quad(M,D,K,x0,nc,emid,ra,rb,eps,maxit;log=false)
    (n,m0)=size(x0)

    #pardiso for sparse systems
    #ps=MKLPardisoSolver()

    Tfm(z)=M*z^2+z*D+K

    #Q=integral (I-T^-1(z)T(X,Lambda))X(zI-Lambda)^-1 dz
    function integrand(z,x,lest,data,resvecs)
            (n,m0)=size(x)

            Txl=M*x*diagm(0 => lest.^2)+D*x*diagm(0 => lest)+K*x
            #Txl=M*x2+D*x1+K*x
            #Txl=resvecs
            Tz=z^2*M+z*D+K

            TinvTxl=\(Tz,Txl)

            gamma=diagm(0 => 1 ./(z.-lest))
            Qk=x*gamma-TinvTxl*gamma
            return Qk
    end

    #solve reduced nonlinear problem with linearization
    #select the m0 solutions that are closest to the center of the contour
    rrsolve(Q)=Quadrrsolve(Q,M,D,K,emid,ra,rb)

    Tf(l,x)=M*x*diagm(0 => l.^2)+D*x*diagm(0 => l)+K*x

    return feast_core(Tf,integrand,rrsolve,x0,nc,emid,ra,rb,eps,maxit,log=log)
end

max_stars_repo_name: rleegates/feastjl

max_issues_repo_name: rleegates/feastjl

max_forks_repo_name: rleegates/feastjl