text:
using BlockArrays, BandedMatrices, BlockBandedMatrices, FillArrays, SparseArrays, Test, ArrayLayouts , LinearAlgebra
import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolsupport, blockrowsupport, colsupport, rowsupport, isbandedblockbanded, bandeddata

@testset "BandedBlockBandedMatrix" begin
    @testset "constructors" begin
        l , u = 1,1
        λ , μ = 1,1
        N = M = 4
        cols = rows = 1:N

        @test Matrix(BandedBlockBandedMatrix(Zeros(sum(rows),sum(cols)), rows,cols, (l,u), (λ,μ))) ==
            Array(BandedBlockBandedMatrix(Zeros(sum(rows),sum(cols)), rows,cols, (l,u), (λ,μ))) ==
            zeros(Float64, 10, 10)

        @test Matrix(BandedBlockBandedMatrix{Int}(Zeros(sum(rows),sum(cols)), rows,cols, (l,u), (λ,μ))) ==
            zeros(Int, 10, 10)

        @test Matrix(BandedBlockBandedMatrix(Eye(sum(rows)), rows,cols, (l,u), (λ,μ))) ==
            Matrix{Float64}(I, 10, 10)

        @test Matrix(BandedBlockBandedMatrix{Int}(Eye(sum(rows)), rows,cols, (l,u), (λ,μ))) ==
            Matrix{Int}(I, 10, 10)

        @test Matrix(BandedBlockBandedMatrix(I, rows,cols, (l,u), (λ,μ))) ==
            Matrix{Float64}(I, 10, 10)

        @test Matrix(BandedBlockBandedMatrix{Int}(I, rows,cols, (l,u), (λ,μ))) ==
            Matrix{Int}(I, 10, 10)

        A = [1 2 3 4 5; 6 7 8 9 10; 11 12 13 14 15; 16 17 18 19 20; 21 22 23 24 25]
        B = BandedBlockBandedMatrix(A, [2,3], [2,3], (0,1), (1,1))
        @test Matrix(B) ==
            [1 2 3 4 0; 6 7 8 9 10; 0 0 13 14 0; 0 0 18 19 20; 0 0 0 24 25]

        l,u,λ,μ = 0,0,4,0
        cols = 1:2:3
        rows = 2*cols
        A = BandedBlockBandedMatrix(fill(1.0,sum(rows),sum(cols)), rows,cols, (l,u), (λ,μ))
        B = [1.0 0.0 0.0 0.0;
             1.0 0.0 0.0 0.0;
             0.0 1.0 0.0 0.0;
             0.0 1.0 1.0 0.0;
             0.0 1.0 1.0 1.0;
             0.0 1.0 1.0 1.0;
             0.0 1.0 1.0 1.0;
             0.0 0.0 1.0 1.0]
        @test Matrix(A) == B

        @test Matrix(sparse(A)) == B
    end

    @testset "col and rowranges" begin
        l , u = 1,1
        λ , μ = 1,1
        N = M = 4
        cols = rows = 1:N
        data = reshape(collect(1:(λ+μ+1)*(l+u+1)*sum(cols)), ((λ+μ+1)*(l+u+1), sum(cols)))
        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))
        @test Matrix(BlockBandedMatrix(A)) == Matrix(A)

        @test blockbandwidths(A) == (l,u)
        @test BlockBandedMatrices.subblockbandwidths(A) == (l,u)

        @test blockrowsupport(A, 1) == Block.(1:2)
        @test blockrowsupport(A, 2) == Block.(1:3)
        @test blockrowsupport(A, 3) == Block.(2:4)

        @test blockcolsupport(A, 1) == Block.(1:2)
        @test blockcolsupport(A, 2) == Block.(1:3)
        @test blockcolsupport(A, 3) == Block.(2:4)

        @test rowsupport(A,1) == 1:3
        @test rowsupport(A,2) == 1:6
        @test rowsupport(A,3) == 1:6
        @test rowsupport(A,4) == 2:10

        @test colsupport(A,1) == 1:3
        @test colsupport(A,2) == 1:6
        @test colsupport(A,3) == 1:6
        @test colsupport(A,4) == 2:10

        l , u = 2,1
        λ , μ = 1,2
        N = M = 4
        cols = rows = 1:N
        data = reshape(Vector(1:(λ+μ+1)*(l+u+1)*sum(cols)), (λ+μ+1, (l+u+1)*sum(cols)))
        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))

        @test A.l == l == blockbandwidth(A,1)
        @test A.u == u == blockbandwidth(A,2)
        @test blockbandwidths(A) == (l, u)

        @test A.λ == λ == subblockbandwidth(A,1)
        @test A.μ == μ == subblockbandwidth(A,2)
        @test subblockbandwidths(A) == (λ, μ)

        @test blockrowsupport(A, 1) == Block.(1:2)
        @test blockrowsupport(A, 2) == Block.(1:3)
        @test blockrowsupport(A, 3) == Block.(1:4)
        @test blockrowsupport(A, 4) == Block.(2:4)

        @test blockcolsupport(A, 1) == Block.(1:3)
        @test blockcolsupport(A, 2) == Block.(1:4)
        @test blockcolsupport(A, 3) == Block.(2:4)
        @test blockcolsupport(A, 4) == Block.(3:4)

        @test rowsupport(A,1) == 1:3
        @test rowsupport(A,2) == 1:6
        @test rowsupport(A,3) == 1:6
        @test rowsupport(A,4) == 1:10

        @test colsupport(A,1) == 1:6
        @test colsupport(A,2) == 1:10
        @test colsupport(A,3) == 1:10
        @test colsupport(A,4) == 2:10

        l , u = -1,1
        λ , μ = 0,1
        rows = 1:5
        cols = 1:6

        data = reshape(Vector{Float64}(1:(λ+μ+1)*(l+u+1)*sum(cols)), (λ+μ+1, (l+u+1)*sum(cols)))
        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))
        @test_throws BandError A[1,1] = 5

        @test A[1,2] == 4
        @test A[Block(2,2)] == [0 0; 0 0]
        @test A[Block(2,3)]  == [8.0 9.0 0.0; 0.0 10.0 11.0]
        @test bandwidths(view(A,Block(2,3))) == bandwidths(A[Block(2,3)]) == (0,1)

        @test blockrowsupport(A, 1) == Block.(2:2)
        @test blockrowsupport(A, 2) == Block.(3:3)
        @test blockrowsupport(A, 3) == Block.(4:4)
        @test blockrowsupport(A, 4) == Block.(5:5)

        @test blockcolsupport(A, 1) == Block.(1:0)
        @test blockcolsupport(A, 2) == Block.(1:1)
        @test blockcolsupport(A, 3) == Block.(2:2)

        @test rowsupport(A,1) == 2:3
        @test rowsupport(A,2) == 4:6
        @test rowsupport(A,3) == 4:6
        @test rowsupport(A,4) == 7:10

        @test colsupport(A,1) == 1:0
        @test colsupport(A,2) == 1:1
        @test colsupport(A,3) == 1:1
        @test colsupport(A,4) == 2:3

        l , u = -1,1
        λ , μ = -1,1
        rows = 1:5
        cols = 1:6

        data = reshape(Vector{Float64}(1:(λ+μ+1)*(l+u+1)*sum(cols)), (λ+μ+1, (l+u+1)*sum(cols)))
        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))
        @test_throws BandError A[1,1] = 5

        @test blockcolsupport(A, 1) == Block.(1:0)
        @test blockcolsupport(A, 2) == Block.(1:1)
        @test blockcolsupport(A, 3) == Block.(2:2)

        @test colsupport(A,1) == 1:0
        @test colsupport(A,2) == 1:1
        @test colsupport(A,3) == 1:1
        @test colsupport(A,4) == 2:3
    end

    @testset "block indexing" begin
        l , u = 1,1
        λ , μ = 1,1
        N = M = 4
        cols = rows = 1:N
        data = reshape(collect(1:(λ+μ+1)*(l+u+1)*sum(cols)), ((λ+μ+1)*(l+u+1), sum(cols)))
        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))

        @test A[Block(1), Block(1)] isa BandedMatrix
        @test A[Block(1), Block(1)] == A[Block(1,1)] == BlockArrays.getblock(A, 1, 1) == BandedMatrix(view(A, Block(1,1)))
        @test A[1,1] == view(A,Block(1),Block(1))[1,1] == view(A,Block(1,1))[1,1] == A[Block(1,1)][1,1]  == A[Block(1),Block(1)][1,1] == 5
        @test A[2,1] == view(A,Block(2),Block(1))[1,1] == view(A,Block(2,1))[1,1] == 8
        @test A[3,1] == view(A,Block(2),Block(1))[2,1] == 9
        @test A[4,1] == 0
        @test A[1,2] == view(A,Block(1,2))[1,1] == 11
        @test A[1,3] == view(A,Block(1,2))[1,2] == view(A,Block(1,2))[2] == 19

        @test bandwidths(view(A, Block(3),Block(1))) == (-720,-720)
        @test isempty(bandeddata(view(A, Block(3),Block(1))))

        @test A[Block(3,1)] == view(A, Block(3),Block(1)) == zeros(3,1)
        @test A[Block(3,1)] ≈ view(A, Block(3),Block(1)) ≈ zeros(3,1)
        @test_throws BandError view(A, Block(3),Block(1))[1,1] = 4
        @test_throws BlockBoundsError view(A, Block(5,1))
    end

    @testset "indexing" begin
        l , u = 1,1
        λ , μ = 1,1
        N = M = 10
        cols = rows = 1:N
        data = reshape(Vector{Float64}(1:(λ+μ+1)*(l+u+1)*sum(cols)), ((λ+μ+1)*(l+u+1),sum(cols)))
        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))

        A[1,1] = 5
        @test A[1,1] == 5

        @test_throws BandError A[1,4] = 5
        A[1,4] = 0
        @test A[1,4] == 0

        # TODO: return a BandedMatrix
        @test A[1:10,1:10] ≈ Matrix(A)[1:10,1:10]
    end

    @testset "banded matrix interface for blocks" begin
        l , u = 1,1
        λ , μ = 1,1
        N = M = 10
        cols = rows = 1:N
        data = reshape(collect(1:(λ+μ+1)*(l+u+1)*sum(cols)), ((λ+μ+1)*(l+u+1), sum(cols)))
        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))

        # test blocks
        V = view(A, Block(1,1))
        @test_throws BoundsError V[2,1]

        # BandedMatrix interface
        @test isbanded(V)
        @test bandwidths(V) == BlockBandedMatrices.subblockbandwidths(A)
        @test V[band(0)] ≈ view(V, band(0)) ≈ A[1:1,1:1]

        # test views of blocks fulfill BnadedMatrix interface
        @test BandedMatrices.inbands_getindex(V, 1, 1) == V[1,1] == 5
        BandedMatrices.inbands_setindex!(V, -1, 1, 1)
        @test A[1,1] == -1
        # these should throw errors but inbands turns it off

        V = view(A, Block(3,4))
        @test V[3,1] == 0
        @test_throws BandError V[3,1] = 5

        view(V, band(0)) .= -3
        @test all(A[Block(3,4)][band(0)] .== -3)

        @test BandedMatrix(V) isa BandedMatrix{Int,Matrix{Int}}
        @test BandedMatrix{Float64}(V) isa BandedMatrix{Float64,Matrix{Float64}}
        @test BandedMatrix{Float64}(BandedMatrix(V)) == BandedMatrix{Float64}(V) ==
                convert(BandedMatrix{Float64}, V)
        @test A[4:6,7:10] ≈ BandedMatrix(V)

        @test A[Block(3,4)].l == A.λ
        @test A[Block(3,4)].u == A.μ

        A[Block(3,4)] = BandedMatrix(Ones{Int}(3,4),(1,1))
        @test A[Block(3,4)] == BandedMatrix(Ones{Int}(3,4),(1,1))

        l , u = 2,1
        λ , μ = 1,2
        N = M = 4
        cols = rows = 1:N
        data = reshape(Vector(1:(λ+μ+1)*(l+u+1)*sum(cols)), (λ+μ+1, (l+u+1)*sum(cols)))
        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))

        @test_throws BandError A[1,4] = 5
        @test_throws BandError view(A, Block(1,3))[2] = 5

        lu = (l , u) = -1,1
        λμ = (λ , μ) = -1,1
        rows = 1:5
        cols = 1:6

        data = reshape(Vector{Float64}(1:(λ+μ+1)*(l+u+1)*sum(cols)), (λ+μ+1, (l+u+1)*sum(cols)))
        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))
        @test_throws BandError A[1,1] = 5

        @test A[1,3] == 3
        @test A[Block(2,2)] == [0 0; 0 0]
        @test A[Block(2,3)]  == [0 5 0; 0 0 6]
        @test bandwidths(A[Block(2,3)]) == (-1,1)
    end

    @testset "banded-block-banded matrix interface for blockranges" begin
        l , u = 1,1
        λ , μ = 1,1
        N = M = 10
        cols = rows = 1:N
        data = reshape(collect(1.0:(λ+μ+1)*(l+u+1)*sum(cols)), ((λ+μ+1)*(l+u+1), sum(cols)))
        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))
        V = view(A, Block.(2:3), Block.(3:4))
        @test isbandedblockbanded(V)

        B = BandedBlockBandedMatrix(V)
        @test B isa BandedBlockBandedMatrix
        @test blockbandwidths(V) == blockbandwidths(B) == (2,0)
        @test subblockbandwidths(V) == subblockbandwidths(B) == subblockbandwidths(A) == (λ,μ)
        @test B == V == A[Block.(2:3), Block.(3:4)]

        @test A[Block.(2:3), Block.(3:4)] isa BandedBlockBandedMatrix

        x = randn(size(B,2))
        y = similar(x, size(B,1))
        @test all((similar(y) .= MulAdd(B, x)) .=== (similar(y) .= MulAdd(V,x)))
    end

    @testset "BLAS arithmetic" begin
        l , u = 1,1
        λ , μ = 1,1
        N = M = 10
        cols = rows = fill(1000,N)
        data = reshape(Vector{Float64}(1:(λ+μ+1)*(l+u+1)*sum(cols)), ((λ+μ+1)*(l+u+1), sum(cols)))
        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))

        V = view(A, Block(N,N))

        AN = A[Block(N,N)]
        BLAS.axpy!(2.0, V, V)
        @test A[Block(N,N)] ≈ 3AN

        Y = zeros(cols[N], cols[N])
        BLAS.axpy!(2.0, V, Y)
        @test Y ≈ 2A[Block(N,N)]

        Y = BandedMatrix(Zeros(cols[N], cols[N]), (λ, μ))
        BLAS.axpy!(2.0, V, Y)
        @test Y ≈ 2A[Block(N,N)]

        Y = BandedMatrix(Zeros(cols[N], cols[N]), (λ+1, μ+1))
        BLAS.axpy!(2.0, V, Y)
        @test Y ≈ 2A[Block(N,N)]

        Y = BandedMatrix(Zeros(cols[N], cols[N]), (0, 0))
        @test_throws BandError BLAS.axpy!(2.0, V, Y)
    end

    @testset "Float32"  begin
        A = BandedBlockBandedMatrix{Float32}(Zeros{Float32}(10,10),
                                    fill(1,10), fill(1,10), (1,1), (1,1))

        @test eltype(A) == Float32

        A = BandedBlockBandedMatrix(Zeros{Float32}(10,10),
                                    fill(1,10), fill(1,10), (1,1), (1,1))

        @test eltype(A) == Float32
    end

    @testset "fill and copy" begin
        l , u = 1,1
        λ , μ = 1,1
        N = M = 10
        cols = rows = 1:N
        data = randn((λ+μ+1)*(l+u+1), sum(cols))
        A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))

        V = view(A,Block(2,3))
        @test_throws BandError fill!(V, 2.0)
        fill!(V, 0)
        @test A[2:3,4:6] == zeros(2,3)

        dataA = randn((λ+μ+1)*(l+u+1), sum(cols))
        A = _BandedBlockBandedMatrix(copy(dataA), rows,cols, (l,u), (λ,μ))

        K,J = 2,1
        V = view(A,Block(K),Block(J))
        fill!(V, 2.0)
        @test Matrix(V) == fill(2.0, 2,1)
        fill!(view(A,Block(K),Block(J)), 0.0)
        @test Matrix(view(A,Block(K),Block(J))) == zeros(2,1)

        K,J = 3,1
        @test_throws BandError fill!(view(A,Block(K),Block(J)), 2.0)
        fill!(view(A,Block(K),Block(J)), 0.0)
        @test Matrix(view(A,Block(K),Block(J))) == zeros(3,1)

        @test_throws BandError fill!(A, 2.0)
        fill!(A, 0.0)
        @test Matrix(A) == zeros(size(A))

        dataA = randn((λ+μ+1)*(l+u+1), sum(cols))
        A = _BandedBlockBandedMatrix(copy(dataA), rows,cols, (l,u), (λ,μ))

        dataB = randn((λ+μ+3)*(l+u+3), sum(cols))
        B = _BandedBlockBandedMatrix(copy(dataB), rows,cols, (l+1,u+1), (λ+1,μ+1))

        B = _BandedBlockBandedMatrix(copy(dataB), rows,cols, (l+1,u+1), (λ+1,μ+1))
        copyto!(view(B, Block(N,N)), view(A, Block(N,N)))
        @test B[Block(N,N)] == A[Block(N,N)]

        B = _BandedBlockBandedMatrix(copy(dataB), rows,cols, (l+1,u+1), (λ+1,μ+1))
        copyto!(B, A)
        @test norm(B[Block(3,1)]) == 0
        @test  B ≈ A
        @test Matrix(B) ≈ Matrix(A)
    end

    @testset "Zero bands" begin
        B = BandedBlockBandedMatrix{Float64}(undef, 1:5,1:5, (-1,-1), (-1,-1))
        @test Matrix(B) == zeros(size(B))

        B = BandedBlockBandedMatrix{Float64}(undef, 1:5,1:5, (-1,-1), (1,-1))
        @test Matrix(B) == zeros(size(B))

        B = BandedBlockBandedMatrix{Float64}(undef, 1:5,1:5, (1,-1), (-1,-1))
        @test Matrix(B) == zeros(size(B))
    end

    @testset "BandedBlockBanded with BlockMatrix" begin
        WithBlockMatrix{T} = BandedBlockBandedMatrix{T, BlockMatrix{T, Matrix{Matrix{T}}}}
        args = [1, 2, 3], [2, 2, 1], (1, 1), (1, 1)
        A = WithBlockMatrix{Int64}(undef, args...)
        B = BandedBlockBandedMatrix{Int64}(undef, axes(A), blockbandwidths(A), subblockbandwidths(A))

        @test eltype(A) === eltype(B) === Int64
        @test typeof(A.data) <: BlockArray
        @test typeof(B.data) <: PseudoBlockArray
        @test size(A) == size(B)
        @test bandrange(A) == bandrange(B)
        @test blockbandwidths(A) == blockbandwidths(B)
        @test blocksize(A) == blocksize(B)

        A = WithBlockMatrix{Int64}(Zeros{Int64}(sum(args[1]),sum(args[2])), args...)
        B = BandedBlockBandedMatrix{Int64}(Zeros{Int64}(sum(args[1]),sum(args[2])), args...)
        @test typeof(A.data) <: BlockArray
        @test typeof(B.data) <: PseudoBlockArray
        @test A == B

        A = WithBlockMatrix{Int64}(Ones{Int64}(sum(args[1]),sum(args[2])), args...)
        B = BandedBlockBandedMatrix{Int64}(Ones{Int64}(sum(args[1]),sum(args[2])), args...)
        @test typeof(A.data) <: BlockArray
        @test typeof(B.data) <: PseudoBlockArray
        @test A == B
        @test (A .+ 1) .* 2 == B .* 2 .+ 2
    end

    @testset "Sparse dimensions" begin
        A = BandedBlockBandedMatrix{Float64}(undef, 1:5,1:5, (-1,1), (-1,1))
        @test size(sparse(A)) == size(A) == (15,15)
    end

    @testset "non-standard blocks" begin
        A = BandedBlockBandedMatrix{Float64}(undef, Int[], 1:5,(-1,1), (-1,1))
        @test BlockBandedMatrices.colstart(A,1) == 1
        A = BandedBlockBandedMatrix{Float64}(undef, 1:2, 1:5,(-1,1), (-1,1))
        A.data .= randn.()
        V = view(A, Block(2,3))
        @test MemoryLayout(typeof(V)) == BandedMatrices.BandedColumnMajor()
        @test isbanded(V)
        @test bandwidths(V) == (-1,1)
        @test BandedMatrix(V) == A[2:3,4:6]
        V2 = view(V, :, 2:3)
        @test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumnMajor()
        @test bandwidths(V2) == (0,0)
        @test BandedMatrix(V2) == A[2:3,5:6]
        V2 = view(V, 2:2, :)
        @test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumnMajor()
        @test bandwidths(V2) == (-2,2)
        @test BandedMatrix(V2) == A[3:3,4:6]
        V2 = view(V, 2:2, 2:3)
        @test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumnMajor()
        @test bandwidths(V2) == (-1,1)
        @test BandedMatrix(V2) == A[3:3,5:6]
    end

    @testset "OneTo axes" begin
        ret = BandedBlockBandedMatrix{Int}(undef, (Base.OneTo(5),Base.OneTo(5)), (0,0), (1,1))
        V = view(ret, Block(1,1))
        @test bandwidths(V) == (1,1)
        A = brand(5,5,1,1)
        B = BandedBlockBandedMatrix(A)
        @test A == B
        @test subblockbandwidths(B) == (1,1)
        @test blockbandwidths(B) == (0,0)
    end
end

if false # turned off since tests have check-bounds=yes
    # test that @inbounds is working properly
    exceed_band(V, k, j) = @inbounds return V[k,j]
    @test exceed_band(V, 2,1) == 8

    @test BandedMatrices.inbands_getindex(V, 2, 1) == 8
    BandedMatrices.inbands_setindex!(V, -2, 5, 1)
    @test A[2,1] == -2
end


text:
export CuSubArray, CuSubVector, CuSubMatrix, CuSubVecOrMat

immutable CuSubArray{T,N,I}
    parent::CuArray{T}
    indexes::I
    offset1::Int
end

#Base.Array(x::CuSubArray) = x.parent[x.indexes...]
#Base.SubArray(x::CuSubArray) = view(Array(x.parent), x.indexes)

function Base.size{T,N,I}(x::CuSubArray{T,N,I}, dim::Int)
    typeof(x.indexes[dim]) == Colon ? size(x.parent,dim) : length(x.indexes[dim])
end
Base.size{T}(x::CuSubArray{T,1}) = (size(x,1),)
Base.size{T}(x::CuSubArray{T,2}) = (size(x,1),size(x,2))
Base.size{T}(x::CuSubArray{T,3}) = (size(x,1),size(x,2),size(x,3))
Base.size{T}(x::CuSubArray{T,4}) = (size(x,1),size(x,2),size(x,3),size(x,4))

Base.strides(x::CuSubArray) = strides(x.parent)
Base.stride(x::CuSubArray, dim::Int) = stride(x.parent, dim)

Base.length(x::CuSubArray) = prod(size(x))
Base.similar{T,N}(x::CuSubArray{T}, dims::NTuple{N,Int}) = CuArray{T}(dims)

function Base.pointer(x::CuSubArray, index::Int=1)
    index == 1 || throw("Not implemented yet.")
    pointer(x.parent, x.offset1+index)
end

function Base.view{T,N,I<:Tuple}(parent::CuArray{T,N}, indexes::I)
    @assert N == length(indexes)
    n = N
    offset1 = 0
    for i = 1:length(indexes)
        index = indexes[i]
        typeof(index) == Int && (n -= 1)
        f = typeof(index) == Colon ? 0 : first(index)-1
        offset1 += f * stride(parent,i)
    end
    CuSubArray{T,n,I}(parent, indexes, offset1)
end
Base.view(parent::CuArray, indexes::Union{UnitRange{Int},Colon,Int}...) = view(parent,indexes)

ext: jl


text:
using MultivariateStats
using LinearAlgebra
using Test
using StableRNGs
using Statistics: mean, cov, cor

@testset "CCA" begin

    rng = StableRNG(34568)

    dx = 5
    dy = 6
    p = 3

    # CCA with zero means

    X = rand(rng, dx, 100)
    Y = rand(rng, dy, 100)

    Px = qr(randn(rng, dx, dx)).Q[:, 1:p]
    Py = qr(randn(rng, dy, dy)).Q[:, 1:p]

    M = CCA(Float64[], Float64[], Px, Py, [0.8, 0.6, 0.4])

    @test size(M)[1] == dx
    @test size(M)[2] == dy
    @test mean(M, :x) == zeros(dx)
    @test mean(M, :y) == zeros(dy)
    @test_throws ArgumentError mean(M, :z)
    @test projection(M, :x) == Px
    @test projection(M, :y) == Py
    @test_throws ArgumentError projection(M, :z)
    @test cor(M) == [0.8, 0.6, 0.4]

    @test predict(M, X, :x) ≈ Px'X
    @test predict(M, Y, :y) ≈ Py'Y

    ## CCA with nonzero means

    ux = randn(rng, dx)
    uy = randn(rng, dy)

    M = CCA(ux, uy, Px, Py, [0.8, 0.6, 0.4])

    @test size(M)[1] == dx
    @test size(M)[2] == dy
    @test mean(M, :x) == ux
    @test mean(M, :y) == uy
    @test projection(M, :x) == Px
    @test projection(M, :y) == Py
    @test cor(M) == [0.8, 0.6, 0.4]

    @test predict(M, X, :x) ≈ Px' * (X .- ux)
    @test predict(M, Y, :y) ≈ Py' * (Y .- uy)

    ## prepare data

    n = 1000
    dg = 10
    G = randn(rng, dg, n)

    X = randn(rng, dx, dg) * G + 0.2 * randn(rng, dx, n)
    Y = randn(rng, dy, dg) * G + 0.2 * randn(rng, dy, n)
    xm = vec(mean(X, dims=2))
    ym = vec(mean(Y, dims=2))
    Zx = X .- xm
    Zy = Y .- ym

    Cxx = cov(X, dims=2)
    Cyy = cov(Y, dims=2)
    Cxy = cov(X, Y, dims=2)
    Cyx = Cxy'

    ## ccacov

    # X ~ Y
    M = fit(CCA, X, Y; method=:cov, outdim=p)
    Px = projection(M, :x)
    Py = projection(M, :y)
    rho = cor(M)
    @test size(M)[1] == dx
    @test size(M)[2] == dy
    @test size(M)[3] == p
    @test mean(M, :x) == xm
    @test mean(M, :y) == ym
    @test issorted(rho; rev=true)

    @test Px' * Cxx * Px ≈ Matrix(I, p, p)
    @test Py' * Cyy * Py ≈ Matrix(I, p, p)
    @test Cxy * (Cyy \ Cyx) * Px ≈ Cxx * Px * Diagonal(rho.^2)
    @test Cyx * (Cxx \ Cxy) * Py ≈ Cyy * Py * Diagonal(rho.^2)
    @test Px ≈ MultivariateStats.qnormalize!(Cxx \ (Cxy * Py), Cxx)
    @test Py ≈ MultivariateStats.qnormalize!(Cyy \ (Cyx * Px), Cyy)

    # Y ~ X
    M = fit(CCA, Y, X; method=:cov, outdim=p)
    Py = projection(M, :x)
    Px = projection(M, :y)
    rho = cor(M)
    @test size(M)[1] == dy
    @test size(M)[2] == dx
    @test size(M)[3] == p
    @test mean(M, :x) == ym
    @test mean(M, :y) == xm
    @test issorted(rho; rev=true)

    @test Px' * Cxx * Px ≈ Matrix(I, p, p)
    @test Py' * Cyy * Py ≈ Matrix(I, p, p)
    @test Cxy * (Cyy \ Cyx) * Px ≈ Cxx * Px * Diagonal(rho.^2)
    @test Cyx * (Cxx \ Cxy) * Py ≈ Cyy * Py * Diagonal(rho.^2)
    @test Px ≈ MultivariateStats.qnormalize!(Cxx \ (Cxy * Py), Cxx)
    @test Py ≈ MultivariateStats.qnormalize!(Cyy \ (Cyx * Px), Cyy)

    ## ccasvd

    # n > d
    M = fit(CCA, X, Y; method=:svd, outdim=p)
    Px = projection(M, :x)
    Py = projection(M, :y)
    rho = cor(M)
    @test size(M)[1] == dx
    @test size(M)[2] == dy
    @test size(M)[3] == p
    @test mean(M, :x) == xm
    @test mean(M, :y) == ym
    @test issorted(rho; rev=true)

    @test Px' * Cxx * Px ≈ Matrix(I, p, p)
    @test Py' * Cyy * Py ≈ Matrix(I, p, p)
    @test Cxy * (Cyy \ Cyx) * Px ≈ Cxx * Px * Diagonal(rho.^2)
    @test Cyx * (Cxx \ Cxy) * Py ≈ Cyy * Py * Diagonal(rho.^2)
    @test Px ≈ MultivariateStats.qnormalize!(Cxx \ (Cxy * Py), Cxx)
    @test Py ≈ MultivariateStats.qnormalize!(Cyy \ (Cyx * Px), Cyy)

    # different input types
    XX = convert.(Float32, X)
    YY = convert.(Float32, Y)

    MM = fit(CCA, view(XX, :, 1:400), view(YY, :, 1:400); method=:svd, outdim=p)

    # test that mixing types doesn't error
    predict(M, XX, :x)
    predict(M, YY, :y)
    predict(MM, XX, :x)
    predict(MM, YY, :y)

    # type stability
    for func in (M->mean(M, :x), M->mean(M, :y),
                 M->projection(M, :x),
                 M->projection(M, :y), cor)
        @test eltype(func(M)) == Float64
        @test eltype(func(MM)) == Float32
    end
end


text:
# Model structure
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

const IntTuple = Union{Tuple{Integer,Integer}, Tuple{Integer,Integer,Integer},Array{Int64,1},Array{Int32,1}}
const RealTuple = Union{Tuple{AbstractFloat,AbstractFloat}, Tuple{AbstractFloat,AbstractFloat,AbstractFloat},Array{AbstractFloat,1}}

export Model, PhysicalParameter, get_dt

###################################################################################################
# PhysicalParameter abstract vector

"""
    PhysicalParameter
        n::IntTuple
        d::RealTuple
        o::RealTuple
        data::Union{Array, Number}

PhysicalParameter structure for physical space parameter.

`n`: number of gridpoints in (x,y,z) for 3D or (x,z) for 2D

`d`: grid spacing in (x,y,z) or (x,z) (in meters)

`o`: origin of coordinate system in (x,y,z) or (x,z) (in meters)

`data`: the array of the parameter values of size n

Constructor
===========

A `PhysicalParameter` can be constructed in various ways but always require the origin `o` and grid spacing `d` that
cannot be infered from the array.

    PhysicalParameter(v::Array{vDT}, d, o) where `v` is an n-dimensional array and n=size(v)

    PhysicalParameter(n, d, o; vDT=Float32) Creates a zero PhysicalParameter

    PhysicalParameter(v::Array{vDT}, A::PhysicalParameter) Creates a PhysicalParameter from the Array `v` with n, d, o from `A`

    PhysicalParameter(v::Array{vDT, N}, n::Tuple, d::Tuple, o::Tuple) where `v` is a vector or nd-array that is reshaped into shape `n`

    PhysicalParameter(v::vDT, n::Tuple, d::Tuple, o::Tuple) Creates a constant (single number) PhyicalParameter

"""
mutable struct PhysicalParameter{vDT} <: AbstractVector{vDT}
    n::Tuple
    d::Tuple
    o::Tuple
    data::Union{Array{vDT}, vDT}
end

mutable struct PhysicalParameterException <: Exception
    msg :: String
end

PhysicalParameter(v::BitArray{N}, args...) where N = v

function PhysicalParameter(v::Array{vDT}, d::Tuple, o::Tuple) where {vDT}
    length(size(v)) != length(o) && throw(PhysicalParameterException("Input array should be $(length(o))-dimensional"))
    return PhysicalParameter{vDT}(size(v), d, o, v)
end

function PhysicalParameter(n::Tuple, d::Tuple, o::Tuple; vDT=Float32)
    return PhysicalParameter{vDT}(n, d, o, zeros(vDT, n))
end

function PhysicalParameter(v::Array{vDT}, A::PhysicalParameter) where {vDT}
    return PhysicalParameter{vDT}(A.n, A.d, A.o, reshape(v, A.n))
end

function PhysicalParameter(v::Array{vDT, N}, n::Tuple, d::Tuple, o::Tuple) where {vDT, N}
    length(v) != prod(n) && throw(PhysicalParameterException("Incompatible number of element in input $(length(v)) with n=$(n)"))
    N == 1 && (v = reshape(v, n))
    return PhysicalParameter{vDT}(n, d, o, v)
end

function PhysicalParameter(v::vDT, n::Tuple, d::Tuple, o::Tuple) where vDT<:Number
    return PhysicalParameter{vDT}(n, d, o, v)
end

PhysicalParameter(p::PhysicalParameter, n::Tuple, d::Tuple, o::Tuple) = p
PhysicalParameter(p::PhysicalParameter) = p

# transpose and such.....
conj(x::PhysicalParameter{vDT}) where vDT = x
transpose(x::PhysicalParameter{vDT}) where vDT = PhysicalParameter{vDT}(x.n[end:-1:1], x.d[end:-1:1], x.o[end:-1:1], permutedims(x.data, length(x.n):-1:1))
adjoint(x::PhysicalParameter{vDT}) where vDT = transpose(x)

# Basic overloads
size(A::PhysicalParameter) = (prod(A.n), 1)
length(A::PhysicalParameter) = prod(A.n)

function norm(A::PhysicalParameter, order::Real=2)
    return norm(vec(A.data), order)
end

dot(A::PhysicalParameter, B::PhysicalParameter) = dot(vec(A.data), vec(B.data))
dot(A::PhysicalParameter, B::Array) = dot(vec(A.data), vec(B))
dot(A::Array, B::PhysicalParameter) = dot(vec(A), vec(B.data))

display(A::PhysicalParameter) = println("$(typeof(A)) of size $(A.n) with origin $(A.o) and spacing $(A.d)")
show(io::IO, A::PhysicalParameter) = print(io, "$(typeof(A)) of size $(A.n) with origin $(A.o) and spacing $(A.d)")
summary(io::IO, A::PhysicalParameter) = print(io, "$(typeof(A)) of size $(A.n) with origin $(A.o) and spacing $(A.d)")
showarg(io::IO, A::PhysicalParameter, toplevel) = print(io, typeof(A), " with size $(A.n), spacing $(A.d) and origin $(A.o)")
show(io::IO, ::MIME{Symbol("text/plain")}, A::PhysicalParameter) = println(io, "$(typeof(A)) of size $(A.n) with origin $(A.o) and spacing $(A.d)")

# Indexing
firstindex(A::PhysicalParameter) = 1
lastindex(A::PhysicalParameter) = length(A)
lastindex(A::PhysicalParameter, dim::Int) = A.n[dim]

function promote_shape(p::PhysicalParameter, A::Array{vDT, N}) where {vDT, N}
    (size(A) != p.n && N>1) && return promote_shape(p.data, A)
    (length(A) == prod(p.n) && N==1) && return size(A)
    return promote_shape(A, A)
end

promote_shape(A::Array{vDT, N}, p::PhysicalParameter) where {vDT, N} = promote_shape(p, A)
reshape(p::PhysicalParameter, n::Tuple{Vararg{Int64,N}}) where N = (prod(n)==prod(p.n) && return p)

dotview(A::PhysicalParameter{vDT}, I::Vararg{Union{Function, Int, UnitRange{Int}}, N}) where {vDT, N} = dotview(A.data, I...)
Base.dotview(m::PhysicalParameter, i) = Base.dotview(m.data, i)

getindex(A::PhysicalParameter, i::Int) = A.data[i]
getindex(A::PhysicalParameter, i::Colon) = A.data[:]

get_step(r::StepRange) = r.step
get_step(r) = 1

function getindex(A::PhysicalParameter{T}, I::Vararg{Union{Int, BitArray, Function, StepRange{Int}, UnitRange{Int}}, N}) where {N, T}
    new_v = getindex(A.data, I...)
    length(size(new_v)) != length(A.n) && (return new_v)
    s = [i == (:) ? 0 : i[1]-1 for i=I]
    st = [get_step(i) for i=I]
    new_o = [ao+i*d for (ao, i, d)=zip(A.o, s, A.d)]
    new_d = [d*s for (d, s)=zip(A.d, st)]
    PhysicalParameter{T}(size(new_v), tuple(new_d...), tuple(new_o...), new_v)
end

setindex!(A::PhysicalParameter, v, I::Vararg{Union{Int, Function, UnitRange{Int}}, N}) where {N} = setindex!(A.data, v, I...)
setindex!(A::PhysicalParameter, v, i::Int) = (A.data[i] = v)

# Constructiors by copy
similar(x::PhysicalParameter{vDT}) where {vDT} = vDT(0) .* x
copy(x::PhysicalParameter{vDT}) where {vDT} = PhysicalParameter{vDT}(x.n, x.d, x.o, x.data)
unsafe_convert(::Type{Ptr{T}}, p::PhysicalParameter{T}) where {T} = unsafe_convert(Ptr{T}, p.data)

# Equality
==(A::PhysicalParameter, B::PhysicalParameter) = (A.data == B.data && A.o == B.o && A.d == B.d)
isapprox(A::PhysicalParameter, B::PhysicalParameter; kwargs...) = (isapprox(A.data, B.data) && A.o == B.o && A.d == B.d)
isapprox(A::PhysicalParameter, B::AbstractArray; kwargs...) = isapprox(A.data, B)
isapprox(A::AbstractArray, B::PhysicalParameter; kwargs...) = isapprox(A, B.data)

# # Arithmetic operations
function compare(A::PhysicalParameter, B::PhysicalParameter)
    A.o != B.o && throw(PhysicalParameterException("Incompatible origins $(A.o) and $(B.o)"))
    A.d != B.d && throw(PhysicalParameterException("Incompatible spacing $(A.d) and $(B.d)"))
    A.n != B.n && throw(PhysicalParameterException("Incompatible sizes $(A.n) and $(B.n)"))
end

for op in [:+, :-, :*, :/]
    @eval function $(op)(A::PhysicalParameter{T}, B::PhysicalParameter{T}) where T
        compare(A, B)
        return PhysicalParameter(broadcast($(op), A.data, B.data), A)
    end
end

# Brodacsting
BroadcastStyle(::Type{<:PhysicalParameter}) = Broadcast.ArrayStyle{PhysicalParameter}()

function similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{PhysicalParameter}}, ::Type{ElType}) where ElType
    # Scan the inputs for the ArrayAndChar:
    A = find_pm(bc)
    # Use the char field of A to create the output
    ElType == Bool ? Ad = similar(Array{ElType}, axes(A.data)) : Ad = similar(A.data)
    PhysicalParameter(Ad, A.d, A.o)
end

"`A = find_pm(As)` returns the first PhysicalParameter among the arguments."
find_pm(bc::Base.Broadcast.Broadcasted) = find_pm(bc.args)
find_pm(args::Tuple) = find_pm(find_pm(args[1]), Base.tail(args))
find_pm(x) = x
find_pm(::Tuple{}) = nothing
find_pm(a::PhysicalParameter, rest) = a
find_pm(::Any, rest) = find_pm(rest)

for op in [:+, :-, :*, :/]
    for (T1, T2) in ([DenseArray, PhysicalParameter], [PhysicalParameter, DenseArray],
                     [PhysicalParameter, PhysicalParameter])
        @eval function broadcasted(::typeof($op), A::$T1, B::$T2)
            pm = find_pm(A, B)
            return PhysicalParameter(materialize(broadcasted($(op), A[:], B[:])), pm)
        end
    end
    @eval broadcasted(::typeof($op), p::PhysicalParameter, bc::Base.Broadcast.Broadcasted) = broadcasted($(op), p, materialize(bc))
    @eval broadcasted(::typeof($op), bc::Base.Broadcast.Broadcasted, p::PhysicalParameter) = broadcasted($(op), materialize(bc), p)
end

function *(A::Union{joMatrix, joLinearFunction, joLinearOperator, joCoreBlock}, p::PhysicalParameter{RDT}) where {RDT}
    @warn "JOLI linear operator, returning julia Array"
    return A*vec(p.data)
end

materialize!(p::PhysicalParameter{RDT}, b::Array{<:Number, N}) where{RDT, N} = (p.data[:] .= b)
materialize!(p::PhysicalParameter{RDT}, b::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{N}}) where{RDT, N} = materialize!(p, collect(b))
materialize!(p::Array, b::Broadcast.Broadcasted{Broadcast.ArrayStyle{PhysicalParameter}}) = materialize!(p, reshape(collect(b), size(p)))

# Linalg Extras
mul!(x::PhysicalParameter, F::Union{joAbstractLinearOperator, joLinearFunction, Array}, y::Array) = mul!(view(x.data, :), F, y)
mul!(x::PhysicalParameter, F::Union{joAbstractLinearOperator, joLinearFunction, Array}, y::PhysicalParameter) = mul!(view(x.data, :), F, view(y.data, :))
mul!(x::Array, F::Union{joAbstractLinearOperator, joLinearFunction, Array}, y::PhysicalParameter) = mul!(x, F, y[1:end])

# For ploting
array2py(p::PhysicalParameter{vDT}, i::Int64, I::CartesianIndex{N}) where {vDT, N} = array2py(p.data, i, I)

###################################################################################################
# Isotropic acoustic

"""
    Model
        n::IntTuple
        d::RealTuple
        o::RealTuple
        nb::Integer
        params::Dict
        rho::Array

Model structure for seismic velocity models.

`n`: number of gridpoints in (x,y,z) for 3D or (x,z) for 2D

`d`: grid spacing in (x,y,z) or (x,z) (in meters)

`o`: origin of coordinate system in (x,y,z) or (x,z) (in meters)

`nb`: number of absorbing boundary points in each direction

`params`: Physical parameters such has squared slowness, denisty or THomesne parameters

Constructor
===========

The parameters `n`, `d`, `o` and `m` are mandatory, whith `nb` and other physical parameters being optional input arguments.

    Model(n, d, o, m; nb=40, rho=1, epsilon=0, delta=0, theta=0, phi=0)

where

`m`: velocity model in slowness squared (s^2/km^2)

`epsilon`: Epsilon thomsen parameter ( between -1 and 1)

`delta`: Delta thomsen parameter ( between -1 and 1 and delta < epsilon)

`theta`: Anisotopy dip in radian

`phi`: Anisotropy asymuth in radian

`rho`: density (g / m^3)

"""
mutable struct Model
    n::IntTuple
    d::RealTuple
    o::RealTuple
    nb::Integer # number of absorbing boundaries points on each side
    params::Dict
end

###################################################################################################
# Constructors

function Model(n::IntTuple, d::RealTuple, o::RealTuple, m;
               epsilon=nothing, delta=nothing, theta=nothing,
               phi=nothing, rho=nothing, qp=nothing, nb=40)

    params = Dict(:m => PhysicalParameter(Float32.(m), d, o))
    for (name, val)=zip([:rho, :qp, :epsilon, :delta, :theta, :phi], [rho, qp, epsilon, delta, theta, phi])
        ~isnothing(val) && (params[name] = PhysicalParameter(Float32.(val), n, d, o))
    end

    return Model(n, d, o, nb, params)
end

function Model(n::IntTuple, d::RealTuple, o::RealTuple, m, rho; nb=40)
    return Model(n, d, o, m; rho=rho, nb=nb)
end

function Model(n::IntTuple, d::RealTuple, o::RealTuple, m, rho, qp; nb=40)
    return Model(n, d, o, m; rho=rho, qp=qp, nb=nb)
end

get_dt(m::Model; dt=nothing) = calculate_dt(m; dt=dt)
getindex(m::Model, sym::Symbol) = m.params[sym]

Base.setproperty!(m::Model, s::Symbol, p::PhysicalParameter{Float32}) = (m.params[s] = p)

function Base.getproperty(obj::Model, sym::Symbol)
    if sym == :params
        return getfield(obj, sym)
    elseif sym in keys(obj.params)
        return obj.params[sym]
    else
         # fallback to getfield
        return getfield(obj, sym)
    end
end

similar(::PhysicalParameter{vDT}, m::Model) where {vDT} = PhysicalParameter(m.n, m.d, m.o; vDT=vDT)
similar(x::Array, m::Model) = similar(x, m.n)

ndims(m::Model) = ndims(m.m.data)

display(m::Model) = println("Model (n=$(m.n), d=$(m.d), o=$(m.o)) with parameters $(keys(m.params))")
show(io::IO, m::Model) = print(io, "Model (n=$(m.n), d=$(m.d), o=$(m.o)) with parameters $(keys(m.params))")
show(io::IO, ::MIME{Symbol("text/plain")}, m::Model) = print(io, "Model (n=$(m.n), d=$(m.d), o=$(m.o)) with parameters $(keys(m.params))")


text:
export params, Linear

abstract type AbstractLayer <: Func end

function (lyr::AbstractLayer)(args...)
    outputs = forward(lyr, args...)
    outputs isa Tuple || (outputs = (outputs,))
    lyr.args = Tuple(WeakRef(x) for x in args)
    lyr.outputs = Tuple(WeakRef(y) for y in outputs)
    return length(outputs) > 1 ? outputs : outputs[1]
end

function params(lyr::AbstractLayer)
    pnames = propertynames(lyr)
    ps = []
    for pname in pnames
        pname in (:args, :outputs) && continue
        p = isdefined(lyr, pname) ? getfield(lyr, pname) : nothing
        if isparam(p)
            push!(ps, p)
        elseif p isa AbstractLayer
            append!(ps, params(p))
        end
    end
    return ps
end

function cleargrad!(x::AbstractLayer)
    for p in params(x)
        cleargrad!(p)
    end
end

cleargrad!(x::AbstractLayer...) = cleargrad!.(x)

"""
    Linear <: AbstractLayer
"""
@func mutable struct Linear <: AbstractLayer
    weight::Var
    bias::Union{Var,Nothing}

    function Linear(in, out, nobias=false, dtype=Float32)
        w = Param(_initW(in, out, dtype))
        b = nobias ? nothing : Param(zeros(dtype, out))
        return new(w, b)
    end
end

_initW(in, out, dtype) = randn(dtype, out, in) .* dtype(sqrt(1 / in))

function forward(lyr::Linear, x)
    return linear(lyr.weight, x, lyr.bias)
end

function Base.show(io::IO, l::Linear)
  print(io, "Linear(", size(l.weight, 2), ", ", size(l.weight, 1))
  l.bias isa Nothing && print(io, "; bias=false")
  print(io, ")")
end


text:
export stoichiometriccoeffs,
    stoichiometricsources,
    stoichiometrictargets,
    kineticorder,
    stoichiometricmatrix,
    conservativelaws,
    Velocity,
    Velocities,
    columns,
    VectorVelocities,
    SpeciesformationTotalRing,
    SpeciesformationKsRing,
    SpeciesformationRings,
    Speciesformationrate,
    SpeciesformationrateInStclass!,
    SpeciesformationrateInStclass,
    SystemF

function stoichiometricsources(net::AbstractMatrix, xs::AbstractVector, getcoefficient)
    return Int.([getcoefficient(p, x) for x in xs, p in net[:,1]])
end

function stoichiometrictargets(net::AbstractMatrix, xs::AbstractVector, getcoefficient)
    return Int.([getcoefficient(p, x) for x in xs, p in net[:,2]])
end

function stoichiometriccoeffs(net::AbstractMatrix, xs::AbstractVector, getcoefficient)
    return hcat(stoichiometricsources(net, xs, getcoefficient), stoichiometrictargets(net, xs, getcoefficient))
end

function stoichiometriccoeffs(net::AbstractMatrix{T}, xs::AbstractVector) where {T <: RingElem}
    return hcat(stoichiometricsources(net, xs, Nemo.coeff), stoichiometrictargets(net, xs, Nemo.coeff))
end

function stoichiometriccoeffs(net::AbstractMatrix, xs::AbstractVector)
    R = parent(xs[1])
    return stoichiometriccoeffs(R.(net), xs)
end

function stoichiometricsources(net::AbstractMatrix)
    return net[:,1:div(end, 2)]
end

function stoichiometrictargets(net::AbstractMatrix)
    return net[:,div(end, 2) + 1:end]
end

# Matrix Y
function kineticorder(net::AbstractMatrix)
    return dropzeroslices(stoichiometricsources(net))
end

# Matrix N
function stoichiometricmatrix(net::AbstractMatrix)
    return dropzeroslices(stoichiometrictargets(net) - stoichiometricsources(net))
end

function conservativelaws(N::AbstractMatrix{T}) where {T <: Integer}
    Nnemo = matrix(FlintIntegerRing(), N)
    nTs, W = left_kernel(Nnemo)
    W = hnf(W)
    return nTs, T.(Array(W))
end

##############################################
##############################################
##############################################
############ Building matrices in Nemo for CRN

################### With k and x
## Velocities
## Input: Array with 'x' variables xs, and the kinetics matrix
## Output: Vector of velocities of the reaction network
## VectorVelocities(gens(R), Y)
function Velocity(R::MPolyRing, v::AbstractVector)
    newV = MPolyBuildCtx(R)
    push_term!(newV, one(R), v)
    return finish(newV)
end

function Velocities(R::MPolyRing)
    return v -> Velocity(R, v)
end

function columns(A::AbstractMatrix)
    return [A[:, i] for i in 1:size(A, 2)]
end

function VectorVelocities(R::MPolyRing, Y::AbstractMatrix)
    return Velocities(R).(columns(Y))
end

function SpeciesformationTotalRing(net, nxs::Integer)
    nks = size(net, 1)
    return PolynomialRing(FlintRationalField(), vcat(["x$i" for i in 1:nxs], ["k$i" for i in 1:nks]))
end

function SpeciesformationKsRing(net)
    nks = size(net, 1)
    return PolynomialRing(FlintRationalField(), ["k$i" for i in 1:nks])
end

function SpeciesformationRings(net, nxs::Integer)
    return (SpeciesformationKsRing(net)..., PolynomialRing(SpeciesformationKsRing(net)[1], ["x$i" for i in 1:nxs])...)
end

function Speciesformationrate(Rx::MPolyRing, Rk::MPolyRing,
                              N::AbstractMatrix, Y::AbstractMatrix)
    V = VectorVelocities(Rx, Y)
    K = Diagonal(gens(Rk))
    return N * K * V
end

function SpeciesformationrateInStclass!(f::AbstractVector,
                                        R::MPolyRing,
                                        W::AbstractMatrix,
                                        xs::UnitRange{T}=1:size(W, 2),
                                        ts::UnitRange{T}=.+(length(xs), 1:size(W, 1))) where {T <: Integer}
    Wx = W * gens(R)[xs] - gens(R)[ts]
    for (i, p) in enumerate(findpivots(W))
        f[p] = Wx[i]
    end
    return f
end

function SpeciesformationrateInStclass(f::AbstractVector,
                                       R::MPolyRing,
                                       W::AbstractMatrix,
                                       xs::UnitRange{T}=1:size(W, 2),
                                       ts::UnitRange{T}=.+(length(xs), 1:size(W, 1))) where {T <: Integer}
    F = f
    return SpeciesformationrateInStclass!(F, R, W, xs, ts)
end

function SystemF(net, nxs)
    Y = StoichiometricCoeffs(net, nxs)
    Ys = StoichiometricSources(Y)
    N = StoichiometricMatrix(Y)
end


text:
"""
    PanelProperties

Panel specific properties calculated during the vortex lattice method analysis.

**Fields**
 - `gamma`: Vortex ring circulation strength, normalized by the reference velocity
 - `velocity`: Local velocity at the panel's bound vortex center, normalized by
    the reference velocity
 - `cfb`: Net force on the panel's bound vortex, as calculated using the
    Kutta-Joukowski theorem, normalized by the reference dynamic pressure and area
 - `cfl`: Force on the left bound vortex from this panel's vortex ring, as
    calculated by the Kutta-Joukowski theorem, normalized by the reference
    dynamic pressure and area
 - `cfr`: Force on the right bound vortex from this panel's vortex ring, as
    calculated by the Kutta-Joukowski theorem, normalized by the reference
    dynamic pressure and area
"""
struct PanelProperties{TF}
    gamma::TF
    velocity::SVector{3, TF}
    cfb::SVector{3, TF}
    cfl::SVector{3, TF}
    cfr::SVector{3, TF}
end

# constructor
function PanelProperties(gamma, velocity, cfb, cfl, cfr)

    TF = promote_type(typeof(gamma), typeof(velocity), eltype(cfb), eltype(cfl),
        eltype(cfr))

    return PanelProperties{TF}(gamma, velocity, cfb, cfl, cfr)
end

Base.eltype(::Type{PanelProperties{TF}}) where TF = TF
Base.eltype(::PanelProperties{TF}) where TF = TF

"""
    System{TF}

Contains pre-allocated storage for internal system variables.

# Fields:
 - `AIC`: Aerodynamic influence coefficient matrix from the surface panels
 - `w`: Normal velocity at the control points from external sources and wakes
 - `Γ`: Circulation strength of the surface panels
 - `V`: Velocity at the wake vertices for each surface
 - `surfaces`: Surfaces, represented by matrices of surface panels
 - `properties`: Surface panel properties for each surface
 - `wakes`: Wake panel properties for each surface
 - `trefftz`: Trefftz panels associated with each surface
 - `reference`: Pointer to reference parameters associated with the system (see [`Reference`](@ref))
 - `freestream`: Pointer to current freestream parameters associated with the system (see [`Freestream`](@ref))
 - `symmetric`: Flags indicating whether each surface is symmetric across the X-Z plane
 - `nwake`: Number of chordwise wake panels to use from each wake in `wakes`,
 - `surface_id`: Surface ID for each surface.  The finite core model is disabled
    when calculating the influence of surfaces/wakes that share the same ID.
 - `trailing_vortices`: Flags to enable/disable trailing vortices
 - `wake_finite_core`: Flag for each wake indicating whether the finite core
    model should be enabled when calculating the wake's influence on itself and
    surfaces/wakes with the same surface ID.
 - `near_field_analysis`: Flag indicating whether a near field analysis has been
    performed for the current system state
 - `derivatives`: Flag indicating whether the derivatives with respect to the
    freestream variables have been calculated
 - `dw`: Derivatives of the R.H.S. with respect to the freestream variables
 - `dΓ`: Derivatives of the circulation strength with respect to the freestream variables
 - `dproperties`: Derivatives of the panel properties with respect to the freestream variables
 - `wake_shedding_locations`: Wake shedding locations for each surface
 - `Vcp`: Velocity due to surface motion at the control points
 - `Vh`: Velocity due to surface motion at the horizontal bound vortex centers
 - `Vv`: Velocity due to surface motion at the vertical bound vortex centers
 - `Vte`: Velocity due to surface motion at the trailing edge vertices
 - `dΓdt`: Derivative of the circulation strength with respect to non-dimensional time
"""
struct System{TF}
    AIC::Matrix{TF}
    w::Vector{TF}
    Γ::Vector{TF}
    V::Vector{Matrix{SVector{3,TF}}}
    surfaces::Vector{Matrix{SurfacePanel{TF}}}
    properties::Vector{Matrix{PanelProperties{TF}}}
    wakes::Vector{Matrix{WakePanel{TF}}}
    trefftz::Vector{Vector{TrefftzPanel{TF}}}
    reference::Array{Reference{TF}, 0}
    freestream::Array{Freestream{TF}, 0}
    symmetric::Vector{Bool}
    nwake::Vector{Int}
    surface_id::Vector{Int}
    wake_finite_core::Vector{Bool}
    trailing_vortices::Vector{Bool}
    xhat::Array{SVector{3, TF}, 0}
    near_field_analysis::Array{Bool, 0}
    derivatives::Array{Bool, 0}
    dw::NTuple{5, Vector{TF}}
    dΓ::NTuple{5, Vector{TF}}
    dproperties::NTuple{5, Vector{Matrix{PanelProperties{TF}}}}
    wake_shedding_locations::Vector{Vector{SVector{3,TF}}}
    previous_surfaces::Vector{Matrix{SurfacePanel{TF}}}
    Vcp::Vector{Matrix{SVector{3, TF}}}
    Vh::Vector{Matrix{SVector{3, TF}}}
    Vv::Vector{Matrix{SVector{3, TF}}}
    Vte::Vector{Vector{SVector{3, TF}}}
    dΓdt::Vector{TF}
end

@inline Base.eltype(::Type{System{TF}}) where TF = TF
@inline Base.eltype(::System{TF}) where TF = TF

"""
    System([TF], surfaces; kwargs...)

Return an object of type `System` with pre-allocated storage for internal system
variables

# Arguments:
 - `TF`: Floating point type, defaults to the floating point type used by `surface`
 - `surfaces`:
        Either:
         - One or more grids of shape (3, nc+1, ns+1) which represents lifting surfaces,
         or
         - One or more matrices of surface panels (see [`SurfacePanel`](@ref)) of
           shape (nc, ns)
        where `nc` is the number of chordwise panels and `ns` is the number of
        spanwise panels

# Keyword Arguments:
 - `nw`: Number of chordwise wake panels to initialize for each surface. Defaults to
    zero wake panels for each surface.
"""
System(args...; kwargs...)

# grid inputs, no provided type
function System(grids::AbstractVector{<:AbstractArray{TF, 3}}; kwargs...) where TF

    return System(TF, grids; kwargs...)
end

# grid inputs, provided type
function System(TF::Type, grids::AbstractVector{<:AbstractArray{<:Any, 3}}; kwargs...)

    nc = size.(grids, 2) .- 1
    ns = size.(grids, 3) .- 1

    return System(TF, nc, ns; kwargs...)
end

# surface inputs, no provided type
function System(surfaces::AbstractVector{<:AbstractMatrix{SurfacePanel{TF}}};
    kwargs...) where TF

    return System(TF, surfaces; kwargs...)
end

# surface inputs, provided type
function System(TF::Type, surfaces::AbstractVector{<:AbstractMatrix{<:SurfacePanel}};
    kwargs...)

    nc = size.(surfaces, 1)
    ns = size.(surfaces, 2)

    return System(TF, nc, ns; kwargs...)
end

"""
    System(TF, nc, ns; nw = zero(nc))

Return an object of type `System` with pre-allocated storage for internal system
variables

# Arguments:
 - `TF`: Floating point type.
 - `nc`: Number of chordwise panels on each surface.
 - `ns`: Number of spanwise panels on each surface.

# Keyword Arguments
 - `nw`: Number of chordwise wake panels for each surface. Defaults to zero wake
    panels on each surface
"""
function System(TF::Type, nc, ns; nw = zero(nc))

    @assert length(nc) == length(ns) == length(nw)

    # number of surfaces
    nsurf = length(nc)

    # number of surface panels
    N = sum(nc .* ns)

    AIC = zeros(TF, N, N)
    w = zeros(TF, N)
    Γ = zeros(TF, N)
    V = [fill((@SVector zeros(TF, 3)), nw[i]+1, ns[i]+1) for i = 1:nsurf]
    surfaces = [Matrix{SurfacePanel{TF}}(undef, nc[i], ns[i]) for i = 1:nsurf]
    properties = [Matrix{PanelProperties{TF}}(undef, nc[i], ns[i]) for i = 1:nsurf]
    wakes = [Matrix{WakePanel{TF}}(undef, nw[i], ns[i]) for i = 1:nsurf]
    trefftz = [Vector{TrefftzPanel{TF}}(undef, ns[i]) for i = 1:nsurf]
    reference = Array{Reference{TF}}(undef)
    freestream = Array{Freestream{TF}}(undef)
    symmetric = [false for i = 1:nsurf]
    nwake = [0 for i = 1:nsurf]
    surface_id = [i for i = 1:nsurf]
    wake_finite_core = [true for i = 1:nsurf]
    trailing_vortices = [false for i = 1:nsurf]
    xhat = fill(SVector{3,TF}(1, 0, 0))
    near_field_analysis = fill(false)
    derivatives = fill(false)
    dw = Tuple(zeros(TF, N) for i = 1:5)
    dΓ = Tuple(zeros(TF, N) for i = 1:5)
    dproperties = Tuple([Matrix{PanelProperties{TF}}(undef, nc[i], ns[i]) for i = 1:length(surfaces)] for j = 1:5)
    wake_shedding_locations = [fill((@SVector zeros(TF, 3)), ns[i]+1) for i = 1:nsurf]
    previous_surfaces = [Matrix{SurfacePanel{TF}}(undef, nc[i], ns[i]) for i = 1:nsurf]
    Vcp = [fill((@SVector zeros(TF, 3)), nc[i], ns[i]) for i = 1:nsurf]
    Vh = [fill((@SVector zeros(TF, 3)), nc[i]+1, ns[i]) for i = 1:nsurf]
    Vv = [fill((@SVector zeros(TF, 3)), nc[i], ns[i]+1) for i = 1:nsurf]
    Vte = [fill((@SVector zeros(TF, 3)), ns[i]+1) for i = 1:nsurf]
    dΓdt = zeros(TF, N)

    return System{TF}(AIC, w, Γ, V, surfaces, properties, wakes, trefftz,
        reference, freestream, symmetric, nwake, surface_id, wake_finite_core,
        trailing_vortices, xhat, near_field_analysis, derivatives,
        dw, dΓ, dproperties, wake_shedding_locations, previous_surfaces, Vcp, Vh,
        Vv, Vte, dΓdt)
end

"""
    get_surface_properties(system)

Return a vector of surface panel properties for each surface, stored as matrices
of panel properties (see [`PanelProperties`](@ref)) of shape (nc, ns) where `nc`
is the number of chordwise panels and `ns` is the number of spanwise panels
"""
get_surface_properties(system) = system.properties


text:
#Kernel functions for use with hacvariance estimator
abstract type KernelFunction end
struct KernelEpanechnikov <: KernelFunction ; end
struct KernelGaussian <: KernelFunction ; end
struct KernelUniform <: KernelFunction ; end
struct KernelBartlett <: KernelFunction ; end
const KERNEL_FUNCTION_DICT = Dict{Symbol,KernelFunction}(
	:epanechnikov => KernelEpanechnikov(),
	:gaussian => KernelGaussian(),
	:uniform => KernelUniform(),
	:bartlett => KernelBartlett()
)::Dict{Symbol,KernelFunction}
(get_kernel_type(kf::Tk)::Tk) where {Tk<:KernelFunction} = kf
(get_kernel_type(kf::Symbol)) = get_kernel_type(KERNEL_FUNCTION_DICT[kf])
(get_kernel_type(kf::String)) = get_kernel_type(Symbol(kf))
(get_kernel_type(kf)) = error("Input type $(typeof(kf)) not able to be converted to kernel function type: $(kf)")
kernel_eval(x::Number, bw::Number, kf::KernelEpanechnikov)::Float64 = (-1.0 <= x / bw <= 1.0) ? (0.75 / bw) * (1 - (x^2 / bw^2)) : 0.0
kernel_eval(x::Number, bw::Number, kf::KernelGaussian)::Float64 = (1 / bw) * (1 / sqrt(2*pi)) * exp(-1 * (x^2 / (2 * bw^2)))
kernel_eval(x::Number, bw::Number, kf::KernelUniform)::Float64 = (-1.0 <= x / bw <= 1.0) ? 1 / 2*bw : 0.0
kernel_eval(x::Number, bw::Number, kf::KernelBartlett)::Float64 = (-1.0 <= x / bw <= 1.0) ? (1 - abs(x / bw)) : 0.0

"""
	hacvariance{T<:Number}(x::AbstractVector{T} ; kf::Symbol=:epanechnikov, bw::Int=-1)::Tuple{Float64, Int}

Get the heteroskedasticity and autocorrelation consistent variance estimator of data vector x.
The function has the following keyword arguments:
	kf <- Kernel function used in estimator. Valid values are :epanechnikov, :gaussian, :uniform, :bartlett
	bw <- Bandwidth used in estimator. If <= -1, then estimate bandwidth using Politis (2003) "Adaptive Bandwidth Choice"
"""
function hacvariance(x::AbstractVector{<:Number}, kf::Tk, bw::Int)::Tuple{Float64, Int} where {Tk<:KernelFunction}
	if bw <= -1 ; (bw, v, xCov) = DependentBootstrap.bandwidth_politis_2003(x)
	else ; (v, xCov) = (var(x), Float64[]) ; end
	length(xCov) < bw && append!(xCov, autocov(x, length(xCov)+1:bw)) #Get any additional autocovariances that we might need
	kernelAdjTerm = 1 / kernel_eval(0.0, bw, kf) #Used to scale kernel functions that don't satisfy k(0) = 1
	for m = 1:bw
		v += 2 * kernelAdjTerm * kernel_eval(m, bw, kf) * xCov[m]
	end
	return (max(v, 0.0), bw)
end
hacvariance(x::AbstractVector{<:Number} ; kf=KernelEpanechnikov(), bw::Int=-1)::Tuple{Float64, Int} = hacvariance(x, get_kernel_type(kf), bw)


text:
function diag(r,i)
    while !isborder(r,HorizonSide(i % 4)) && !isborder(r,HorizonSide((i+1) % 4))
        move!(r,HorizonSide(i % 4))
        move!(r,HorizonSide((i+1) % 4))
        putmarker!(r)
    end
    while ismarker(r)
        move!(r,HorizonSide((i+2) % 4))
        move!(r,HorizonSide((i+3) % 4))
    end
end

function task_13(r)
    for i in 0:3
        diag(r,i)
    end
    putmarker!(r)
end

task_13(r)


text:
function _point_to_cartesian(layer::T, c::Point; side=:center) where {T <: SimpleSDMLayer}
    lonside, latside = true, true
    if side == :bottomleft
        lonside, latside = false, false
    end
    if side == :topright
        lonside, latside = true, true
    end
    lon = SimpleSDMLayers._match_longitude(layer, c[1]; lower=lonside)
    lat = SimpleSDMLayers._match_latitude(layer, c[2]; lower=latside)
    isnothing(lon) && return nothing
    isnothing(lat) && return nothing
    return CartesianIndex(lat, lon)
end

function _match_latitude(layer::T, lat::K; lower::Bool=true) where {T <: SimpleSDMLayer, K <: AbstractFloat}
    layer.bottom <= lat <= layer.top || return nothing
    lat == layer.bottom  && return 1
    lat == layer.top  && return size(layer, 1)
    relative = (lat - layer.bottom)/(layer.top - layer.bottom)
    fractional = relative * size(layer, 1)+1
    if lat in layer.bottom:2stride(layer,2):layer.top
        f = lower ? floor : ceil
        d = lower ? 1 : 0
        return min(f(Int64, fractional-d), size(layer, 1))
    else
        return min(floor(Int, fractional), size(layer, 1))
    end
end

function _match_longitude(layer::T, lon::K; lower::Bool=true) where {T <: SimpleSDMLayer, K <: AbstractFloat}
    layer.left <= lon <= layer.right || return nothing
    lon == layer.left  && return 1
    lon == layer.right  && return size(layer, 2)
    relative = (lon - layer.left)/(layer.right - layer.left)
    fractional = relative * size(layer, 2)+1
    if lon in layer.left:2stride(layer,1):layer.right
        f = lower ? floor : ceil
        d = lower ? 1 : 0
        return min(f(Int64, fractional-d), size(layer, 2))
    else
        return min(floor(Int, fractional), size(layer, 2))
    end
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()
module Nettle
using Compat

# Load libnettle from BinDeps
const depfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depfile)
    include(depfile)
else
    error("libnettle not properly installed. Please run Pkg.build(\"Nettle\")")
end

include( "hash_common.jl" )
include( "hash.jl" )
include( "hmac.jl" )
include( "cipher.jl" )

function get_libnettle_version()
    # Current version (3.1.1) of nettle doesn't provide a runtime API to query
    # the version. Using the existance of a symbol added in 3.0 to determine
    # which version is loaded at runtime. See #46 and
    # http://upstream-tracker.org/versions/nettle.html
    hdl = @compat Libdl.dlopen_e(nettle)
    @compat(Libdl.dlsym_e(hdl, :nettle_aes192_invert_key)) == C_NULL ? 2 : 3
end

function __init__()
    global const nettle_major_version = get_libnettle_version()
end

if VERSION >= v"0.4.0-dev+6521"
    include("precompile.jl")
    _precompile_()
end

end # module


text:
depths = parse.(Int, readlines("day1/depths.txt"))

p1 = sum(diff(depths) .> 0)

@show p1

p2 = sum(diff([sum(depths[i:i+2]) for i in 1:length(depths)-2]) .> 0)

@show p2


text:
using PMDPs
using Test

@testset "CountingProcesses.jl" begin
    include("bernoulli_process.jl")
end


text:
using Thebes, Luxor, Random

function cuboid(cpos::Point3D, w, h, d)
    pin(cpos + (w/2, h/2, d/2), cpos + (-w/2, h/2, d/2))
end

function frame(scene, framenumber)
    background("midnightblue")
    setlinejoin("bevel")
    setline(5.0)

    eased_n = scene.easingfunction(framenumber, 0, 1, scene.framerange.stop)
    eased_n_2π = eased_n * 2π

    perspective(0)
    eyepoint(Point3D(100 * cos(eased_n_2π), 100 * sin(eased_n_2π), 50 + 20sin(eased_n_2π)))

    carpet()
    axes3D()
    temp = Point3D[]
    k = 19
    for z in -50:k:50
        for y in -50:k:50
            for x in -50:k:50
                push!(temp, Point3D(x, y, z))
            end
        end
    end
    # sort by distance
    sort!(temp, lt = (a, b) -> distance(a, eyepoint()) > distance(b, eyepoint()))
    for pt in temp
        sethue(rescale.([pt.x, pt.y, pt.z], -50, 50, 0, 1)...)
        pin(pt, gfunction = (pt3, pt2) -> begin
            d = rescale(distance(pt3, eyepoint()), 1, 200, 10, 5)
            circle(pt2, d, :fill)
        end)
    end
end

animation = Movie(400, 400, "thebes")

Random.seed!(42)

animate(animation, [Scene(animation, frame, 1:100)],
    framerate=15,
    creategif=true,
    pathname="animation.gif")


text:
@testset "Segment Intersection" begin
    # generate n end-to-end segments
    t = SegmentMap{Int,Int}()
    segments = Any[]
    a = 1
    for i in 1:10000
        b = a + rand(0:100)
        t[(a, b)] = i
        push!(segments, (a, b))
        a = b + 1
    end

    # one
    x = rand(0:a - 1)
    @test length(collect(intersect(t, (x, x)))) === 1
    @test hasintersection(t, x)

    # nothing
    @test length(collect(intersect(t, (a, a)))) === 0
    @test !hasintersection(t, a)

    # everything
    @test length(collect(intersect(t, (segments[1][1], segments[end][2])))) == length(t)

    # always false predicate
    @test length(collect(intersect(t, (segments[1][1], segments[end][2]), (a, b) -> false))) == 0

    # half true
    @test length(collect(intersect(t, (segments[1][1], segments[end][2]), (a, b) -> isodd(a.value)))) == div(length(segments), 2)

    @test all(Bool[hasintersection(t, segment[1]) for segment in segments])
    @test all(Bool[hasintersection(t, segment[end]) for segment in segments])

    # some random intersection queries
    function random_intersection_query()
        i = rand(1:length(segments))
        j = rand(i:length(segments))
        return length(collect(intersect(t, (segments[i][1], segments[j][2])))) == j - i + 1
    end

    @test all(Bool[random_intersection_query() for i in 1:1000])

    # segments separated by 1
    t = SegmentMap{Int,Int}()

    @test !hasintersection(t, 1)

    segments = Any[]
    gaps = Any[]
    a = 1
    for i in 1:10000
        b = a + rand(0:100)
        if iseven(b)
            t[a, b] = i
            push!(segments, (a, b))
        else
            push!(gaps, (a, b))
        end
        a = b + 1
    end

    @test all(Bool[hasintersection(t, segment[1]) for segment in segments])
    @test all(Bool[hasintersection(t, segment[end]) for segment in segments])
    @test !all(Bool[hasintersection(t, segment[1]) for segment in gaps])
    @test !all(Bool[hasintersection(t, segment[end]) for segment in gaps])

    # firstintersection with lower bound
    T = SegmentTrees.SegmentBTree{Int,Segment{Int},4}

    xs = [
        # leaf node 1
        Segment(1, 1),
        Segment(1, 4),
        Segment(2, 2),
        Segment(2, 2),

        # leaf node 2
        Segment(3, 3),
        Segment(3, 3),
        Segment(3, 3),
        Segment(3, 3),

        # leaf node 3
        Segment(4, 4),
        Segment(4, 4),
        Segment(4, 4),
        Segment(4, 4),
    ]

    t = T([Segment(4, 4)])

    @test !(SegmentTrees.firstintersection(t.root, Segment(4, 4), Segment(2, 2)) === nothing)

    x = SegmentTrees.Intersection{Int,Segment{Int},4}()
    SegmentTrees.firstintersection!(t, Segment(4, 4), Segment(2, 2),
                                     x, SegmentTrees.true_cmp)
    @test x.index != 0
end

@testset "Nonunique" begin
    t = SegmentTree{Int,SegmentValue{Int,Int}}()
    n = 1000
    for i in 1:n
        push!(t, SegmentValue{Int,Int}(50, 100, i))
    end
    push!(t, SegmentValue{Int,Int}(1, 200, n + 1))
    push!(t, SegmentValue{Int,Int}(1, 1, n + 2))
    push!(t, SegmentValue{Int,Int}(200, 200, n + 3))

    @test (length(collect(t)) == length(t) == n + 3)
    @test issorted(collect(keys(t)))
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)
end

@testset "Tree Intersection" begin
    n = 10000
    t1 = SegmentTrees.SegmentMap{Int,Int}()
    t2 = SegmentTrees.SegmentMap{Int,Int}()
    global maxend = 1000000
    for k in 1:n
        # generate small-ish segments so we avoid the worst case
        # of O(n^2) intersecting pairs and this test runs quickly
        u = rand(1:maxend)
        v = rand(u:u + 1000)
        t1[(u, v)] = k

        u = rand(1:maxend)
        v = rand(u:u + 1000)
        t2[(u, v)] = k
    end

    # Since SegmentTrees has two tree intersection algorithms, I'm
    # testing by checking that they are in agreement.
    a = intersect(t1, t2, method=:successive)
    b = intersect(t1, t2, method=:iterative)
    @test collect(a) == collect(b)
    @test length(collect(a)) == length(collect(b)) > 0

    # test filter predicate
    count = 0
    for (a, b) in intersect(t1, t2)
        if a.value == b.value
            count += 1
        end
    end
    @test length(collect(intersect(t1, t2, (a, b) -> a.value == b.value))) == count

    ## handle a particular intersection case that may not otherwise get hit
    t1 = SegmentMap{Int,Int}()
    t1[(1, 2)] = 1
    t2 = SegmentMap{Int,Int}()
    t2[(1001, 1002)] = 2
    @test isempty(collect(intersect(t1, t2, method=:iterative)))
    @test isempty(collect(intersect(t1, t2, method=:successive)))
    @test isempty(collect(intersect(t1, t2)))
end


text:
#using Revise, Test
using BifurcationKit, LinearAlgebra, Setfield
const BK = BifurcationKit

function test_newton(x0)
	Ty = eltype(x0)

	F(x, p) = x.^3 .- 1
	Jac(x, p) = diagm(0 => 3 .* x.^2)

	opts = NewtonPar( tol = Ty(1e-8))
	sol, hist, flag, _ = newton(F, Jac, x0, nothing, opts, normN = x->norm(x,Inf))
end
####################################################################################################
# we test the regular newton algorithm
# simple case
test_newton(ones(10) .+ rand(10) * 0.1)

# test types for newton
sol, = test_newton(Float16.(ones(10) .+ rand(10) * 0.1))
@test eltype(sol) == Float16
####################################################################################################
function test_newton_palc(x0, p0)
	Ty = eltype(x0)
	N = length(x0)

	θ = Ty(0.2)
	dotθ = BK.DotTheta(dot)

	F(x, p) = x.^3 .- p
	Jac(x, p) = diagm(0 => 3 .* x.^2)

	z0 = BorderedArray(x0, p0)
	τ0 = BorderedArray(rand(Ty, N), convert(typeof(p0), 0.2))
	zpred = BorderedArray(x0, convert(typeof(p0), 0.3))
	optn = NewtonPar{Ty, DefaultLS, DefaultEig}()
	optc = ContinuationPar{Ty, DefaultLS, DefaultEig}(newtonOptions = optn)
	sol, hist, flag, _ = newtonPALC(F, Jac, p0, (@lens _), z0, τ0, zpred, Ty(0.02), θ, optc, dotθ)
end

sol, = test_newton_palc(ones(10) .+ rand(10) * 0.1, 1.)

#test type
sol, = test_newton_palc(Float32.(ones(10) .+ rand(10) * 0.1), Float32(1.))
@test typeof(sol) == BorderedArray{Vector{Float32}, Float32}
####################################################################################################
# test of  deflated problems
_T = Float32
F4def = (x,p) -> (x-1) * (x-2)
deflationOp = DeflationOperator(_T(2), dot, _T(1), [[_T(1)]])
@test eltype(deflationOp) == _T
@test deflationOp(rand(_T,1)) isa _T


text:
using DelimitedFiles
using IterTools

const data = dropdims(readdlm("day10.txt", Int); dims=2)

function wrangle!(data::Vector{Int})
    pushfirst!(data, 0)
    sort!(data)
    return data
end

wrangle!(data)

function count_diff(data::Vector{Int})
    diff_data = diff(data)
    # +1 because there is a fake final point that is 3 above the max value as
    # stated in the instructions.
    return count(==(1), diff_data) * (count(==(3), diff_data) + 1)
end

println("Part 1: $(count_diff(data))")

function part2(data::Vector{Int})
    total = 1
    for consecutive in groupby(==(1), diff(data))
        if length(consecutive) > 1 && consecutive[1] == 1
            total *= 1 + binomial(length(consecutive), 2)
        end
    end
    return total
end

println("Part 2: $(part2(data))")


text:
using Test
#import Pkg; Pkg.add("LinearAlgebraicRepresentation")

include("./Box.jl")
include("./Frustum.jl")
include("./Geometry.jl")
#include("./GLColorBuffer.jl")
#include("./GLMesh.jl")
#include("./GLPhongShader.jl")
#include("./GLShader.jl")
#include("./GLText.jl")
#include("./GLUtils.jl")
#include("./GLVertexArray.jl")
#include("./GLVertexBuffer.jl")
#include("./Matrix.jl")
#include("./Point.jl")
#include("./Quaternion.jl")
#include("./Viewer.jl")
#include("./ViewerGL.jl")


text:
@testset "$TEST $G" begin

g5 = DG(4)
add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)
z = katz_centrality(g5, 0.4)
@test round.(z, digits=2) == [0.32, 0.44, 0.62, 0.56]

end # testset


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to general functions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Backend inference
# ==============================================================================

# minimal other table interface
struct TestVec{T} <: AbstractArray{T,1}
    data::Array{T,1}
end
Base.IndexStyle(::Type{A}) where {A<:TestVec} = Base.IndexCartesian()
Base.size(A::TestVec) = size(getfield(A, :data))
Base.getindex(A::TestVec, index::Int) = getindex(getfield(A, :data), index)
Base.collect(::Type{T}, itr::TestVec) where {T} = TestVec(collect(T, getfield(itr, :data)))

struct MinimalTable
    data::Matrix
    colnames::TestVec
end

Tables.istable(x::MinimalTable) = true
Tables.columnaccess(::MinimalTable) = true
Tables.columnnames(x::MinimalTable) = getfield(x, :colnames)
Tables.columns(x::MinimalTable) = x
Base.getindex(x::MinimalTable, i1, i2) = getindex(getfield(x, :data), i1, i2)
Base.getproperty(x::MinimalTable, s::Symbol) = getindex(x, :, findfirst(==(s), Tables.columnnames(x)))
Base.convert(::Type{<:TestVec}, x::Array) = TestVec(x)

@testset "Back-end inference" begin
    data = rand(3,3)
    mt = MinimalTable(data, [Symbol("Col. 1"), Symbol("Col. 2"), Symbol("Col. 3")])

    # Text
    # ==========================================================================

    auto   = pretty_table(String, data, tf = unicode)
    manual = pretty_table(String, data, tf = unicode, backend = :text)
    mintable = pretty_table(String, mt, tf = unicode, backend = :text)

    @test auto == manual == mintable

    # HTML
    # ==========================================================================

    auto   = pretty_table(String, data, tf = html_simple)
    manual = pretty_table(String, data, tf = html_simple, backend = :html)
    mintable = pretty_table(String, mt, tf = html_simple, backend = :html)

    @test auto == manual == mintable

    # LaTeX
    # ==========================================================================

    auto   = pretty_table(String, data, tf = latex_default)
    manual = pretty_table(String, data, tf = latex_default, backend = :latex)
    mintable = pretty_table(String, mt, tf = latex_default, backend = :latex)

    @test auto == manual == mintable

    # Error
    # ==========================================================================

    @test_throws TypeError pretty_table(data, tf = [])
end

# Include table in file
# ==============================================================================

@testset "Include Pretty Table to file" begin

    # Text
    # ==========================================================================

    path = "test.txt"

    orig = """
    This is one line.

    This is "another" line.

    % <PrettyTables Table 1> This should be deleted.
    This should be deleted.
    This should be deleted.
    This should be deleted.
    This should be deleted.

    % </PrettyTables>

    <PrettyTables Table 2></PrettyTables>

    <PrettyTables Table 3>This should not be deleted.
    This should not be deleted.
    This should not be deleted.
    This should not be deleted."""

    open(path,"w") do f
        write(f, orig)
    end

    data_table_1 = [1 2 3
                    4 5 6]

    data_table_2 = [7 8 9
                    1 2 3]

    include_pt_in_file(path, "Table 2", data_table_2, tf = mysql,
                       body_hlines = [1])
    include_pt_in_file(path, "Table 1", data_table_1, alignment = :c,
                       show_row_number = true, backup_file = false)
    include_pt_in_file(path, "Table 3", data_table_2, alignment = :c,
                       show_row_number = true, backup_file = false)

    result = read(path, String)
    backup = read(path * "_backup", String)

    expected = """
    This is one line.

    This is "another" line.

    % <PrettyTables Table 1>
    ┌─────┬────────┬────────┬────────┐
    │ Row │ Col. 1 │ Col. 2 │ Col. 3 │
    ├─────┼────────┼────────┼────────┤
    │   1 │   1    │   2    │   3    │
    │   2 │   4    │   5    │   6    │
    └─────┴────────┴────────┴────────┘
    % </PrettyTables>

    <PrettyTables Table 2>
    +--------+--------+--------+
    | Col. 1 | Col. 2 | Col. 3 |
    +--------+--------+--------+
    |      7 |      8 |      9 |
    +--------+--------+--------+
    |      1 |      2 |      3 |
    +--------+--------+--------+
    </PrettyTables>

    <PrettyTables Table 3>This should not be deleted.
    This should not be deleted.
    This should not be deleted.
    This should not be deleted."""

    @test backup == orig
    @test result == expected

    # HTML (`tag_append` option)
    # ==========================================================================

    data_table_1 = [1 2 3
                    4 5 6]

    path = "test.html"

    orig = """

This is a table:

    """

    open(path,"w") do f
        write(f, orig)
    end

    data_table_1 = [1 2 3
                    4 5 6]

    include_pt_in_file(path, "Table 1", data_table_1, backend = :html,
                       standalone = false, tag_append = " -->")

    result = read(path, String)
    backup = read(path * "_backup", String)

    expected = """

This is a table:

    Col. 1
    Col. 2
    Col. 3

    1
    2
    3

    4
    5
    6

    """

    @test backup == orig
    @test result == expected

    # Markdown (`remove_tags` option)
    # ==========================================================================

    data_table_1 = [1 2 3
                    4 5 6]

    path = "test.md"

    orig = """
    # Markdown

    This is a markdown table.

    <PrettyTables Table 1> This should be removed.
    This should be removed.
    This should be removed.
    This should be removed.
    This should be removed.
    </PrettyTables>
    """

    open(path,"w") do f
        write(f, orig)
    end

    data_table_1 = [1 2 3
                    4 5 6]

    include_pt_in_file(path, "Table 1", data_table_1, tf = markdown,
                       remove_tags = true)

    result = read(path, String)
    backup = read(path * "_backup", String)

    expected = """
    # Markdown

    This is a markdown table.

    | Col. 1 | Col. 2 | Col. 3 |
    |--------|--------|--------|
    |      1 |      2 |      3 |
    |      4 |      5 |      6 |

    """

    @test backup == orig
    @test result == expected
end

# Print table to string
# ==============================================================================

@testset "Print table to string" begin
    data     = rand(10,4)
    header_m = [1 2 3 4]
    header_v = [1 2 3 4]

    expected = sprint(pretty_table, data)
    result   = pretty_table(String, data)
    @test result == expected

    expected = sprint((io,data)->pretty_table(io, data, header_m), data)
    result   = pretty_table(String, data, header_m)
    @test result == expected

    expected = sprint((io,data)->pretty_table(io, data, header_v), data)
    result   = pretty_table(String, data, header_v)
    @test result == expected

    dict = Dict(:a => 1,
                :b => 2,
                :c => 3,
                :d => 4)

    expected = sprint(pretty_table, data)
    result   = pretty_table(String, data)
    @test result == expected
end


text:
module TestGenomeGraphs

using GenomeGraphs, BioSequences, Test

end # module

max_issues_repo_name: BioJulia/BioSequenceGraphs

max_forks_repo_name: BioJulia/BioSequenceGraphs


text:
@testset "Thunk" begin
    @test @thunk(3) isa Thunk

    @testset "show" begin
        rep = repr(Thunk(rand))
        @test occursin(r"Thunk\(.*rand.*\)", rep)
    end

    @testset "Externing" begin
        @test extern(@thunk(3)) == 3
        @test extern(@thunk(@thunk(3))) == 3
    end

    @testset "unthunk" begin
        @test unthunk(@thunk(3)) == 3
        @test unthunk(@thunk(@thunk(3))) isa Thunk
    end

    @testset "calling thunks should call inner function" begin
        @test (@thunk(3))() == 3
        @test (@thunk(@thunk(3)))() isa Thunk
    end

    @testset "erroring thunks should include the source in the backtrack" begin
        expected_line = (@__LINE__) + 2  # for testing it is at right palce
        try
            x = @thunk(error())
            extern(x)
        catch err
            err isa ErrorException || rethrow()
            st = stacktrace(catch_backtrace())
            # Should be 2nd last line, as last line will be the `error` function
            stackframe = st[2]
            @test stackframe.line == expected_line
            @test stackframe.file == Symbol(@__FILE__)
        end
    end
end


text:
module ClusterManagers

import Base: launch, manage

export launch, manage
# PBS doesn't have the same semantics as SGE wrt to file accumulate,
# a different solution will have to be found
include("qsub.jl")
include("scyld.jl")
include("condor.jl")

end


text:
using Stars
using Test
cd("test")

# @testset "st_as_sf"
begin
    file_vi = "data/MOD13A2_Henan_2015_2020_10km.tif"
    ga = rast(file_vi)
    mask_shrink = ga.A[:, :, 1] .!= 0

    data, mask, b = st_as_sf(file_vi)

    ## post processing
    r_mask = rast(mask, b)
    lon, lat = st_dim(r_mask)
    d_coord = maskCoords(r_mask)
    @test length(d_coord[:, 1]) == 737
end

# @testset "st_as_stars"
begin
    function test_rastDiff(r1, r2)
        r_diff = (r1 - r2) |> maximum
        delta = filter(!isnan, r_diff.A) |> maximum
        @test delta <= 1e-9
    end

    file_vi = "data/MOD13A2_Henan_2015_2020_10km.tif"
    data, mask, b = st_as_sf(file_vi)
    r_mask = rast(mask, b)

    r1 = st_as_stars(data, mask, b)           # bbox works
    r2 = st_as_stars(data, mask, bbox2vec(b)) # bbox in vector format also works
    r3 = st_as_stars(data, r_mask)

    test_rastDiff(r1, r2)
    test_rastDiff(r2, r3)
end


text:
module ITensorPySCF

using PyCall
using ITensors

export pyscf, fci, gto, scf, ao2mo

const pyscf = PyNULL()
const fci = PyNULL()
const gto = PyNULL()
const scf = PyNULL()
const ao2mo = PyNULL()

function __init__()
  copy!(pyscf, pyimport_conda("pyscf", "pyscf"))
  copy!(fci, pyimport_conda("pyscf.fci", "pycsf"))
  copy!(gto, pyimport_conda("pyscf.gto", "pycsf"))
  copy!(scf, pyimport_conda("pyscf.scf", "pycsf"))
  copy!(ao2mo, pyimport_conda("pyscf.ao2mo", "pycsf"))
end

end


text:
# Let S₊ be the cone of symmetric semidefinite matrices in
# the n*(n+1)/2 dimensional space of symmetric R^{nxn} matrices.
# It is well known that S₊ is a self-dual proper cone.
# Let P₊ be the cone of symmetric semidefinite matrices in
# the n^2 dimensional space of R^{nxn} matrices and
# let D₊ be the cone of matrices A such that A+Aᵀ ∈ P₊.
# P₊ is not proper since it is not solid (as it is not n^2 dimensional) so it is not ensured that (P₊)** = P₊
# However this is the case since, as we will see, (P₊)* = D₊ and (D₊)* = P₊.
# * Let us first see why (P₊)* = D₊.
#   If B is symmetric, then ⟨A,B⟩ = ⟨Aᵀ,Bᵀ⟩ = ⟨Aᵀ,B⟩ so 2⟨A,B⟩ = ⟨A,B⟩ + ⟨Aᵀ,B⟩ = ⟨A+Aᵀ,B⟩
#   Therefore, ⟨A,B⟩ ⩾ 0 for all B ∈ P₊ if and only if ⟨A+Aᵀ,B⟩ ⩾ 0 for all B ∈ P₊
#   Since A+Aᵀ is symmetric and we know that S₊ is self-dual, we have shown that (P₊)*
#   is the set of matrices A such that A+Aᵀ is PSD
# * Let us now see why (D₊)* = P₊.
#   Since A ∈ D₊ implies that Aᵀ ∈ D₊, B ∈ (D₊)* means that ⟨A+Aᵀ,B⟩ ⩾ 0 for any A ∈ D₊ hence B is positive semi-definite.
#   To see why it should be symmetric, simply notice that if B[i,j] < B[j,i] then ⟨A,B⟩ can be made arbitrarily small by setting
#   A[i,j] += s
#   A[j,i] -= s
#   with s arbitrarilly large, and A stays in D₊ as A+Aᵀ does not change.
#
# Typically, SDP primal/dual are presented as
# min ⟨C, X⟩                                                                max ∑ b_ky_k
# ⟨A_k, X⟩ = b_k ∀k                                                         C - ∑ A_ky_k ∈ S₊
#        X ∈ S₊                                                                      y_k free ∀k
# Here, as we allow A_i to be non-symmetric, we should rather use
# min ⟨C, X⟩                                                                max ∑ b_ky_k
# ⟨A_k, X⟩ = b_k ∀k                                                         C - ∑ A_ky_k ∈ P₊
#        X ∈ D₊                                                                      y_k free ∀k
# which is implemented as
# min ⟨C, Z⟩ + (C[i,j]-C[j-i])s[i,j]                                        max ∑ b_ky_k
# ⟨A_k, Z⟩ + (A_k[i,j]-A_k[j,i])s[i,j] = b_k ∀k                   C+Cᵀ - ∑ (A_k+A_kᵀ)y_k ∈ S₊
#       s[i,j] free  1 ⩽ i,j ⩽ n with i > j     C[i,j]-C[j-i] - ∑ (A_k[i,j]-A_k[j,i])y_k = 0  1 ⩽ i,j ⩽ n with i > j
#        Z ∈ S₊                                                                      y_k free ∀k
# where "∈ S₊" only look at the diagonal and upper diagonal part.
# In the last primal program, we have the variables Z = X + Xᵀ and a upper triangular matrix S such that X = Z + S - Sᵀ

"""
    SquareBridge{T, F<:MOI.AbstractVectorFunction,
                 G<:MOI.AbstractScalarFunction,
                 TT<:MOI.AbstractSymmetricMatrixSetTriangle,
                 ST<:MOI.AbstractSymmetricMatrixSetSquare} <: AbstractBridge

The `SquareBridge` reformulates the constraint of a square matrix to be in `ST`
to a list of equality constraints for pair or off-diagonal entries with
different expressions and a `TT` constraint the upper triangular part of the
matrix.

For instance, the constraint for the matrix
```math
\\begin{pmatrix}
  1      & 1 + x & 2 - 3x\\\\
  1 +  x & 2 + x & 3 -  x\\\\
  2 - 3x & 2 + x &     2x
\\end{pmatrix}
```
to be PSD can be broken down to the constraint of the symmetric matrix
```math
\\begin{pmatrix}
  1      & 1 + x & 2 - 3x\\\\
  \\cdot & 2 + x & 3 -  x\\\\
  \\cdot & \\cdot &    2x
\\end{pmatrix}
```
and the equality constraint between the off-diagonal entries (2, 3) and (3, 2)
``2x == 1``. Note that now symmetrization constraint need to be added between
the off-diagonal entries (1, 2) and (2, 1) or between (1, 3) and (3, 1) since
the expressions are the same.
"""
struct SquareBridge{
    T,
    F<:MOI.AbstractVectorFunction,
    G<:MOI.AbstractScalarFunction,
    TT<:MOI.AbstractSymmetricMatrixSetTriangle,
    ST<:MOI.AbstractSymmetricMatrixSetSquare,
} <: AbstractBridge
    square_set::ST
    triangle::CI{F,TT}
    sym::Vector{Pair{Tuple{Int,Int},CI{G,MOI.EqualTo{T}}}}
end
function bridge_constraint(
    ::Type{SquareBridge{T,F,G,TT,ST}},
    model::MOI.ModelLike,
    f::F,
    s::ST,
) where {T,F,G,TT,ST}
    f_scalars = MOIU.eachscalar(f)
    sym = Pair{Tuple{Int,Int},CI{G,MOI.EqualTo{T}}}[]
    dim = MOI.side_dimension(s)
    upper_triangle_indices = Int[]
    trilen = div(dim * (dim + 1), 2)
    sizehint!(upper_triangle_indices, trilen)
    k = 0
    for j in 1:dim
        for i in 1:j
            k += 1
            push!(upper_triangle_indices, k)
            # We constrain the entries (i, j) and (j, i) to be equal
            upper = f_scalars[i+(j-1)*dim]
            lower = f_scalars[j+(i-1)*dim]
            diff = MOIU.operate!(-, T, upper, lower)
            MOIU.canonicalize!(diff)
            # The value 1e-10 was decided in https://github.com/jump-dev/JuMP.jl/pull/976
            # This avoid generating symmetrization constraints when the
            # functions at entries (i, j) and (j, i) are almost identical
            if !MOIU.isapprox_zero(diff, 1e-10)
                if MOIU.isapprox_zero(diff, 1e-8)
                    @warn "The entries ($i, $j) and ($j, $i) of the" *
                          " positive semidefinite constraint are almost" *
                          " identical but a constraint is added to ensure their" *
                          " equality because the largest difference between the" *
                          " coefficients is smaller than 1e-8 but larger than" *
                          " 1e-10."
                end
                push!(
                    sym,
                    (i, j) => MOIU.normalize_and_add_constraint(
                        model,
                        diff,
                        MOI.EqualTo(zero(T)),
                        allow_modify_function = true,
                    ),
                )
            end
        end
        k += dim - j
    end
    @assert length(upper_triangle_indices) == trilen
    triangle = MOI.add_constraint(
        model,
        f_scalars[upper_triangle_indices],
        MOI.triangular_form(s),
    )
    return SquareBridge{T,F,G,TT,ST}(s, triangle, sym)
end

function MOI.supports_constraint(
    ::Type{SquareBridge{T}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{<:MOI.AbstractSymmetricMatrixSetSquare},
) where {T}
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:SquareBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(
    ::Type{SquareBridge{T,F,G,TT,ST}},
) where {T,F,G,TT,ST}
    return [(F, TT), (G, MOI.EqualTo{T})]
end
function concrete_bridge_type(
    ::Type{<:SquareBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    ST::Type{<:MOI.AbstractSymmetricMatrixSetSquare},
) where {T}
    S = MOIU.scalar_type(F)
    G = MOIU.promote_operation(-, T, S, S)
    TT = MOI.triangular_form(ST)
    return SquareBridge{T,F,G,TT,ST}
end

# Attributes, Bridge acting as a model
function MOI.get(
    ::SquareBridge{T,F,G,TT},
    ::MOI.NumberOfConstraints{F,TT},
) where {T,F,G,TT}
    return 1
end
function MOI.get(
    bridge::SquareBridge{T,F,G},
    ::MOI.NumberOfConstraints{G,MOI.EqualTo{T}},
) where {T,F,G}
    return length(bridge.sym)
end
function MOI.get(
    bridge::SquareBridge{T,F,G,TT},
    ::MOI.ListOfConstraintIndices{F,TT},
) where {T,F,G,TT}
    return [bridge.triangle]
end
function MOI.get(
    bridge::SquareBridge{T,F,G},
    ::MOI.ListOfConstraintIndices{G,MOI.EqualTo{T}},
) where {T,F,G}
    return map(pair -> pair.second, bridge.sym)
end

# Indices
function MOI.delete(model::MOI.ModelLike, bridge::SquareBridge)
    MOI.delete(model, bridge.triangle)
    for pair in bridge.sym
        MOI.delete(model, pair.second)
    end
end

# Attributes, Bridge acting as a constraint
function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::SquareBridge{T},
) where {T}
    tri = MOIU.eachscalar(MOI.get(model, attr, bridge.triangle))
    dim = MOI.side_dimension(bridge.square_set)
    sqr = Vector{eltype(tri)}(undef, dim^2)
    sqrmap(i, j) = (j - 1) * dim + i
    k = 0
    for j in 1:dim
        for i in 1:j
            k += 1
            sqr[sqrmap(i, j)] = tri[k]
            sqr[sqrmap(j, i)] = tri[k]
        end
    end
    for sym in bridge.sym
        i, j = sym.first
        diff = MOI.get(model, attr, sym.second)
        set = MOI.get(model, MOI.ConstraintSet(), sym.second)
        upper = sqr[sqrmap(i, j)]
        lower = MOIU.operate(-, T, upper, diff)
        lower = MOIU.operate!(-, T, lower, MOI.constant(set))
        sqr[sqrmap(j, i)] = MOIU.convert_approx(eltype(tri), lower)
    end
    return MOIU.vectorize(sqr)
end
function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::SquareBridge)
    return bridge.square_set
end
function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimal,
    bridge::SquareBridge{T},
) where {T}
    tri = MOI.get(model, attr, bridge.triangle)
    dim = MOI.side_dimension(bridge.square_set)
    sqr = Vector{eltype(tri)}(undef, dim^2)
    k = 0
    for j in 1:dim
        for i in 1:j
            k += 1
            sqr[i+(j-1)*dim] = sqr[j+(i-1)*dim] = tri[k]
        end
    end
    return sqr
end
function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDual,
    bridge::SquareBridge,
)
    tri = MOI.get(model, attr, bridge.triangle)
    dim = MOI.side_dimension(bridge.square_set)
    sqr = Vector{eltype(tri)}(undef, dim^2)
    k = 0
    for j in 1:dim
        for i in 1:j
            k += 1
            # The triangle constraint uses only the upper triangular part
            if i == j
                sqr[i+(j-1)*dim] = tri[k]
            else
                sqr[i+(j-1)*dim] = 2tri[k]
                sqr[j+(i-1)*dim] = zero(eltype(sqr))
            end
        end
    end
    for pair in bridge.sym
        i, j = pair.first
        dual = MOI.get(model, attr, pair.second)
        sqr[i+(j-1)*dim] += dual
        sqr[j+(i-1)*dim] -= dual
    end
    return sqr
end


text:
function createPointsOcta(n::Int64)
    points,triangulation  = computexyz(n)
	return points,triangulation
end

function linspace(start,stop,n)
	return collect(range(start,stop=stop,length=n))
end

function computexyz(n::Int64,slerpflag=true::Bool)
	pt0 = [0.0, 0.0, 1.0]
	pt1 = [0.0, 1.0, 0.0]
	pt2 = [1.0, 0.0, 0.0]

	if slerpflag
		# use slerp instead of linspace
		pts01 = slerp(pt0,pt1,n)
		pts02 = slerp(pt0,pt2,n)
	else
		pts01 = linspace(pt0,pt1,n) #type is: n-element LinSpace{Array{Float64,1}}:
		pts02 = linspace(pt0,pt2,n)
	end

	nptsoctant = Int64(n*(n+1)/2)
	pts = zeros(3,nptsoctant)

	# Creates points in planar geometry
	counter = 0;
	for i=1:n
		if slerpflag
			if i==1
				tmp = pts01[:,1]
			else
				tmp = slerp(pts01[:,i],pts02[:,i],i)
			end
		else
			tmp =linspace(pts01[i],pts02[i],i)
		end
		for j=1:i
			counter += 1
			if slerpflag
				pts[:,counter] = tmp[:,j]
			else
				pts[:,counter] = tmp[j]
			end
		end

	end

	# Project points onto sphere
	for i=1:nptsoctant
		pts[:,i] = pts[:,i]/norm(pts[:,i])
	end

	# We now want to get the connectivity. Therefore enumerate all
	#the points (overkill from a computational point of view) and then
	#write the correct connectivity.

    # Matrix that assigns an ID to the points
   ids=zeros(Int64,n,n); # Too large, but not important.
   #Matrix that will later contain all triangles
   nTrianglesOctant = Int64(n*n-2*n+1);
   triangles = zeros(Int64,3,nTrianglesOctant);

	counter = 0
	for i=1:n
		for j=1:i
			counter +=1
			ids[i,j] = counter;
		end
	end

	# Now create triangles
	counter = 0
	tmp = zeros(Int64,1)
	for i=1:n
		for j=1:i-1
			tmp  = [ids[i,j],ids[i,j+1],ids[i-1,j]]
			counter += 1
			triangles[:,counter] = tmp;
		end
		if i<n
			for j=1:i-1
				tmp = [ids[i,j],ids[i,j+1],ids[i+1,j+1]]
				counter +=1
				triangles[:,counter] = tmp;
			end
		end
	end

	# Now we have the quadrature points and triangles for a single octant
	ptsAll = deepcopy(pts);
	#mat ptsAll = pts; // the quadrature points for all octants, but with duplicated #entries

	tmp = deepcopy(pts); tmp[1,:] *= -1.0
	ptsAll = deepcopy(hcat(ptsAll,tmp));

	tmp = deepcopy(pts); tmp[2,:] *= -1.0;
	ptsAll = deepcopy(hcat(ptsAll,tmp));

	tmp = deepcopy(pts); tmp[1,:] *= -1.0; tmp[2,:] *= -1.0;
	ptsAll = deepcopy(hcat(ptsAll,tmp));

	tmp = deepcopy(ptsAll); tmp[3,:] *= -1.0;
	ptsAll = deepcopy(hcat(ptsAll,tmp));

	trianglesAll = deepcopy(triangles);
	for i=1:7
		trianglesAll = (hcat(trianglesAll,triangles .+i .*nptsoctant))
	end
	ptsAll,triangulation =unique(ptsAll,trianglesAll)

	ptsAll = permutedims(ptsAll)
	triangulation = permutedims(triangulation)
	return ptsAll, triangulation
end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
using JLD, Plots, LaTeXStrings
font = Plots.font("DejaVu Sans", 24)
pyplot(size=(1200,800),linewidth=4,xscale=:log10,yscale=:log10,
       xtickfont=font,ytickfont=font,markersize=8,
       legendfont=font,guidefont=font,titlefont=font,
       xlabel="Δx",ylabel="||error||",
       title="Stokes Equation (constant μ) Convergence",grid=false)

function main()
    # import jld file
    jldfile = "data/TEMP_stokes2D_const_validation.jld"

    VL1arr   = load(jldfile,"VL1arr")
    VL2arr   = load(jldfile,"VL2arr")
    VLInfarr = load(jldfile,"VLInfarr")
    PL1arr   = load(jldfile,"PL1arr")
    PL2arr   = load(jldfile,"PL2arr")
    PLInfarr = load(jldfile,"PLInfarr")

    harr     = load(jldfile,"harr")
    κarr     = load(jldfile,"κarr")

    timearr  = load(jldfile,"timearr")

    # create plots
    # error plots
    p1 = plot(harr,PL1arr,color=:red,linestyle=:dash,label="L1 (pressure)")
    plot!(p1,harr,PL2arr,color=:blue,linestyle=:dash,label="L2")
    plot!(p1,harr,PLInfarr,color=:green,linestyle=:dash,label="Linf")
    scatter!(p1,harr,PL1arr,color=:red,label="")
    scatter!(p1,harr,PL2arr,color=:blue,label="")
    scatter!(p1,harr,PLInfarr,color=:green,label="")
    plot!(p1,harr,harr.^2,color=:black,linestyle=:dash,label="O(Δx^2)")

    plot!(p1,harr,VL1arr,color=:red,label="L1 (velocity)")
    plot!(p1,harr,VL2arr,color=:blue,label="L2")
    plot!(p1,harr,VLInfarr,color=:green,label="Linf")
    scatter!(p1,harr,VL1arr,color=:red,label="")
    scatter!(p1,harr,VL2arr,color=:blue,label="")
    scatter!(p1,harr,VLInfarr,color=:green,label="")
    plot!(p1,harr,0.01*harr.^4,color=:black,label="O(Δx^4)")

    # time plots
    p2 = plot(title="Stokes (constant μ) Time-Accuracy",xlabel="||error||",ylabel="elapsed time (s)")
    colorArr = [:red,:green,:blue,:orange,:yellow,:purple]
    for i=1:length(timearr)
        scatter!([PLInfarr[i]],[timearr[i]],
                 color=colorArr[i],label="",markershape=:star,
                 markersize=16)
        scatter!([VLInfarr[i]],[timearr[i]],
                 color=colorArr[i],label="",markershape=:circle,
                 markersize=12)
    end

    # condition number plot
    p3 = plot(harr,κarr,label="",linewidth=8,color=:blue)
    scatter!(p3,harr,κarr,label="",ylabel=L"\kappa",
            title="Stokes 2D (constant μ) condition number",markersize=10,
            color=:blue)

    # export plots
    run(`mkdir -p figures`)
    plot(p1)
    savefig(p1,"figures/TEMPstokes2D_const_validation.png")
    plot(p2)
    savefig(p2,"figures/TEMPstokes2D_const_time_accuracy.png")
    plot(p3)
    savefig(p3,"figures/TEMPstokes2D_const_cond.png")

    nothing
end

main()

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
using SurrealNumbers
out_dir = "Data/"

# addition table
x = convert.(SurrealFinite, [0, 1/2, 1, 2])
A = zeros(SurrealFinite, length(x), length(x))
B = zeros(Int, length(x), length(x))
for i=1:length(x)
    for j=1:length(x)
        A[i,j] = x[i] + x[j]
        B[i,j] = generation(x[i]) + generation(x[j])
    end
end
float.(A)
A
generation.(A)
# A .== canonicalise.(A)
iscanonical.(A)
B

a10 = convert(SurrealFinite, 3/2)
a11 = convert(SurrealFinite, 1/2) + convert(SurrealFinite, 1)
a12 = convert(SurrealFinite, 3/4) + convert(SurrealFinite, 3/4)

file = "$(out_dir)addition_ex_10.dot"
FID = open(file, "w")
surreal2dag(FID, a10)
close(FID)
run(`dot -Tpdf -O $file`)

file = "$(out_dir)addition_ex_11.dot"
FID = open(file, "w")
surreal2dag(FID, a11)
close(FID)
run(`dot -Tpdf -O $file`)

file = "$(out_dir)addition_ex_12.dot"
FID = open(file, "w")
surreal2dag(FID, a12)
close(FID)
run(`dot -Tpdf -O $file`)


text:
using Gadfly
using DataFrames

include("../utility.jl")

# ------------------------------------------------------------------------------
# coût de l'allocation aux petits n
# ------------------------------------------------------------------------------

function main1()
    N = [1,2,3,4,5,6,7,8,9,
         10,20,30,40,50,60,70,80,90,100,200,300,400,500,750,
         1_000,2_500,5_000,7_500,10_000,100_000,1_000_000]

    cpu = Vector{Float64}(length(N))
    @inbounds for i in 1:length(N)
        n = N[i]
        a = rand(n)
        # sqrt(a)
        gc()
        tsum = 0.0
        tsum2 = 0.0

        ncycle = 10_000 ÷ n + 1
        nrep = 500

        for j in 1:nrep
            t = @elapsed for k in 1:ncycle
                Vector{Float64}(n)
                # sqrt(a)
            end
            t = t / ncycle / n * 1e9
            tsum += t
            # tsum2 += t^2
        end
        #
        tavg = tsum / nrep
        # tavg = (@timed for j in 1:nrep Vector{Float64}(n) end)[2] / nrep / n * 1e9
        # tavg2 = tsum2 / nrep
        # tunc = sqrt((tavg2 - tavg^2) / (nrep - 1))

        # tavg = mean([@elapsed sqrt(a) for i in 1:100]) * (1e9 / n)
        cpu[i] = tavg
    end
    N, cpu
end

function main2()
    N = [1,2,3,4,5,6,7,8,9,
         10,20,30,40,50,60,70,80,90,100,200,300,400,500,750,
         1_000,2_500,5_000,7_500,10_000,100_000,1_000_000]

    cpu = Vector{Float64}(length(N))
    @inbounds for i in 1:length(N)
        n = N[i]
        a = rand(n)
        sqrt(a)
        gc()
        tsum = 0.0
        tsum2 = 0.0

        ncycle = 10_000 ÷ n + 1
        nrep = 100

        for j in 1:nrep
            t = @elapsed for k in 1:ncycle
                # Vector{Float64}(n)
                sqrt(a)
            end
            t = t / ncycle / n * 1e9
            tsum += t
            # tsum2 += t^2
        end
        #
        tavg = tsum / nrep
        # tavg = (@timed for j in 1:nrep sqrt(a) end)[2] / nrep / n * 1e9
        # tavg2 = tsum2 / nrep
        # tunc = sqrt((tavg2 - tavg^2) / (nrep - 1))

        # tavg = mean([@elapsed sqrt(a) for i in 1:100]) * (1e9 / n)
        cpu[i] = tavg
    end
    N, cpu
end

layers = Layer[]

main1()
main2()

df = DataFrame() ; df[:N], df[:CPU1] = main1()
push!(layers,layer(df, x="N", y="CPU1", Geom.line, Theme(default_color=color("blue")))[1])

df[:CPU2] = main2()[2]
push!(layers,layer(df, x="N", y="CPU2", Geom.line, Theme(default_color=color("yellow")))[1])
println(df)
p = Gadfly.plot(layers,
                # Scale.y_log10,
                Scale.x_log10,
                Guide.xlabel("n-element vector"),
                Guide.ylabel("CPU time in nsec/element"),
                Guide.title("CPU time with n elements"),
                Guide.manual_color_key("",
                ["ALLOC", "Base.sqrt"],
                ["blue","yellow"])
                )

path = Pkg.dir("MKL") * "/benchmark/cpu_regarding_n/"
draw(PNG(path*"sqrt_yichao2.png", 20cm, 20cm), p)

# ------------------------------------------------------------------------------
# comparaison des différentes méthodes de bench : problème pour les petits n
# ------------------------------------------------------------------------------

function main3()
    N = [1,2,3,4,5,6,7,8,9,
         10,20,30,40,50,60,70,80,90,100,200,300,400,500,750,
         1_000,2_500,5_000,7_500,10_000,100_000,1_000_000]

    df = DataFrame(N=Int64[],ALLOC=Float64[],CPU1=Float64[],CPU2=Float64[],CPU3=Float64[],CPU4=Float64[],CPU5=Float64[])

    @inbounds for i in 1:length(N)
        n = N[i]
        a = rand(n)
        sqrt(a)

        gc()

        tavg1 = mean([@elapsed sqrt(a) for i in 1:100]) * (1e9 / n)

        ncycle = 10_000 ÷ n + 1
        nrep = 1000
        tsum = 0.0
        for j in 1:nrep
            t = @elapsed for k in 1:ncycle sqrt(a) end
            tsum+=  t / ncycle / n * 1e9
        end
        tavg2 = tsum / nrep

        tavg3 = (@timed for j in 1:nrep sqrt(a) end)[2]
        tavg3 = tavg3 / nrep / n * 1e9

        tavg4 = ibench2(sqrt,(a,),n,0.1)[1][2]*1e9/n

        tavg5 = ibench(sqrt,(a,),0.1)[1][2]*1e9/n

        talloc = 0.0

        # talloc = (@timed for j in 1:nrep Array(Float64,n) end)[2] / nrep / n * 1e9

        push!(df,[n,talloc, tavg1,tavg2,tavg3,tavg4,tavg5])
    end
    df
end

df = main3();
println(df)

# ------------------------------------------------------------------------------
# comparaison allocation vs fonction sur sqrt
# ------------------------------------------------------------------------------

function main4()
    N = [1,2,3,4,5,6,7,8,9,
         10,20,30,40,50,60,70,80,90,100,200,300,400,500,750,
         1_000,2_500,5_000,7_500,10_000,100_000,1_000_000]

    df = DataFrame(N=Int64[],ALLOC=Float64[],SQRT=Float64[])
    @inbounds for i in 1:length(N)
        T = Float64
        n = N[i]
        a = rand(n)
        sqrt(a)
        nrep = 1000

        gc()
        # tcpu = (@timed for j in 1:nrep sqrt(a) end)[2]
        # tcpu = tcpu / nrep / n * 1e9
        ncycle = 10_000 ÷ n + 1
        tcpu = 0.0
        for j in 1:nrep
            tcpu += @elapsed for k in 1:ncycle sqrt(a) end
        end
        tcpu = tcpu / nrep / ncycle / n * 1e9

        gc()
        # talloc = (@timed for j in 1:nrep Vector{T}(n) end)[2]
        # talloc = talloc / nrep / n * 1e9
        ncycle = 10_000 ÷ n + 1
        talloc = 0.0
        for j in 1:nrep
            talloc += @elapsed for k in 1:ncycle Vector{T}(n) end
        end
        talloc = talloc / nrep / ncycle / n * 1e9

        push!(df,[n,talloc, tcpu])
    end
    df
end

df = main4();
df[:diff] = df[:SQRT]-df[:ALLOC]
println(df)