text:
@testset "Neighbors Full" begin

    @test hasmethod(AgentBasedModels.argumentsFull_!,(AgentBasedModels.Program_,String))
    @test hasmethod(AgentBasedModels.loopFull_,(AgentBasedModels.Program_,Expr,String))

    @test AgentBasedModels.argumentsFull_!(AgentBasedModels.Program_(@agent(3)),"gpu") == Nothing

    @test_nowarn begin

        m = @agent(
            0,

            [l1,l2]::Local,
            [i1,i2]::Identity,

            UpdateLocalInteraction = begin
                if sqrt((l1.j-l1.i)^2+(l2.j-l2.i)^2) < 5
                    l1.i += 1
                end
            end
        )
        m = compile(m)
    end

    m = @agent(
        0,

        [l1,l2,intLocal,int]::Local,
        [i1,i2]::Identity,

        Equation = begin
            nothing
        end,

        UpdateInteraction = begin
            if sqrt((l1.j-l1.i)^2+(l2.j-l2.i)^2) < 5
                int.i += 1
            end
        end,

        UpdateLocalInteraction = begin
            if sqrt((l1.j-l1.i)^2+(l2.j-l2.i)^2) < 5
                intLocal.i += 1
            end
        end
    )
    m = compile(m,debug=false)

    @test begin
        com = Community(m,N=4)
        com.l1 = [1,1,0,0]
        com.l2 = [0,1,1,0]

        comt = m.evolve(com,dt=0.1,tMax=1)
        comt[1].int == [4,4,4,4]
    end

    @test begin
        com = Community(m,N=4)
        com.l1 = [1,1,0,0]
        com.l2 = [0,1,1,0]

        comt = m.evolve(com,dt=0.1,tMax=1)
        comt[1].intLocal == [4,4,4,4]
    end

    @test begin
        com = Community(m,N=4)
        com.l1 = [7,7,0,0]
        com.l2 = [0,7,7,0]

        comt = m.evolve(com,dt=0.1,tMax=1)
        comt[1].int == [1,1,1,1]
    end

    @test begin
        com = Community(m,N=4)
        com.l1 = [7,7,0,0]
        com.l2 = [0,7,7,0]

        comt = m.evolve(com,dt=0.1,tMax=1)
        comt[1].intLocal == [1,1,1,1]
    end

end

@testset "Neighbors Grid" begin

    @test hasmethod(AgentBasedModels.argumentsGrid_!,(AgentBasedModels.Program_, String))
    @test hasmethod(AgentBasedModels.loopGrid_,(AgentBasedModels.Program_, Expr, String))

    @test begin
        xx = []
        for x in -2.:1.:2.
            push!(xx,AgentBasedModels.position2gridVectorPosition_(x,-1.5,1.,5))
        end
        xx == Array(1:5)
    end

    @test begin
        xx = []
        for y in -3.:1.:3.
            for x in -2.:1.:2.
                push!(xx,AgentBasedModels.position2gridVectorPosition_(x,-1.5,1.,5,y,-2.5,1.,7))
            end
        end
        xx == Array(1:5*7)
    end

    @test begin
        xx = []
        for z in -4.:1.:4.
            for y in -3.:1.:3.
                for x in -2.:1.:2.
                    push!(xx,AgentBasedModels.position2gridVectorPosition_(x,-1.5,1.,5,y,-2.5,1.,7,z,-3.5,1.,9))
                end
            end
        end
        xx == Array(1:5*7*9)
    end

    @test begin
        xx = []
        for x in 1:3
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(1,x,5,false))
        end
        xx == [-1,1,2]
    end

    @test begin
        xx = []
        for x in 1:3
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(2,x,5,false))
        end
        xx == [1,2,3]
    end

    @test begin
        xx = []
        for x in 1:3
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(2,x,5,true))
        end
        xx == [4,2,3]
    end

    @test begin
        xx = []
        for x in 1:3
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(4,x,5,true))
        end
        xx == [3,4,2]
    end

    @test begin
        xx = []
        for x in 1:9
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(1,x,5,false,7,false))
        end
        xx == [-1,-1,-1,-1,1,2,-1,6,7]
    end

    @test begin
        xx = []
        for x in 1:9
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(7,x,5,false,7,false))
        end
        xx == [1,2,3,6,7,8,11,12,13]
    end

    @test begin
        xx = []
        for x in 1:9
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(25,x,5,false,5,false))
        end
        xx == [19,20,-1,24,25,-1,-1,-1,-1]
    end

    @test begin
        xx = []
        for x in 1:9
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(2,x,5,true,7,false))
        end
        xx == [-1,-1,-1,4,2,3,9,7,8]
    end

    @test begin
        xx = []
        for x in 1:9
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(4,x,5,true,7,false))
        end
        xx == [-1,-1,-1,3,4,2,8,9,7]
    end

    @test begin
        xx = []
        for x in 1:9
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(6,x,5,false,7,true))
        end
        xx == [-1,26,27,-1,6,7,-1,11,12]
    end

    @test begin
        xx = []
        for x in 1:9
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(19,x,5,true,5,true))
        end
        xx == [13,14,12,18,19,17,8,9,7]
    end

    @test begin
        xx = []
        for x in 1:27
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(1,x,5,false,5,false,5,false))
        end
        xx == [-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,1,2,-1,6,7,-1,-1,-1,-1,26,27,-1,31,32]
    end

    @test begin
        xx = []
        for x in 1:27
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(125,x,5,false,5,false,5,false))
        end
        xx == [94,95,-1,99,100,-1,-1,-1,-1,119,120,-1,124,125,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1]
    end

    @test begin
        xx = []
        for x in 1:27
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(32,x,5,false,5,false,5,false))
        end
        xx == [1,2,3,6,7,8,11,12,13,26,27,28,31,32,33,36,37,38,51,52,53,56,57,58,61,62,63]
    end

    @test begin
        xx = []
        for x in 1:27
            push!(xx,AgentBasedModels.gridVectorPositionNeighbour_(32,x,5,true,5,true,5,true))
        end
        xx == [94,92,93,84,82,83,89,87,88,44,42,43,34,32,33,39,37,38,69,67,68,59,57,58,64,62,63]
    end

    m = @agent 3 [w]::Local;

    @test_nowarn AgentBasedModels.argumentsGrid_!(AgentBasedModels.Program_(m),"cpu")

    mo = compile(m,neighbors="grid",debug=false)
    @test_throws ErrorException begin #Error if box not specified
        com = Community(mo,N=12)
        com.x .= [-11,-7,-2,2,7,11,-11,-7,-2,2,7,11]
        com.y .= 0

        comt = mo.evolve(com,dt=0.1,tMax=1)
    end

    for platform in testplatforms

        #Check 1D
        m = @agent(
            1,

            [intLocal,int]::Local,

            Equation = begin
                d(l2) = 0.
            end,

            UpdateInteraction = begin
                if abs(x.i - x.j) <= .9
                    int.i += 1
                end
            end,

            UpdateLocalInteraction = begin
                intLocal.i += 1
            end,
        )

        mo = compile(m,platform=platform,neighbors="grid",debug=false)
        # println(mo.program)
        @test begin
            com = Community(mo,N=4)
            com.x .= [-1.5,-.5,.5,1.5]
            com.simulationBox = [-1. 1.]
            com.radiusInteraction = 1.1

            comt = mo.evolve(com,dt=0.1,tMax=1)
            (comt[1].intLocal == [2,3,3,2]) && (comt[1].int == [1,1,1,1])
        end

        #Check 2D
        m = @agent(
            2,

            [intLocal,int]::Local,

            Equation = begin
                d(l2) = 0.
            end,

            UpdateInteraction = begin
                if sqrt((x.i - x.j)^2+(y.i - y.j)^2) <= 1.1
                    int.i += 1
                end
            end,

            UpdateLocalInteraction = begin
                intLocal.i += 1
            end,
        )

        mo = compile(m,platform=platform,neighbors="grid",debug=false)
        # println(mo.program)
        @test begin
            com = Community(mo,N=16)
            com.x .= [-1.5,-.5,.5,1.5,-1.5,-.5,.5,1.5,-1.5,-.5,.5,1.5,-1.5,-.5,.5,1.5]
            com.y .= [-1.5,-1.5,-1.5,-1.5,-.5,-.5,-.5,-.5,.5,.5,.5,.5,1.5,1.5,1.5,1.5]
            com.simulationBox = [-1. 1.;-1. 1.]
            com.radiusInteraction = 1.1

            comt = mo.evolve(com,dt=0.1,tMax=1)
            (comt[1].intLocal == [4.,6.,6.,4.,6.,9.,9.,6.,6.,9.,9.,6.,4.,6.,6.,4.]) && (comt[1].int == [3.,4.,4.,3.,4.,5.,5.,4.,4.,5.,5.,4.,3.,4.,4.,3.])
        end

        #Check 3D
        m = @agent(
            3,

            [intLocal,int]::Local,

            Equation = begin
                d(l2) = 0.
            end,

            UpdateInteraction = begin
                if sqrt((x.i - x.j)^2+(y.i - y.j)^2+(z.i - z.j)^2) <= 1.1
                    int.i += 1
                end
            end,

            UpdateLocalInteraction = begin
                intLocal.i += 1
            end,
        )

        mo = compile(m,platform=platform,neighbors="grid",debug=false)
        # println(mo.program)
        @test begin
            com = Community(mo,N=64)
            intLocal = zeros(64)
            intt = zeros(64)
            for (i,x) in enumerate([-1.5,-.5,.5,1.5])
                for (j,y) in enumerate([-1.5,-.5,.5,1.5])
                    for (k,z) in enumerate([-1.5,-.5,.5,1.5])
                        com.x[i+4*(j-1)+16*(k-1)] = x
                        com.y[i+4*(j-1)+16*(k-1)] = y
                        com.z[i+4*(j-1)+16*(k-1)] = z

                        intLocal[i+4*(j-1)+16*(k-1)] = 27.
                        intt[i+4*(j-1)+16*(k-1)] = 7.
                        if abs(x) == 1.5 || abs(y) == 1.5 || abs(z) == 1.5
                            intLocal[i+4*(j-1)+16*(k-1)] = 18.
                            intt[i+4*(j-1)+16*(k-1)] = 6.
                        end
                        if (abs(x) == 1.5 && abs(y) == 1.5) || (abs(x) == 1.5 && abs(z) == 1.5) || (abs(y) == 1.5 && abs(z) == 1.5)
                            intLocal[i+4*(j-1)+16*(k-1)] = 12.
                            intt[i+4*(j-1)+16*(k-1)] = 5.
                        end
                        if abs(x) == 1.5 && abs(y) == 1.5 && abs(z) == 1.5
                            intLocal[i+4*(j-1)+16*(k-1)] = 8.
                            intt[i+4*(j-1)+16*(k-1)] = 4.
                        end
                    end
                end
            end
            com.simulationBox = [-1. 1.;-1. 1.;-1. 1]
            com.radiusInteraction = 1.1
            comt = mo.evolve(com,dt=0.1,tMax=1)
            (comt[1].intLocal == intLocal) && (comt[1].int == intt)
        end

        #Check 1D Periodic
        m = @agent(
            1,

            [intLocal,int]::Local,

            Equation = begin
                d(l2) = 0.
            end,

            UpdateInteraction = begin
                if abs(x.i - x.j) <= .9
                    int.i += 1
                end
            end,

            UpdateLocalInteraction = begin
                intLocal.i += 1
            end,

            Boundary = BoundaryFlat(1,Periodic())
        )

        mo = compile(m,platform=platform,neighbors="grid",debug=false)
        # println(mo.program)
        @test begin
            com = Community(mo,N=4)
            com.x .= [-1.5,-.5,.5,1.5]
            com.simulationBox = [-1.6 1.6]
            com.radiusInteraction = 1.

            comt = mo.evolve(com,dt=0.1,tMax=1)
            (comt[1].intLocal == [3,3,3,3]) && (comt[1].int == [1,1,1,1])
        end

        #Check 2D Periodic
        m = @agent(
            2,

            [intLocal,int]::Local,

            Equation = begin
                d(l2) = 0.
            end,

            UpdateInteraction = begin
                if sqrt((x.i - x.j)^2+(y.i - y.j)^2) <= 1.1
                    int.i += 1
                end
            end,

            UpdateLocalInteraction = begin
                intLocal.i += 1
            end,

            Boundary = BoundaryFlat(2,Periodic())
        )

        mo = compile(m,platform=platform,neighbors="grid",debug=false)
        # println(mo.program)
        @test begin
            com = Community(mo,N=16)
            com.x .= [-1.5,-.5,.5,1.5,-1.5,-.5,.5,1.5,-1.5,-.5,.5,1.5,-1.5,-.5,.5,1.5]
            com.y .= [-1.5,-1.5,-1.5,-1.5,-.5,-.5,-.5,-.5,.5,.5,.5,.5,1.5,1.5,1.5,1.5]
            com.simulationBox = [-2. 2.;-2. 2.]
            com.radiusInteraction = 1.1

            comt = mo.evolve(com,dt=0.1,tMax=1)
            (comt[1].intLocal == [6.,6.,6.,6.,9.,9.,9.,9.,9.,9.,9.,9.,6.,6.,6.,6.]) && (comt[1].int == [3.,4.,4.,3.,4.,5.,5.,4.,4.,5.,5.,4.,3.,4.,4.,3.])
        end

        m = @agent(
            2,

            [intLocal,int]::Local,

            Equation = begin
                d(l2) = 0.
            end,

            UpdateInteraction = begin
                if sqrt((x.i - x.j)^2+(y.i - y.j)^2) <= 1.1
                    int.i += 1
                end
            end,

            UpdateLocalInteraction = begin
                intLocal.i += 1
            end,

            Boundary = BoundaryFlat(2,Free(),Periodic())
        )

        mo = compile(m,platform=platform,neighbors="grid",debug=false)
        # println(mo.program)
        @test begin
            com = Community(mo,N=16)
            com.x .= [-1.5,-.5,.5,1.5,-1.5,-.5,.5,1.5,-1.5,-.5,.5,1.5,-1.5,-.5,.5,1.5]
            com.y .= [-1.5,-1.5,-1.5,-1.5,-.5,-.5,-.5,-.5,.5,.5,.5,.5,1.5,1.5,1.5,1.5]
            com.simulationBox = [-2. 2.;-2. 2.]
            com.radiusInteraction = 1.1

            comt = mo.evolve(com,dt=0.1,tMax=1)
            (comt[1].intLocal == [6.,9.,9.,6.,6.,9.,9.,6.,6.,9.,9.,6.,6.,9.,9.,6.]) && (comt[1].int == [3.,4.,4.,3.,4.,5.,5.,4.,4.,5.,5.,4.,3.,4.,4.,3.])
        end

        m = @agent(
            2,

            [intLocal,int]::Local,

            Equation = begin
                d(l2) = 0.
            end,

            UpdateInteraction = begin
                if sqrt((x.i - x.j)^2+(y.i - y.j)^2) <= 1.1
                    int.i += 1
                end
            end,

            UpdateLocalInteraction = begin
                intLocal.i += 1
            end,

            Boundary = BoundaryFlat(2,Periodic(),Periodic())
        )

        mo = compile(m,platform=platform,neighbors="grid",debug=false)
        # println(mo.program)
        @test begin
            com = Community(mo,N=16)
            com.x .= [-1.5,-.5,.5,1.5,-1.5,-.5,.5,1.5,-1.5,-.5,.5,1.5,-1.5,-.5,.5,1.5]
            com.y .= [-1.5,-1.5,-1.5,-1.5,-.5,-.5,-.5,-.5,.5,.5,.5,.5,1.5,1.5,1.5,1.5]
            com.simulationBox = [-2. 2.;-2. 2.]
            com.radiusInteraction = 1.1

            comt = mo.evolve(com,dt=0.1,tMax=1)
            (comt[1].intLocal == [9.,9.,9.,9.,9.,9.,9.,9.,9.,9.,9.,9.,9.,9.,9.,9.]) && (comt[1].int == [3.,4.,4.,3.,4.,5.,5.,4.,4.,5.,5.,4.,3.,4.,4.,3.])
        end

        #Check 3D Periodic
        m = @agent(
            3,

            [intLocal,int]::Local,

            Equation = begin
                d(l2) = 0.
            end,

            UpdateInteraction = begin
                if sqrt((x.i - x.j)^2+(y.i - y.j)^2+(z.i - z.j)^2) <= 1.1
                    int.i += 1
                end
            end,

            UpdateLocalInteraction = begin
                intLocal.i += 1
            end,

            Boundary = BoundaryFlat(3,Periodic(),Free(),Free())
        )

        mo = compile(m,platform=platform,neighbors="grid",debug=false)
        # println(mo.program)
        @test begin
            com = Community(mo,N=64)
            intLocal = zeros(64)
            intt = zeros(64)
            for (i,x) in enumerate([-1.5,-.5,.5,1.5])
                for (j,y) in enumerate([-1.5,-.5,.5,1.5])
                    for (k,z) in enumerate([-1.5,-.5,.5,1.5])
                        com.x[i+4*(j-1)+16*(k-1)] = x
                        com.y[i+4*(j-1)+16*(k-1)] = y
                        com.z[i+4*(j-1)+16*(k-1)] = z

                        intLocal[i+4*(j-1)+16*(k-1)] = 27.
                        intt[i+4*(j-1)+16*(k-1)] = 7.

                        if (abs(x) == 1.5 && abs(y) == 1.5) || (abs(x) == 1.5 && abs(z) == 1.5) || abs(y) == 1.5 || abs(z) == 1.5
                            intLocal[i+4*(j-1)+16*(k-1)] = 18.
                        end
                        if (abs(y) == 1.5 && abs(z) == 1.5)
                            intLocal[i+4*(j-1)+16*(k-1)] = 12.
                        end

                        if abs(x) == 1.5 || abs(y) == 1.5 || abs(z) == 1.5
                            intt[i+4*(j-1)+16*(k-1)] = 6.
                        end
                        if (abs(x) == 1.5 && abs(y) == 1.5) || (abs(x) == 1.5 && abs(z) == 1.5) || (abs(y) == 1.5 && abs(z) == 1.5)
                            intt[i+4*(j-1)+16*(k-1)] = 5.
                        end
                        if abs(x) == 1.5 && abs(y) == 1.5 && abs(z) == 1.5
                            intt[i+4*(j-1)+16*(k-1)] = 4.
                        end
                    end
                end
            end
            com.simulationBox = [-2. 2.;-2. 2.;-2. 2]
            com.radiusInteraction = 1.1
            comt = mo.evolve(com,dt=0.1,tMax=1)
            (comt[1].intLocal == intLocal) && (comt[1].int == intt)
        end

        m = @agent(
            3,

            [intLocal,int]::Local,

            Equation = begin
                d(l2) = 0.
            end,

            UpdateInteraction = begin
                if sqrt((x.i - x.j)^2+(y.i - y.j)^2+(z.i - z.j)^2) <= 1.1
                    int.i += 1
                end
            end,

            UpdateLocalInteraction = begin
                intLocal.i += 1
            end,

            Boundary = BoundaryFlat(3,Free(),Periodic(),Free())
        )

        mo = compile(m,platform=platform,neighbors="grid",debug=false)
        # println(mo.program)
        @test begin
            com = Community(mo,N=64)
            intLocal = zeros(64)
            intt = zeros(64)
            for (i,x) in enumerate([-1.5,-.5,.5,1.5])
                for (j,y) in enumerate([-1.5,-.5,.5,1.5])
                    for (k,z) in enumerate([-1.5,-.5,.5,1.5])
                        com.x[i+4*(j-1)+16*(k-1)] = x
                        com.y[i+4*(j-1)+16*(k-1)] = y
                        com.z[i+4*(j-1)+16*(k-1)] = z

                        intLocal[i+4*(j-1)+16*(k-1)] = 27.
                        intt[i+4*(j-1)+16*(k-1)] = 7.

                        if (abs(x) == 1.5 && abs(y) == 1.5) || (abs(y) == 1.5 && abs(z) == 1.5) || abs(x) == 1.5 || abs(z) == 1.5
                            intLocal[i+4*(j-1)+16*(k-1)] = 18.
                        end
                        if (abs(x) == 1.5 && abs(z) == 1.5)
                            intLocal[i+4*(j-1)+16*(k-1)] = 12.
                        end

                        if abs(x) == 1.5 || abs(y) == 1.5 || abs(z) == 1.5
                            intt[i+4*(j-1)+16*(k-1)] = 6.
                        end
                        if (abs(x) == 1.5 && abs(y) == 1.5) || (abs(x) == 1.5 && abs(z) == 1.5) || (abs(y) == 1.5 && abs(z) == 1.5)
                            intt[i+4*(j-1)+16*(k-1)] = 5.
                        end
                        if abs(x) == 1.5 && abs(y) == 1.5 && abs(z) == 1.5
                            intt[i+4*(j-1)+16*(k-1)] = 4.
                        end
                    end
                end
            end
            com.simulationBox = [-2. 2.;-2. 2.;-2. 2]
            com.radiusInteraction = 1.1
            comt = mo.evolve(com,dt=0.1,tMax=1)
            (comt[1].intLocal == intLocal) && (comt[1].int == intt)
        end

        m = @agent(
            3,

            [intLocal,int]::Local,

            Equation = begin
                d(l2) = 0.
            end,

            UpdateInteraction = begin
                if sqrt((x.i - x.j)^2+(y.i - y.j)^2+(z.i - z.j)^2) <= 1.1
                    int.i += 1
                end
            end,

            UpdateLocalInteraction = begin
                intLocal.i += 1
            end,

            Boundary = BoundaryFlat(3,Free(),Free(),Periodic())
        )

        mo = compile(m,platform=platform,neighbors="grid",debug=false)
        # println(mo.program)
        @test begin
            com = Community(mo,N=64)
            intLocal = zeros(64)
            intt = zeros(64)
            for (i,x) in enumerate([-1.5,-.5,.5,1.5])
                for (j,y) in enumerate([-1.5,-.5,.5,1.5])
                    for (k,z) in enumerate([-1.5,-.5,.5,1.5])
                        com.x[i+4*(j-1)+16*(k-1)] = x
                        com.y[i+4*(j-1)+16*(k-1)] = y
                        com.z[i+4*(j-1)+16*(k-1)] = z

                        intLocal[i+4*(j-1)+16*(k-1)] = 27.
                        intt[i+4*(j-1)+16*(k-1)] = 7.

                        if (abs(x) == 1.5 && abs(z) == 1.5) || (abs(y) == 1.5 && abs(z) == 1.5) || abs(x) == 1.5 || abs(y) == 1.5
                            intLocal[i+4*(j-1)+16*(k-1)] = 18.
                        end
                        if (abs(x) == 1.5 && abs(y) == 1.5)
                            intLocal[i+4*(j-1)+16*(k-1)] = 12.
                        end

                        if abs(x) == 1.5 || abs(y) == 1.5 || abs(z) == 1.5
                            intt[i+4*(j-1)+16*(k-1)] = 6.
                        end
                        if (abs(x) == 1.5 && abs(y) == 1.5) || (abs(x) == 1.5 && abs(z) == 1.5) || (abs(y) == 1.5 && abs(z) == 1.5)
                            intt[i+4*(j-1)+16*(k-1)] = 5.
                        end
                        if abs(x) == 1.5 && abs(y) == 1.5 && abs(z) == 1.5
                            intt[i+4*(j-1)+16*(k-1)] = 4.
                        end
                    end
                end
            end
            com.simulationBox = [-2. 2.;-2. 2.;-2. 2]
            com.radiusInteraction = 1.1
            comt = mo.evolve(com,dt=0.1,tMax=1)
            (comt[1].intLocal == intLocal) && (comt[1].int == intt)
        end

        m = @agent(
            3,

            [intLocal,int]::Local,

            Equation = begin
                d(l2) = 0.
            end,

            UpdateInteraction = begin
                if sqrt((x.i - x.j)^2+(y.i - y.j)^2+(z.i - z.j)^2) <= 1.1
                    int.i += 1
                end
            end,

            UpdateLocalInteraction = begin
                intLocal.i += 1
            end,

            Boundary = BoundaryFlat(3,Periodic(),Periodic(),Periodic())
        )

        mo = compile(m,platform=platform,neighbors="grid",debug=false)
        # println(mo.program)
        @test begin
            com = Community(mo,N=64)
            intLocal = zeros(64)
            intt = zeros(64)
            for (i,x) in enumerate([-1.5,-.5,.5,1.5])
                for (j,y) in enumerate([-1.5,-.5,.5,1.5])
                    for (k,z) in enumerate([-1.5,-.5,.5,1.5])
                        com.x[i+4*(j-1)+16*(k-1)] = x
                        com.y[i+4*(j-1)+16*(k-1)] = y
                        com.z[i+4*(j-1)+16*(k-1)] = z

                        intLocal[i+4*(j-1)+16*(k-1)] = 27.
                        intt[i+4*(j-1)+16*(k-1)] = 7.

                        if abs(x) == 1.5 || abs(y) == 1.5 || abs(z) == 1.5
                            intt[i+4*(j-1)+16*(k-1)] = 6.
                        end
                        if (abs(x) == 1.5 && abs(y) == 1.5) || (abs(x) == 1.5 && abs(z) == 1.5) || (abs(y) == 1.5 && abs(z) == 1.5)
                            intt[i+4*(j-1)+16*(k-1)] = 5.
                        end
                        if abs(x) == 1.5 && abs(y) == 1.5 && abs(z) == 1.5
                            intt[i+4*(j-1)+16*(k-1)] = 4.
                        end
                    end
                end
            end
            com.simulationBox = [-2. 2.;-2. 2.;-2. 2]
            com.radiusInteraction = 1.1
            comt = mo.evolve(com,dt=0.1,tMax=1)
            (comt[1].intLocal == intLocal) && (comt[1].int == intt)
        end

    end

end


text:
# Autogenerated wrapper script for JSBSim_jll for x86_64-apple-darwin
export libJSBSimJL

using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("JSBSim")
JLLWrappers.@declare_library_product(libJSBSimJL, "@rpath/libJSBSimJL.dylib")
function __init__()
    JLLWrappers.@generate_init_header(libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        libJSBSimJL,
        "lib/libJSBSimJL.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
using CombinatorialBandits
using Test

using LightGraphs
using Random
using Distributions

is_travis = "TRAVIS_JULIA_VERSION" in keys(ENV) || true
# TODO: need for optional dependencies for this to work, I suppose.
# https://github.com/JuliaLang/Pkg.jl/issues/1285
if ! is_travis
  using JuMP
  # Why Gurobi?
  # - Need support for lazy constraints (elementary paths):
  #   - Gurobi.jl does not seem to work with lazy constraints
  #   - Mosek does not support them
  # - Need support for MISOCP
  #   - Pajarito is not yet ported to MOI
  using Gurobi
end

function test_policy_trace(instance::CombinatorialInstance{T}, policy::Policy,
                           n_rounds::Int, total_arm_count::Int,
                           final_regret::Float64, final_reward::Float64) where T
  Random.seed!(1)
  s, t = simulate(instance, policy, n_rounds, with_trace=true)

  # State should still be consistent.
  @test s.round == n_rounds
  @test sum(values(s.arm_counts)) <= total_arm_count
  @test s.regret ≈ final_regret atol=1.e-9
  @test s.reward ≈ final_reward
  @test s.regret + s.reward ≈ final_regret + final_reward

  # Now, the trace.
  @test length(t.states) == n_rounds
  @test length(t.arms) == n_rounds
  @test length(t.reward) == n_rounds
  @test length(t.policy_details) == n_rounds

  @test t.states[end].regret ≈ s.regret
  @test t.states[end].reward ≈ s.reward

  for i in 1:n_rounds
    @test t.states[i].round == i
    if i > 1
      @test t.states[i].reward >= t.states[i - 1].reward

      for arm in keys(t.states[i].arm_counts)
        @test t.states[i].arm_counts[arm] >= t.states[i - 1].arm_counts[arm]
        @test t.states[i].arm_reward[arm] >= t.states[i - 1].arm_reward[arm]
      end

      @test is_feasible(instance, t.arms[i])
    end

    for arm in keys(t.states[i].arm_counts)
      if t.states[i].arm_counts[arm] > 0
        @test t.states[i].arm_average_reward[arm] ≈ t.states[i].arm_reward[arm] / t.states[i].arm_counts[arm]
      end
    end
  end
end

import CombinatorialBandits: solve_linear

struct PerfectBipartiteMatchingNoSolver <: PerfectBipartiteMatchingSolver end
solve_linear(::PerfectBipartiteMatchingNoSolver, ::Dict{Tuple{Int64, Int64}, Float64}) = Tuple{Int, Int}[(1, 1)]

struct ElementaryPathNoSolver <: ElementaryPathSolver end
solve_linear(::ElementaryPathNoSolver, ::Dict{Tuple{Int64, Int64}, Float64}) = Tuple{Int, Int}[(1, 2)]

struct SpanningTreeNoSolver <: SpanningTreeSolver end
solve_linear(::SpanningTreeNoSolver, ::Dict{Tuple{Int64, Int64}, Float64}) = Tuple{Int, Int}[(1, 2)]

@testset "CombinatorialBandits.jl" begin
  @testset "Instances" begin
    include("instances_elementarypath.jl")
    include("instances_matching.jl")
    include("instances_spanningtree.jl")
  end

  @testset "Policies" begin
    include("policies_ts.jl")
    include("policies_llr.jl")
    include("policies_cucb.jl")
    include("policies_escb2.jl")

    if ! is_travis
      include("policies_graves_lai.jl")
    end
  end

  include("main.jl")
end


text:
"""
A MultiplicationOperator is defined by a (matrix-like) object that multiplies
coefficients. The multiplication is in-place if type parameter INPLACE is true,
otherwise it is not in-place.
"""
struct MultiplicationOperator{T,ARRAY,INPLACE} <: DictionaryOperator{T}
    src     ::  Dictionary
    dest    ::  Dictionary
    object  ::  ARRAY

    function MultiplicationOperator{T,ARRAY,INPLACE}(src, dest, object) where {T,ARRAY,INPLACE}
        # @assert size(object,1) == length(dest)
        # @assert size(object,2) == length(src)
        new(src, dest, object)
    end
end

object(op::MultiplicationOperator) = op.object

MultiplicationOperator(args...; options...) =
    MultiplicationOperator{deduce_eltype(args...)}(args...; options...)

MultiplicationOperator{T}(src::Dictionary, dest::Dictionary, object; inplace = false) where {T} =
    MultiplicationOperator{T,typeof(object),inplace}(src, dest, object)

MultiplicationOperator(matrix::AbstractMatrix{T}) where {T <: Number} =
    MultiplicationOperator(DiscreteVectorDictionary{T}(size(matrix, 2)), DiscreteVectorDictionary{T}(size(matrix, 1)), matrix)

similar_operator(op::MultiplicationOperator{S,ARRAY,INPLACE}, src, dest) where {S,ARRAY,INPLACE} =
    MultiplicationOperator(src, dest, object(op); inplace=INPLACE)

unsafe_wrap_operator(src, dest, op::MultiplicationOperator) =
    similar_operator(op, src, dest)

isinplace(op::MultiplicationOperator{T,ARRAY,INPLACE}) where {T,ARRAY,INPLACE} = INPLACE

# We pass the object as an additional variable so we can dispatch on it.
# We only intercept non-inplace operators, as the in-place operators are
# intercept higher up where apply_inplace! ends up being called instead
apply!(op::MultiplicationOperator{T,ARRAY,false}, coef_dest, coef_src) where {T,ARRAY} =
    _apply!(op, coef_dest, coef_src, op.object)

# General definition
function _apply!(op::MultiplicationOperator{T,ARRAY,false}, coef_dest, coef_src, object) where {T,ARRAY}
    # Note: this is very likely to allocate memory in the right hand side
    coef_dest[:] = object * coef_src
    coef_dest
end

# Definition in terms of mul! for vectors
_apply!(op::MultiplicationOperator{ELT2,Array{ELT1,2},false}, coef_dest::AbstractVector{T}, coef_src::AbstractVector{T}, object) where {ELT1,ELT2,T} =
    mul!(coef_dest, object, coef_src)

# Be forgiving for matrices: if the coefficients are multi-dimensional, reshape to a linear array first.
_apply!(op::MultiplicationOperator{ELT2,Array{ELT1,2},false}, coef_dest::AbstractArray{T,N1}, coef_src::AbstractArray{T,N2}, object) where {ELT1,ELT2,T,N1,N2} =
    apply!(op, reshape(coef_dest, length(coef_dest)), reshape(coef_src, length(coef_src)))

# In-place definition
# We pass the object as an additional variable so we can dispatch on it
apply_inplace!(op::MultiplicationOperator{T,ARRAY,true}, coef_srcdest) where {T,ARRAY} =
    _apply_inplace!(op, coef_srcdest, op.object)

_apply_inplace!(op::MultiplicationOperator{T,ARRAY,true}, coef_srcdest, object) where {T,ARRAY} =
    object * coef_srcdest

# TODO: introduce unsafe_matrix here, and make a copy otherwise
# matrix(op::MatrixOperator) = op.object
#
# matrix!(op::MatrixOperator, a::Array) = (a[:] = op.object)

adjoint(op::MultiplicationOperator) = adjoint_multiplication(op, object(op))
# This can be overriden for types of objects that do not support adjoint
adjoint_multiplication(op::MultiplicationOperator, object) =
    MultiplicationOperator(dest(op), src(op), adjoint(object))

adjoint_multiplication(op::MultiplicationOperator, object::Array{T,2}) where {T} =
    # We copy the adjoint in order to avoid storing an Adjoint type
    MultiplicationOperator(dest(op), src(op), copy(adjoint(object)))

inv(op::MultiplicationOperator) = inv_multiplication(op, object(op))

# This can be overriden for types of objects that do not support inv
inv_multiplication(op::MultiplicationOperator, object) = MultiplicationOperator(dest(op), src(op), inv(object))

# Use QR for matrices by default
# inv_multiplication(op::MatrixOperator, matrix) = QR_solver(op)

"""
A FunctionOperator applies a given function to the set of coefficients and
returns the result.
"""
struct FunctionOperator{T,F} <: DictionaryOperator{T}
    src     ::  Dictionary
    dest    ::  Dictionary
    fun     ::  F
end

FunctionOperator(args...) = FunctionOperator{deduce_eltype(args...)}(args...)

FunctionOperator{T}(src, dest, fun) where {T} =
    FunctionOperator{T,typeof(fun)}(src, dest, fun)

similar_operator(op::FunctionOperator{T,F}, src, dest) where {T,F} =
    FunctionOperator{promote_type(T,operatoreltype(src,dest)),F}(src, dest, op.fun)

unsafe_wrap_operator(src, dest, op::FunctionOperator) =
    similar_operator(op, src, dest)

# Warning: this very likely allocates memory
apply!(op::FunctionOperator, coef_dest, coef_src) = apply_fun!(op, op.fun, coef_dest, coef_src)

function apply_fun!(op::FunctionOperator, fun, coef_dest, coef_src)
    coef_dest[:] = fun(coef_src)
end

adjoint(op::FunctionOperator) = adjoint_function(op, op.fun)
# This can be overriden for types of functions that do not support adjoint
adjoint_function(op::FunctionOperator, fun) =
    FunctionOperator(dest(op), src(op), adjoint(fun))

inv(op::FunctionOperator) = inv_function(op, op.fun)

# This can be overriden for types of functions that do not support inv
inv_function(op::FunctionOperator, fun) = FunctionOperator(dest(op), src(op), inv(fun))

string(op::FunctionOperator) = "Function " * string(op.fun)

"A linear combination of operators: val1 * op1 + val2 * op2."
struct OperatorSum{T,OP1 <: DictionaryOperator,OP2 <: DictionaryOperator,S} <: DictionaryOperator{T}
    op1         ::  OP1
    op2         ::  OP2
    val1        ::  T
    val2        ::  T
    scratch     ::  S

    function OperatorSum{T,OP1,OP2,S}(op1::OP1, op2::OP2, val1::T, val2::T, scratch::S) where {OP1 <: DictionaryOperator,OP2 <: DictionaryOperator,T,S}
        # We don't enforce that source and destination of op1 and op2 are the same, but at least
        # their sizes and coefficient types must match.
        @assert size(src(op1)) == size(src(op2))
        @assert size(dest(op1)) == size(dest(op2))
        @assert coefficienttype(src(op1)) == coefficienttype(src(op2))
        @assert coefficienttype(dest(op1)) == coefficienttype(dest(op2))

        new(op1, op2, val1, val2, scratch)
    end
end

function OperatorSum{T}(op1::DictionaryOperator, op2::DictionaryOperator, val1::Number, val2::Number) where {T}
    scratch = zeros(T, dest(op1))
    OperatorSum{T,typeof(op1),typeof(op2),typeof(scratch)}(op1, op2, T(val1), T(val2), scratch)
end

OperatorSum(op1::DictionaryOperator, op2::DictionaryOperator, val1::Number, val2::Number;
            T=promote_type(eltype(op1), eltype(op2), typeof(val1), typeof(val2))) =
    OperatorSum{T}(op1, op2, val1, val2)

src(op::OperatorSum) = src(op.op1)

dest(op::OperatorSum) = dest(op.op1)

adjoint(op::OperatorSum) = OperatorSum(adjoint(op.op1), adjoint(op.op2), conj(op.val1), conj(op.val2))

pinv(op::OperatorSum) = OperatorSum(pinv(op.op1), pinv(op.op2), op.val1, op.val2)

isdiag(op::OperatorSum) = isdiag(op.op1) && isdiag(op.op2)

apply_inplace!(op::OperatorSum, coef_srcdest) =
    apply_sum_inplace!(op, op.op1, op.op2, coef_srcdest)

function apply_sum_inplace!(op::OperatorSum, op1::DictionaryOperator, op2::DictionaryOperator, coef_srcdest)
    scratch = op.scratch

    apply!(op1, scratch, coef_srcdest)
    apply!(op2, coef_srcdest)

    for i in eachindex(coef_srcdest)
        coef_srcdest[i] = op.val1 * scratch[i] + op.val2 * coef_srcdest[i]
    end
    coef_srcdest
end

apply!(op::OperatorSum, coef_dest, coef_src) = apply_sum!(op, op.op1, op.op2, coef_dest, coef_src)

function apply_sum!(op::OperatorSum, op1::DictionaryOperator, op2::DictionaryOperator, coef_dest, coef_src)
    scratch = op.scratch

    apply!(op1, scratch, coef_src)
    apply!(op2, coef_dest, coef_src)

    for i in eachindex(coef_dest)
        coef_dest[i] = op.val1 * scratch[i] + op.val2 * coef_dest[i]
    end
    coef_dest
end

components(op::OperatorSum) = (op.op1,op.op2)

(+)(op1::DictionaryOperator, op2::DictionaryOperator) = OperatorSum(op1, op2, 1, 1)
(-)(op1::DictionaryOperator, op2::DictionaryOperator) = OperatorSum(op1, op2, 1, -1)

(+)(I1::UniformScaling, op2::DictionaryOperator) = (+)(ScalingOperator(I1.λ, src(op2), dest(op2)), op2)
(-)(I1::UniformScaling, op2::DictionaryOperator) = (-)(ScalingOperator(I1.λ, src(op2), dest(op2)), op2)
(+)(op1::DictionaryOperator, I2::UniformScaling) = (+)(op1, ScalingOperator(I2.λ, src(op1), dest(op1)))
(-)(op1::DictionaryOperator, I2::UniformScaling) = (-)(op1, ScalingOperator(I2.λ, src(op1), dest(op1)))

function stencilarray(op::OperatorSum)
    local s1, s2
    if op.val1 == -1
        s1 = "-"
    elseif op.val1 == 1
        s1 = ""
    else
        s1 = op1.val1
    end
    if op.val2 == -1
        s2 = " - "
    elseif op.val2 == 1
        s2 = " + "
    else
        s2 = " + " * string(op.val2) * " "
    end
    [s1, op.op1, s2, op.op2]
end

stencil_parentheses(op::OperatorSum) = true
object_parentheses(op::OperatorSum) = true

"""
An operator that calls linearize on a native representation of a set, returning
a Vector with the length of the set. For example, one can not apply a matrix
to a non-arraylike representation, hence the representation has to be linearized
first.
"""
struct LinearizationOperator{T} <: DictionaryOperator{T}
    src         ::  Dictionary
    dest        ::  Dictionary
end

LinearizationOperator(dicts::Dictionary...) =
    LinearizationOperator{operatoreltype(dicts...)}(dicts...)

LinearizationOperator{T}(src::Dictionary) where {T} =
    LinearizationOperator{T}(src, DiscreteVectorDictionary{coefficienttype(src)}(length(src)))

similar_operator(::LinearizationOperator{T}, src, dest) where {T} =
    LinearizationOperator{promote_type(T,operatoreltype(src,dest))}(src, dest)

apply!(op::LinearizationOperator, coef_dest, coef_src) =
    linearize_coefficients!(coef_dest, coef_src)

isdiag(op::LinearizationOperator) = true

Base.adjoint(op::LinearizationOperator{T}) where {T} =
    DelinearizationOperator{T}(dest(op), src(op))

"The inverse of a LinearizationOperator."
struct DelinearizationOperator{T} <: DictionaryOperator{T}
    src         ::  Dictionary
    dest        ::  Dictionary
end

DelinearizationOperator(dicts::Dictionary...) =
    DelinearizationOperator{operatoreltype(dicts...)}(dicts...)

DelinearizationOperator{T}(dest::Dictionary) where {T} =
    DelinearizationOperator{T}(DiscreteVectorDictionary{T}(length(dest)), dest)

similar_operator(::DelinearizationOperator{T}, src, dest) where {T} =
    DelinearizationOperator{promote_type(T,operatoreltype(src,dest))}(src, dest)

apply!(op::DelinearizationOperator, coef_dest, coef_src) =
    delinearize_coefficients!(coef_dest, coef_src)

isdiag(op::DelinearizationOperator) = true

Base.adjoint(op::DelinearizationOperator{T}) where {T} =
    LinearizationOperator{T}(dest(op), src(op))

const AlternatingSignOperator{T} = DiagonalOperator{T,AlternatingSigns{T}}

AlternatingSignOperator(src::Dictionary) = AlternatingSignOperator{operatoreltype(src)}(src)

function AlternatingSignOperator{T}(src::Dictionary) where {T}
    diag = Diagonal(AlternatingSigns{T}(length(src)))
    DiagonalOperator{T}(diag, src, src)
end

strings(op::AlternatingSignOperator) = tuple("Alternating sign operator of length $(size(op,1))")

const CoefficientScalingOperator{T} = DiagonalOperator{T,ScaledEntry{T}}

CoefficientScalingOperator(src::Dictionary, index::Int, scalar) =
	CoefficientScalingOperator{promote_type(typeof(scalar),operatoreltype(src))}(src, index, scalar)

function CoefficientScalingOperator{T}(src::Dictionary, index::Int, scalar) where {T}
    diag = Diagonal(ScaledEntry{T}(length(src), index, scalar))
	DiagonalOperator{T}(diag, src, src)
end

strings(op::CoefficientScalingOperator) =
    tuple("Diagonal operator of length $(size(op,1)) that scales coefficient $(op.A.diag.index) by $(op.A.diag.scalar)")


text:
using Makie
import Makie.plot
using Statistics
using SparseArrays
function plot(X::Unfold.DesignMatrix;standardize=true,sort=false)

    designmat = Unfold.get_Xs(X);
    if standardize
        designmat = designmat ./ std(designmat,dims=1)
        designmat[isinf.(designmat)] .= 1.
    end
    if sort
        designmat = Base.sortslices(designmat,dims=1)
    end
    labels = Unfold.get_coefnames(X)

    if isa(designmat, SparseMatrixCSC)
        @assert(!sort,"Sorting does not make sense for timeexpanded designmatrices")
        designmat = Matrix(designmat[end÷2-2000:end÷2+2000,:])
    end
    # plot Designmatrix
    fig, ax, hm = heatmap(designmat',axis=(xticks=(1:length(labels),labels),xticklabelrotation = pi/8),)
    Colorbar(fig[1,2],hm)

    if isa(designmat, SparseMatrixCSC)
        ax.yreversed = true
    end
    return fig
end


text:
##
##!------------------------------------------------------------------------------------------------!
##  graph.jl
##
##  Graph (undirected) types and associated functions.
##  All functions except those with a "ß_" prefix are unsafe.
##
##  Creation date: 16/08/15
##  Author: Arindam Biswas <arindam dot b at eml dot cc>
##!------------------------------------------------------------------------------------------------!
##

# The base graph type.
abstract Graph

abstract BareGraph

type GraphException <: Exception
end

type GraphDynamic <: Graph
    order::UInt32
    size::UInt64
    adjacencies::Vector{Vector{UInt32}}
    vertex_info::Vector{Int64}
    edge_info::Vector{Vector{Int64}}
    vertex_indicator::Vector{Bool} # vertex_indicator[u] is true if u is in the graph.
end

type GraphStatic <: Graph
    order::UInt32
    size::UInt64
    adjacencies::Vector{Vector{UInt32}}
    vertex_info::Vector{Int64}
    edge_info::Vector{Vector{Int64}}
end

type BareGraphDynamic <: BareGraph
    order::UInt32
    size::UInt64
    adjacencies::Vector{Vector{UInt32}}
    vertex_indicator::Vector{Bool} # vertex_indicator[u] is true if u is in the graph.
end

type BareGraphStatic <: BareGraph
    order::UInt32
    size::UInt64
    adjacencies::Vector{Vector{UInt32}}
end

function GraphDynamic(order::UInt32)
    adjacencies = Vector{Vector{UInt32}}(order)
    vertex_info = Vector{Int64}(order)
    edge_info = Vector{Vector{Int64}}(order)
    vertex_indicator = Vector{Bool}(order)

    for i in 1:order
        adjacencies[i] = Vector{UInt32}()
        edge_info[i] = Vector{Int64}()
        vertex_indicator[i] = true
    end
    return GraphDynamic(order, UInt64(0), adjacencies, vertex_info, edge_info, vertex_indicator)
end

function GraphStatic(order::UInt32)
    adjacencies = Vector{Vector{UInt32}}(order)
    vertex_info = Vector{Int64}(order)
    edge_info = Vector{Vector{Int64}}(order)

    for i in 1:order
        adjacencies[i] = Vector{UInt32}()
        edge_info[i] = Vector{Int64}()
    end
    return GraphStatic(order, UInt64(0), adjacencies, vertex_info, edge_info)
end

function BareGraphDynamic(order::UInt32)
    adjacencies = Vector{Vector{UInt32}}(order)
    vertex_indicator = Vector{Bool}(order)

    for i in 1:order
        adjacencies[i] = Vector{UInt32}()
        vertex_indicator[i] = true
    end
    return BareGraphDynamic(order, UInt64(0), adjacencies, vertex_indicator)
end

function BareGraphStatic(order::UInt32)
    adjacencies = Vector{Vector{UInt32}}(order)

    for i in 1:order
        adjacencies[i] = Vector{UInt32}()
    end
    return BareGraphStatic(order, UInt64(0), adjacencies)
end

## Properties
order(g::Union{Graph, BareGraph}) = g.order
size(g::Union{Graph, BareGraph}) = g.size
degree(g::Union{Graph, BareGraph}, u::UInt32) = length(g.adjacencies[u])

## Access
function vertices(g::Union{GraphDynamic, BareGraphDynamic})
    vertices = Vector{UInt32}(g.order)
    # order keeps changing; use length of vertex_indicator.
    i, j = length(g.vertex_indicator), g.order
    while i >= 1
        if g.vertex_indicator[i]
            vertices[j] = i
            j -= 1
        end
        i -= 1
    end
    return vertices
end

function vertices(g::Union{GraphStatic, BareGraphStatic})
    return Vector{UInt32}(1:g.order)
end

function edges(g::Union{GraphDynamic, BareGraphDynamic})
    edges = Vector{Tuple{UInt32, UInt32}}(g.size)
    j = 1
    l = length(g.vertex_indicator)
    for i in 1:l
        if !g.vertex_indicator[i]
            continue
        end
        for t in g.adjacencies[i]
            if i < t
                edges[j] = (i, t)
                j += 1
            end
        end
    end
    return edges
end

function edges(g::Union{GraphStatic, BareGraphStatic})
    edges = Vector{Tuple{UInt32, UInt32}}(g.size)
    j = 1
    for i in 1:g.order
        for t in g.adjacencies[i]
            if i < t
                edges[j] = (i, t)
                j += 1
            end
        end
    end
    return edges
end

function exists(g::Union{GraphDynamic, BareGraphDynamic}, u::UInt32)
    return g.vertex_indicator[u]
end

function exists(g::Union{GraphStatic, BareGraphStatic}, u::UInt32)
    return 0 < u <= g.order
end

# Finds the index of u in l or returns 0, if u is not in l.
function _find_index(l::Vector{UInt32}, u::UInt32)
    for i in 1:length(l)
        if u == l[i]
            return i
        end
    end
    return 0
end

function adjacent(g::Union{GraphDynamic, BareGraphDynamic}, u::UInt32, v::UInt32)
    return _find_index(g.adjacencies[u], v) != 0
end

# Does a binary search on l[i:j] and returns the index if u is found, or 0 if not.
function _binary_search(l::Vector{UInt32}, i::UInt32, j::UInt32, u::UInt32)
    # l contains tuples of the form (target, edgeinfo).
    if i == j
        return l[i] == u ? i : 0
    end

    k = uint32((i + j) / 2)
    v = l[k]
    if u < v
        return _binary_search(l, i, k, u)
    elseif u > v
        return _binary_search(l, k + 1, j, u)
    else
        return k
    end
end

function adjacent(g::Union{GraphStatic, BareGraphStatic}, u::UInt32, v::UInt32)
    # For static graphs, the adjacency list is sorted.
    return _binary_search(g.adjacencies[u], 1, length(g.adjacencies[u]), v) != 0
end

function neighbours(g::Union{Graph, BareGraph}, u::UInt32)
    return g.adjacencies[u]
end

function vertex_info(g::Union{GraphDynamic, GraphStatic}, u::UInt32)
    return g.vertex_info[u]
end

function edge_info(g::GraphDynamic, u::UInt32, v::UInt32)
    l = g.edge_info[u]
    return l[_find_index(g.adjacencies[u], v)]
end

function edge_info(g::GraphStatic, u::UInt32, v::UInt32)
    l = g.edge_info[u]
    return l[_binary_search(g.adjacencies[u], 1, length(g.adjacencies[u]), v)]
end

## Mutation
function add_vertex!(g::GraphDynamic)
    # The vertex is added at the end.
    push!(g.adjacencies, Vector{UInt32}())
    push!(g.vertex_info, Int64(0))
    push!(g.edge_info, Vector{Int64}())
    push!(g.vertex_indicator, true)
    g.order += 1
end

function add_vertex!(g::GraphDynamic, u::UInt32)
    # Assume u was previously deleted.
    g.vertex_indicator[u] = true
    g.order += 1
end

function add_vertex!(g::BareGraphDynamic)
    # The vertex is always added at the end.
    push!(g.adjacencies, Vector{UInt32}())
    push!(g.vertex_indicator, true)
    g.order += 1
end

function add_vertex!(g::BareGraphDynamic, u::UInt32)
    # Assume u was previously deleted.
    g.vertex_indicator[u] = true
    g.order += 1
end

function add_edge!(g::GraphDynamic, u::UInt32, v::UInt32)
    push!(g.adjacencies[u], v)
    push!(g.edge_info[u], v)
    # Store duplicates for the "reverse edge" too. neighbours() needs this.
    push!(g.adjacencies[v], u)
    push!(g.edge_info[v], u)
    g.size += 1
end

function add_edge!(g::BareGraphDynamic, u::UInt32, v::UInt32)
    push!(g.adjacencies[u], v)
    # Store a duplicate for the "reverse edge" too. neighbours() needs this.
    push!(g.adjacencies[v], u)
    g.size += 1
end

function remove_vertex!(g::Union{GraphDynamic, BareGraphDynamic}, u::UInt32)
    # Do not actually remove the vertex, only set vertex_indicator[u] to false.
    g.vertex_indicator[u] = false
    # Remove all edges incident with u.
    n = copy(g.adjacencies[u])
    for v in n
        remove_edge!(g, u, v)
    end
    g.order -= 1
end

function remove_edge!(g::GraphDynamic, u::UInt32, v::UInt32)
    n = g.adjacencies[u]
    l = g.edge_info[u]
    # Find the index of v in the adjacency list.
    i = _find_index(n, v)
    n[i] = n[end]
    pop!(n)
    l[i] = l[end]
    pop!(l)

    # Do the same thing for the "reverse edge".
    n = g.adjacencies[v]
    l = g.edge_info[v]
    i = _find_index(n, u)
    n[i] = n[end]
    pop!(n)
    l[i] = l[end]
    pop!(l)
    g.size -= 1
end

function remove_edge!(g::BareGraphDynamic, u::UInt32, v::UInt32)
    n = g.adjacencies[u]
    # Find the index of v in the adjacency list.
    i = _find_index(n, v)
    n[i] = n[end]
    pop!(n)

    # Do the same thing for the "reverse edge".
    n = g.adjacencies[v]
    i = _find_index(n, u)
    n[i] = n[end]
    pop!(n)
    g.size -= 1
end

function vertex_info!(g::Union{GraphDynamic, GraphStatic}, u::UInt32, info::Int64)
    g.vertex_info[u] = info
end

function edge_info!(g::Union{GraphDynamic, GraphStatic}, u::UInt32, v::UInt32, info::Int64)
    i = _find_index(g.adjacencies[u], v)
    g.edge_info[u][i] = info
    # Do the same for the "reverse edge".
    i = _find_index(g.adjacencies[v], u)
    g.edge_info[v][i] = info
end


text:
module CUBLAS

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cuComplex, cuDoubleComplex, libraryPropertyType, cudaDataType
using ..CUDA: libcublas, unsafe_free!, @retry_reclaim, isdebug, @sync

using GPUArrays

using LinearAlgebra

using BFloat16s

using CEnum

# core library
include("libcublas_common.jl")
include("error.jl")
include("libcublas.jl")
include("libcublas_deprecated.jl")

# low-level wrappers
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("linalg.jl")

# thread cache for task-local library handles
const thread_handles = Vector{Union{Nothing,cublasHandle_t}}()
const thread_xt_handles = Vector{Union{Nothing,cublasXtHandle_t}}()

function math_mode!(handle, mode)
    flags = 0

    # https://github.com/facebookresearch/faiss/issues/1385
    if version(handle) > v"11"
        flags = CUBLAS_MATH_DISALLOW_REDUCED_PRECISION_REDUCTION
    end

    flags |= if mode == CUDA.PEDANTIC_MATH
        # prevent use of tensor cores
        if VERSION < v"11"
            CUBLAS_DEFAULT_MATH
        else
            CUBLAS_PEDANTIC_MATH
        end
    elseif mode == CUDA.DEFAULT_MATH
        # use tensor cores, but don't reduce precision
        if VERSION < v"11"
            CUBLAS_TENSOR_OP_MATH
        else
            CUBLAS_DEFAULT_MATH
        end
    elseif mode == CUDA.FAST_MATH
        # we'll additionally select a compute-mode with reduced precision whenever possible
        if VERSION < v"11"
            CUBLAS_TENSOR_OP_MATH
        else
            CUBLAS_TF32_TENSOR_OP_MATH
        end
    end

    cublasSetMathMode(handle, flags)

    return
end

function handle()
    tid = Threads.threadid()
    if @inbounds thread_handles[tid] === nothing
        ctx = context()
        thread_handles[tid] = get!(task_local_storage(), (:CUBLAS, ctx)) do
            handle = cublasCreate()
            finalizer(current_task()) do task
                CUDA.isvalid(ctx) || return
                context!(ctx) do
                    cublasDestroy_v2(handle)
                end
            end

            cublasSetStream_v2(handle, CuStreamPerThread())

            if version(handle) >= v"11.2"
                workspace = CuArray{UInt8}(undef, 4*1024*1024)  # TODO: 256-byte aligned
                cublasSetWorkspace_v2(handle, workspace, sizeof(workspace))
                task_local_storage((:CUBLAS, :workspace), workspace)
            end

            math_mode!(handle, CUDA.math_mode())

            handle
        end
    end
    something(@inbounds thread_handles[tid])
end

function xt_handle()
    tid = Threads.threadid()
    if @inbounds thread_xt_handles[tid] === nothing
        ctx = context()
        thread_xt_handles[tid] = get!(task_local_storage(), (:CUBLASxt, ctx)) do
            handle = cublasXtCreate()
            finalizer(current_task()) do task
                CUDA.isvalid(ctx) || return
                context!(ctx) do
                    cublasXtDestroy(handle)
                end
            end

            # select the devices
            # TODO: this is weird, since we typically use a single device per thread/context
            devs = convert.(Cint, CUDA.devices())
            cublasXtDeviceSelect(handle, length(devs), devs)

            handle
        end
    end
    something(@inbounds thread_xt_handles[tid])
end

function __init__()
    resize!(thread_handles, Threads.nthreads())
    fill!(thread_handles, nothing)

    resize!(thread_xt_handles, Threads.nthreads())
    fill!(thread_xt_handles, nothing)

    CUDA.atdeviceswitch() do
        tid = Threads.threadid()
        thread_handles[tid] = nothing
        thread_xt_handles[tid] = nothing
    end

    CUDA.attaskswitch() do
        tid = Threads.threadid()
        thread_handles[tid] = nothing
        thread_xt_handles[tid] = nothing
    end
end

function log_message(cstr)
    str = unsafe_string(cstr)
    print(str)
    # NOTE: we can't `@debug` these messages, because the logging function is called for
    #       every line... also not sure what the i!/I! prefixes mean (info?)
    return
end

function __runtime_init__()
    # enable library logging when launched with JULIA_DEBUG=CUBLAS
    if isdebug(:init, CUBLAS)
        callback = @cfunction(log_message, Nothing, (Cstring,))
        cublasSetLoggerCallback(callback)
    end
end

end


text:
const IMAGEOFFSET = 16
const LABELOFFSET = 8

const NROWS = 28
const NCOLS = 28

const TRAINIMAGES = joinpath(@__DIR__, "..", "data", "train-images.idx3-ubyte")
const TRAINLABELS = joinpath(@__DIR__, "..", "data", "train-labels.idx1-ubyte")
const TESTIMAGES  = joinpath(@__DIR__, "..", "data", "t10k-images.idx3-ubyte")
const TESTLABELS  = joinpath(@__DIR__, "..", "data", "t10k-labels.idx1-ubyte")

function imageheader(io::IO)
    magic_number = bswap(read(io, UInt32))
    total_items = bswap(read(io, UInt32))
    nrows = bswap(read(io, UInt32))
    ncols = bswap(read(io, UInt32))
    return magic_number, Int(total_items), Int(nrows), Int(ncols)
end

function labelheader(io::IO)
    magic_number = bswap(read(io, UInt32))
    total_items = bswap(read(io, UInt32))
    return magic_number, Int(total_items)
end

function getimage(io::IO)
    img = Array{UInt8}(undef, NCOLS, NROWS)
    for i in 1:NCOLS, j in 1:NROWS
        img[i, j] = read(io, UInt8)
    end
    return img
end

function getimage(io::IO, index::Integer)
    seek(io, IMAGEOFFSET + NROWS * NCOLS * (index - 1))
    return getimage(io)
end

function getlabel(io::IO)
    i = read(io, UInt8)
    v = zeros(10)
    v[i+1] = 1.
    return v
end

function getlabel(io::IO, index::Integer)
    seek(io, LABELOFFSET + (index - 1))
    return getlabel(io)
end

function trainimages(i::Integer)
    getimage(IOBuffer(read(TRAINIMAGES)), i)
end

function testimages(i::Integer)
    getimage(IOBuffer(read(TESTIMAGES)), i)
end

function trainlabels(i::Integer)
    getlabel(IOBuffer(read(TRAINLABELS)), i)
end

function testlabels(i::Integer)
    getlabel(IOBuffer(read(TESTLABELS)), i)
end

function trainimages()
    io = IOBuffer(read(TRAINIMAGES))
    _, nimages, nrows, ncols = imageheader(io)
    [Int.(getimage(io)) for _ in 1:nimages]
end

function testimages()
    io = IOBuffer(read(TESTIMAGES))
    _, nimages, nrows, ncols = imageheader(io)
    [Int.(getimage(io)) for _ in 1:nimages]
end

function trainlabels()
    io = IOBuffer(read(TRAINLABELS))
    _, nlabels = labelheader(io)
    [getlabel(io) for _ in 1:nlabels]
end

function testlabels()
    io = IOBuffer(read(TESTLABELS))
    _, nlabels = labelheader(io)
    [getlabel(io) for _ in 1:nlabels]
end

traindata() = collect(zip(vec.(trainimages()), trainlabels()))
testdata()  = collect(zip(vec.(testimages()), testlabels()))


text:
module Maps # end

using Makie

using FileIO, Colors, JSON

earth_img() = "https://svs.gsfc.nasa.gov/vis/a000000/a002900/a002915/bluemarble-2048.png"
mars_img() = "https://www.jpl.nasa.gov/spaceimages/images/largesize/PIA23518_hires.jpg"

"""
    globe(img = earth_img(), args...)
Display a 3D globe in a Makie plot.

Any extra plotting arguments will be passed-through to Makie.
See [Makie.plot](@ref) arguments for more details.
"""
function globe(image=earth_img(), args...; kwargs...)
    earth = try
        load(download(image))
    catch e
        @warn("Downloading the image failed. Using random image, so this test will fail! (error: $e)")
        rand(RGBAf0, 100, 100) # don't error test when e.g. offline
    end;
    m = GLNormalUVMesh(Sphere(Point3f0(0), 1f0), 60);
    scene = mesh(m, args...; color = earth, shading = false, kwargs...);
    Makie.AbstractPlotting.cam3d_cad!(scene)

    scene
end

"""
    usa(resolution = 0.025, args...)
Display a 2D map of the United States, on an Equirectangular projection.

Any extra plotting arguments will be passed-through to Makie.
See [Makie.plot](@ref) arguments for more details.

Data comes from http://openstreetmap.org.
"""
function usa(resolution = 0.015, args...; kwargs...)
    allregions = usa_polygons()
    scene = plot_regions(allregions, args...; resolution = resolution, kwargs...)
    text!(scene,
          "© OpenStreetMap contributors",
          position = (0,0),
          color = :red,
         )
    scene
end
function usa_polygons()
    USA_json = JSON.parse(open("data/nominatim/USA.json", "r"))
    allregions = USA_json[1]["geojson"]["coordinates"]
    # Flatten the regions (for some reason they're all in a nested vector)
    allregions = [Point2f0.(map.(Float32, c))
                     for r in allregions
                     for c in r]
    sort(allregions, by=length, rev=true)
end

function scale_region(region, resolution)
    scale = (1/resolution)
    orig_size = length(region)
    if scale >= orig_size÷2
        scale = orig_size÷2 - 1
    end

    # Skip every nth element, incrementing counter by floating point
    out_size = Int(ceil(1/scale * orig_size))
    out = similar(region, out_size)
    r_idx, out_idx = 1,1
    while floor(Int, r_idx) <= orig_size
        out[out_idx] = region[floor(Int, r_idx)]
        r_idx += scale
        out_idx += 1
    end
    out
end
function plot_regions(regions, args...; kwargs...)
    scene = Scene()
    for r in regions
        plot_region!(scene, r, args...; kwargs...)
    end
    scene
end
function plot_region!(scene, region::Array{Point2f0}, pltargs...; resolution = 0.01, pltkwargs...)
    scaled = scale_region(region, resolution)
    # NOTE: there is an error in `poly` that it sometimes deletes a point from inputs! D:
    # https://github.com/JuliaPlots/Makie.jl/issues/417
    # So for now, we copy() the input and wrap it in a try/catch
    try
        Makie.poly!(scene, copy(scaled), pltargs...; pltkwargs...)
    catch e
        e isa ErrorException && return scene
        rethrow()
    end
end

end # module


text:
# Autogenerated wrapper script for CUDA_Toolkit_jll for x86_64-w64-mingw32-cuda+10.2
export libcublas, libcudadevrt, libcufft, libcupti, libcurand, libcusolver, libcusolverMg, libcusparse, libdevice, libnvtoolsext, libnvvm, nvdisasm, nvlink, ptxas

JLLWrappers.@generate_wrapper_header("CUDA_Toolkit")
JLLWrappers.@declare_library_product(libcublas, "cublas64_10.dll")
JLLWrappers.@declare_file_product(libcudadevrt)
JLLWrappers.@declare_library_product(libcufft, "cufft64_10.dll")
JLLWrappers.@declare_library_product(libcupti, "cupti64_102.dll")
JLLWrappers.@declare_library_product(libcurand, "curand64_10.dll")
JLLWrappers.@declare_library_product(libcusolver, "cusolver64_10.dll")
JLLWrappers.@declare_library_product(libcusolverMg, "cusolverMg64_10.dll")
JLLWrappers.@declare_library_product(libcusparse, "cusparse64_10.dll")
JLLWrappers.@declare_file_product(libdevice)
JLLWrappers.@declare_library_product(libnvtoolsext, "nvToolsExt64_1.dll")
JLLWrappers.@declare_library_product(libnvvm, "nvvm64_33_0.dll")
JLLWrappers.@declare_executable_product(nvdisasm)
JLLWrappers.@declare_executable_product(nvlink)
JLLWrappers.@declare_executable_product(ptxas)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcublas,
        "bin\\cublas64_10.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        libcudadevrt,
        "lib\\cudadevrt.lib",
    )

    JLLWrappers.@init_library_product(
        libcufft,
        "bin\\cufft64_10.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcupti,
        "bin\\cupti64_102.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcurand,
        "bin\\curand64_10.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcusolver,
        "bin\\cusolver64_10.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcusolverMg,
        "bin\\cusolverMg64_10.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcusparse,
        "bin\\cusparse64_10.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        libdevice,
        "share\\libdevice\\libdevice.10.bc",
    )

    JLLWrappers.@init_library_product(
        libnvtoolsext,
        "bin\\nvToolsExt64_1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnvvm,
        "bin\\nvvm64_33_0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        nvdisasm,
        "bin\\nvdisasm.exe",
    )

    JLLWrappers.@init_executable_product(
        nvlink,
        "bin\\nvlink.exe",
    )

    JLLWrappers.@init_executable_product(
        ptxas,
        "bin\\ptxas.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export CrossNode

import ..Nodes.draw

mutable struct CrossNode <: Ranger.AbstractNode
    base::Nodes.NodeData
    transform::Nodes.TransformProperties{Float64}

    mesh::Geometry.Mesh

    color::Rendering.Palette

    function CrossNode(world::Ranger.World, name::String, parent::Ranger.AbstractNode)
        o = new()

        o.base = Nodes.NodeData(name, parent, world)
        o.transform = Nodes.TransformProperties{Float64}()
        o.mesh = Geometry.Mesh()
        o.color = Rendering.White()

        # horizontal
        Geometry.add_vertex!(o.mesh, -0.5, 0.0)
        Geometry.add_vertex!(o.mesh, 0.5, 0.0)

        # vertical
        Geometry.add_vertex!(o.mesh, 0.0, -0.5)
        Geometry.add_vertex!(o.mesh, 0.0, 0.5)

        Geometry.build!(o.mesh)

        o
    end
end

function Nodes.draw(node::CrossNode, context::Rendering.RenderContext)
    if Nodes.is_dirty(node)
        Rendering.transform!(context, node.mesh)
        Nodes.set_dirty!(node, false)
    end

    Rendering.set_draw_color(context, node.color)

    v1 = node.mesh.bucket[1]
    v2 = node.mesh.bucket[2]
    v3 = node.mesh.bucket[3]
    v4 = node.mesh.bucket[4]
    Rendering.draw_horz_line(context, v1.x, v2.x, v1.y)

    Rendering.draw_vert_line(context, v3.x, v3.y, v4.y);
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
module Rematch

import MacroTools
import MacroTools: @capture
import Compat: occursin, fieldcount

macro splice(iterator, body)
  @assert iterator.head == :call
  @assert iterator.args[1] == :in
  Expr(:..., :(($(esc(body)) for $(esc(iterator.args[2])) in $(esc(iterator.args[3])))))
end

struct MatchFailure
    value
end

function assert_num_fields(::Type{T}, matched::Integer) where {T}
    actual = fieldcount(T)
    @assert actual == matched "Tried to match $matched fields of $actual field struct $T"
end

"""
Statically get the fieldcount of a type. Useful to avoid runtime calls to
fieldcount.
"""
@generated function evaluated_fieldcount(t::Type{T}) where T
    # TODO(nathan.daly) Does this need to be @generated, or is this function even necessary?
    fieldcount(T)
end

function handle_destruct_fields(value::Symbol, pattern, subpatterns, len, get::Symbol, bound::Set{Symbol}, asserts::Vector{Expr}; allow_splat=true)
    # NOTE we assume `len` is cheap
    fields = []
    seen_splat = false
    for (i,subpattern) in enumerate(subpatterns)
        if (subpattern isa Expr) && (subpattern.head == :(...))
            @assert allow_splat && !seen_splat "Too many ... in pattern $pattern"
            @assert length(subpattern.args) == 1
            seen_splat = true
            push!(fields, (:($i:($len-$(length(subpatterns)-i))), subpattern.args[1]))
        elseif seen_splat
            push!(fields, (:($len-$(length(subpatterns)-i)), subpattern))
        elseif (subpattern isa Expr) && (subpattern.head == :kw)
            field_symbol = Meta.quot(subpattern.args[1])
            push!(fields, (field_symbol, subpattern.args[2]))
        else
            push!(fields, (i, subpattern))
        end
    end
    Expr(:&&,
        if seen_splat
            :($len >= $(length(subpatterns)-1))
        else
            :($len == $(length(subpatterns)))
        end,
        @splice (i, (field, subpattern)) in enumerate(fields) quote
            $(Symbol("$(value)_$i")) = $get($value, $field)
            $(handle_destruct(Symbol("$(value)_$i"), subpattern, bound, asserts))
        end)
end

function handle_destruct(value::Symbol, pattern, bound::Set{Symbol}, asserts::Vector{Expr})
    if pattern == :(_)
        # wildcard
        true
    elseif !(pattern isa Expr || pattern isa Symbol) ||
           pattern == :nothing ||
           @capture(pattern, _quote_macrocall) ||
           @capture(pattern, Symbol(_))
        # constant
        quote
            $value == $pattern
        end
    elseif (pattern isa Expr && pattern.head == :$)
        # interpolated value
        quote
            $value == $(esc(pattern.args[1]))
        end
    elseif @capture(pattern, subpattern_Symbol)
        # variable
        # if the pattern doesn't match, we don't want to set the variable
        # so for now just set a temp variable
        our_sym = Symbol("variable_$pattern")
        if pattern in bound
            # already bound, check that this value matches
            quote
                $our_sym == $value
            end
        else
            # bind
            push!(bound, pattern)
            quote
                $our_sym = $value;
                true
            end
        end
    elseif @capture(pattern, subpattern1_ || subpattern2_) || (@capture(pattern, f_(subpattern1_, subpattern2_)) && f == :|)
        # disjunction
        # need to only bind variables which exist in both branches
        bound1 = copy(bound)
        bound2 = copy(bound)
        body1 = handle_destruct(value, subpattern1, bound1, asserts)
        body2 = handle_destruct(value, subpattern2, bound2, asserts)
        union!(bound, intersect(bound1, bound2))
        quote
            $body1 || $body2
        end
    elseif @capture(pattern, subpattern1_ && subpattern2_) || (@capture(pattern, f_(subpattern1_, subpattern2_)) && f == :&)
        # conjunction
        body1 = handle_destruct(value, subpattern1, bound, asserts)
        body2 = handle_destruct(value, subpattern2, bound, asserts)
        quote
            $body1 && $body2
        end
    elseif @capture(pattern, _where)
        # guard
        @assert length(pattern.args) == 2
        subpattern = pattern.args[1]
        guard = pattern.args[2]
        quote
            $(handle_destruct(value, subpattern, bound, asserts)) &&
            let $(bound...)
                # bind variables locally so they can be used in the guard
                $(@splice variable in bound quote
                    $(esc(variable)) = $(Symbol("variable_$variable"))
                end)
                $(esc(guard))
            end
        end
    elseif @capture(pattern, T_(subpatterns__))
        # struct
        len = length(subpatterns)
        named_fields = [pat.args[1] for pat in subpatterns if (pat isa Expr) && pat.head == :kw]
        named_count = length(named_fields)
        @assert named_count == length(unique(named_fields)) "Pattern $pattern has duplicate named arguments ($(named_fields))."
        @assert named_count == 0 || named_count == len "Pattern $pattern mixes named and positional arguments."
        if named_count == 0
            # Pattern uses positional arguments to refer to fields e.g. Foo(0, true)
            expected_fieldcount = gensym("$(T)_expected_fieldcount")
            actual_fieldcount = gensym("$(T)_actual_fieldcount")
            push!(asserts, quote
                if typeof($(esc(T))) <: Function
                    throw(LoadError("Attempted to match on a function", @__LINE__, AssertionError("Incorrect match usage")))
                end
                if !(isstructtype(typeof($(esc(T)))) || issabstracttype(typeof($(esc(T)))))
                    throw(LoadError("Attempted to match on a pattern that is not a struct", @__LINE__, AssertionError("Incorrect match usage")))
                end
                # This assertion is necessary:
                # If $expected_fieldcount < $actual_fieldcount, to catch missing fields.
                # If $expected_fieldcount > $actual_fieldcount, to avoid a BoundsError.
                $expected_fieldcount = evaluated_fieldcount($(esc(T)))
                $actual_fieldcount = $(esc(len))
                if $expected_fieldcount != $actual_fieldcount
                    error("Pattern field count is $($actual_fieldcount) expected $($expected_fieldcount)")
                end

            end)
        else
            # Pattern uses named arguments to refer to fields e.g. Foo(x=0, y=true)
            # Could assert that the expected field names are a subset of the actual field names.
            # Can omit the assertion because if the field doesn't exist getfield() will fail with "type $T has no field $field".
        end
        quote
            # I would prefer typeof($value) == $(esc(T)) but this doesn't convey type information in Julia 0.6
            $value isa $(esc(T)) &&
            $(handle_destruct_fields(value, pattern, subpatterns, length(subpatterns), :getfield, bound, asserts; allow_splat=false))
        end
    elseif @capture(pattern, (subpatterns__,))
        # tuple
        quote
            ($value isa Tuple) &&
            $(handle_destruct_fields(value, pattern, subpatterns, :(length($value)), :getindex, bound, asserts; allow_splat=true))
        end
    elseif @capture(pattern, [subpatterns__])
        # array
        quote
            ($value isa AbstractArray) &&
            $(handle_destruct_fields(value, pattern, subpatterns, :(length($value)), :getindex, bound, asserts; allow_splat=true))
        end
    elseif @capture(pattern, subpattern_::T_)
        # typeassert
        quote
            ($value isa $(esc(T))) &&
            $(handle_destruct(value, subpattern, bound, asserts))
        end
    else
        error("Unrecognized pattern syntax: $pattern")
    end
end

function handle_match_eq(expr)
    if @capture(expr, pattern_ = value_)
        asserts = Expr[]
        bound = Set{Symbol}()
        body = handle_destruct(:value, pattern, bound, asserts)
        quote
            $(asserts...)
            value = $(esc(value))
            $body || throw(MatchFailure(value))
            $(@splice variable in bound quote
                $(esc(variable)) = $(Symbol("variable_$variable"))
            end)
            value
        end
    else
        error("Unrecognized match syntax: $expr")
    end
end

function handle_match_case(value, case, tail, asserts)
    if @capture(case, pattern_ => result_)
        bound = Set{Symbol}()
        body = handle_destruct(:value, pattern, bound, asserts)
        quote
            if $body
                let $(bound...)
                    # export bindings
                    $(@splice variable in bound quote
                      $(esc(variable)) = $(Symbol("variable_$variable"))
                      end)
                    $(esc(result))
                end
            else
                $tail
            end
        end
    else
        error("Unrecognized case syntax: $case")
    end
end

function handle_match_cases(value, match)
    tail = :(throw(MatchFailure(value)))
    if @capture(match, begin cases__ end)
        asserts = Expr[]
        for case in reverse(cases)
            tail = handle_match_case(:value, case, tail, asserts)
        end
        quote
            $(asserts...)
            value = $(esc(value))
            $tail
        end
    else
        error("Unrecognized match syntax: $value $match")
    end
end

"""
    @match pattern = value

If `value` matches `pattern`, bind variables and return `value`. Otherwise, throw `MatchFailure`.
"""
macro match(expr)
    handle_match_eq(expr)
end

"""
    @match value begin
        pattern1 => result1
        pattern2 => result2
        ...
    end

Return `result` for the first matching `pattern`. If there are no matches, throw `MatchFailure`.
"""
macro match(value, cases)
    handle_match_cases(value, cases)
end

"""
Patterns:

  * `_` matches anything
  * `foo` matches anything, binds value to `foo`
  * `Foo(x,y,z)` matches structs of type `Foo` with fields matching `x,y,z`
  * `Foo(y=1)` matches structs of type `Foo` whose `y` field equals `1`
  * `[x,y,z]` matches `AbstractArray`s with 3 entries matching `x,y,z`
  * `(x,y,z)` matches `Tuple`s with 3 entries matching `x,y,z`
  * `[x,y...,z]` matches `AbstractArray`s with at least 2 entries, where `x` matches the first entry, `z` matches the last entry and `y` matches the remaining entries.
  * `(x,y...,z)` matches `Tuple`s with at least 2 entries, where `x` matches the first entry, `z` matches the last entry and `y` matches the remaining entries.
  * `_::T` matches any subtype (`isa`) of T
  * `x || y` matches values which match either `x` or `y` (only variables which exist in both branches will be bound)
  * `x && y` matches values which match both `x` and `y`
  * `x where condition` matches only if `condition` is true (`condition` may use any variables that occur earlier in the pattern eg `(x, y, z where x + y > z)`)
  * Anything else is treated as a constant and tested for equality
  * Expressions can be interpolated in as constants via standard interpolation syntax `\$(x)`

Patterns can be nested arbitrarily.

Repeated variables only match if they are equal (`==`). For example `(x,x)` matches `(1,1)` but not `(1,2)`.
"""
:(@match)

export @match

end


text:
function set_boundary(Qbase,cellxmax,cellymax,vecAx,vecAy,bdcon)
    """bdcon[i][j]
        i : 境界番号(x-,x+ ,y-,y+)
        j=1-6 : "bd1_con":"2",
                "bd1_rho":"1.0",
                "bd1_u":"300.0",
                "bd1_v":"0.0",
                "bd1_p":"1.0",
                "bd1_T":"300.0",
    """
    # bd1 = x-
    if Int(bdcon[1][1]) == 0
        for j in 1:cellymax
            for k in 1:4
                Qbase[1,j,k] = bdcon[1][k+1]
            end
        end
    elseif Int(bdcon[1][1]) == 1
        for j in 1:cellymax
            for k in 1:4
                Qbase[1,j,k] = Qbase[2,j,k]
            end
        end
    elseif Int(bdcon[1][1]) == 2
        for j in 1:cellymax
            xvec = vecAx[2,j,1]
            yvec = vecAx[2,j,2]
            u = Qbase[2,j,2]
            v = Qbase[2,j,3]

            Qbase[1,j,1] = Qbase[2,j,1]
            Qbase[1,j,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[1,j,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
            Qbase[1,j,4] = Qbase[2,j,4]
        end
    elseif Int(bdcon[1][1]) == 3
        for j in 1:cellymax
            u = Qbase[2,j,2]
            v = Qbase[2,j,3]

            Qbase[1,j,1] = Qbase[2,j,1]
            Qbase[1,j,2] = -u
            Qbase[1,j,3] = -v
            Qbase[1,j,4] = Qbase[2,j,4]
        end
    else
        println("\n check boundary condition ! \n")
        throw(UndefVarError(:x))
    end

    # bd2 = x+
    if Int(bdcon[2][1]) == 0
        for j in 1:cellymax
            for k in 1:4
                Qbase[cellxmax,j,k] = bdcon[2][k+1]
            end
        end
    elseif Int(bdcon[2][1]) == 1
        for j in 1:cellymax
            for k in 1:4
                Qbase[cellxmax,j,k] = Qbase[cellxmax-1,j,k]
            end
        end
    elseif Int(bdcon[2][1]) == 2
        for j in 1:cellymax
            xvec = vecAx[cellymax,j,1]
            yvec = vecAx[cellymax,j,2]
            u = Qbase[cellxmax-1,j,2]
            v = Qbase[cellxmax-1,j,3]

            Qbase[cellxmax,j,1] = Qbase[cellxmax-1,j,1]
            Qbase[cellxmax,j,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[cellxmax,j,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
            Qbase[cellxmax,j,4] = Qbase[cellxmax-1,j,4]
        end
    elseif Int(bdcon[2][1]) == 3
        for j in 1:cellymax
            u = Qbase[cellxmax-1,j,2]
            v = Qbase[cellxmax-1,j,3]

            Qbase[cellxmax,j,1] = Qbase[cellxmax-1,j,1]
            Qbase[cellxmax,j,2] = -u
            Qbase[cellxmax,j,3] = -v
            Qbase[cellxmax,j,4] = Qbase[cellxmax-1,j,4]
        end
    else
        println("\n check boundary condition ! \n")
        throw(UndefVarError(:x))
    end

    # bd3 = y-
    if Int(bdcon[3][1]) == 0
        for i in 1:cellxmax
            for k in 1:4
                Qbase[i,1,k] = bdcon[3][k+1]
            end
        end
    elseif Int(bdcon[3][1]) == 1
        for i in 1:cellxmax
            for k in 1:4
                Qbase[i,1,k] = Qbase[i,2,k]
            end
        end
    elseif Int(bdcon[3][1]) == 2
        for i in 1:cellxmax
            xvec = vecAy[i,2,1]
            yvec = vecAy[i,2,2]
            u = Qbase[i,2,2]
            v = Qbase[i,2,3]

            Qbase[i,1,1] = Qbase[i,2,1]
            Qbase[i,1,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[i,1,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
            Qbase[i,1,4] = Qbase[i,2,4]
        end
    elseif Int(bdcon[3][1]) == 3
        for i in 1:cellxmax
            u = Qbase[i,2,2]
            v = Qbase[i,2,3]

            Qbase[i,1,1] = Qbase[i,2,1]
            Qbase[i,1,2] = -u
            Qbase[i,1,3] = -vecAy
            Qbase[i,1,4] = Qbase[i,2,4]
        end
    else
        println("\n check boundary condition ! \n")
        throw(UndefVarError(:x))
    end

    # bd4 = y+
    if Int(bdcon[4][1]) == 0
        for i in 1:cellxmax
            for k in 1:4
                Qbase[i,cellymax,k] = bdcon[4][k+1]
            end
        end
    elseif Int(bdcon[4][1]) == 1
        for i in 1:cellxmax
            for k in 1:4
                Qbase[i,cellymax,k] = Qbase[i,cellxmax-1,k]
            end
        end
    elseif Int(bdcon[4][1]) == 2
        for i in 1:cellxmax
            xvec = vecAy[i,cellymax,1]
            yvec = vecAy[i,cellymax,2]
            u = Qbase[i,cellymax-1,2]
            v = Qbase[i,cellymax-1,3]

            Qbase[i,cellymax,1] = Qbase[i,cellymax-1,1]
            Qbase[i,cellymax,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[i,cellymax,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
            Qbase[i,cellymax,4] = Qbase[i,cellymax-1,4]
        end
    elseif Int(bdcon[4][1]) == 3
        for i in 1:cellxmax
            u = Qbase[i,cellymax-1,2]
            v = Qbase[i,cellymax-1,3]

            Qbase[i,cellymax,1] = Qbase[i,cellymax-1,1]
            Qbase[i,cellymax,2] = -u
            Qbase[i,cellymax,3] = -v
            Qbase[i,cellymax,4] = Qbase[i,cellymax-1,4]
        end
    else
        println("\n check boundary condition ! \n")
        throw(UndefVarError(:x))
    end

    return Qbase
end


text:
import Base: eachindex

eachindex(X::AbstractMatrix, dim::Int) = 1:size(X, dim)
function eachindex(X::AbstractBlockMatrix, dim::Int)
    if dim == 1
        return BlockVector(eachindex.(X.blocks[:, 1], dim))
    elseif dim == 2
        return BlockVector(eachindex.(X.blocks[1, :], dim))
    else
        throw(ErrorException("Bad dimension."))
    end
end
eachindex(X::BS, dim::Int) = eachindex(X.data, dim)
eachindex(X::BlockTri, dim::Int) = eachindex(X.data, dim)
eachindex(X::LazyPDMat, dim::Int) = eachindex(X.Σ, dim)


text:
# TO DO: Comment parse_chstr.jl, even I can't figure out what I've done here anymore
# function parse_charr(chan_in::Array{String,1}; d='.'::Char, fdsn=false::Bool)
function parse_charr(chan_in::Array{String, 1}, d::Char, fdsn::Bool)
    N = length(chan_in)
    chan_data = Array{String,2}(undef, 0, 5)

    # Initial pass to parse to string array
    for i = 1:N
      chan_line = [strip(String(j)) for j in split(chan_in[i], d, keepempty=true, limit=5)]
      L = length(chan_line)
      if L < 5
        resize!(chan_line, 5)
        chan_line[L+1:5] .= ""
      end
      chan_data = vcat(chan_data, reshape(chan_line, 1, 5))
    end

    return fdsn == true ? minreq(String.(chan_data)) : map(String, chan_data)
end

function parse_chstr(chan_in::String, d::Char, fdsn::Bool, SL::Bool)
# function parse_chstr(chan_in::String; d=','::Char, fdsn=false::Bool, SL=false::Bool)
  chan_out = Array{String, 2}(undef, 0, 5)
  if safe_isfile(chan_in)
    return parse_chstr(join([strip(j, ['\r','\n']) for j in
      filter(i -> !startswith(i, ['#','*']), open(readlines, chan_in))],','), ',', false, false)
  else
    chan_data = [strip(String(j)) for j in split(chan_in, d)]
    for j = 1:length(chan_data)

      # Build array
      tmp_data = map(String, split(chan_data[j], '.'))
      L = length(tmp_data)
      if L < 5
        append!(tmp_data, Array{String,1}(undef, 5-L))
        tmp_data[L+1:5] .= ""
      end
      chan_out = cat(chan_out, reshape(tmp_data, 1, 5), dims=1)
    end
  end
  if fdsn == true
    minreq!(chan_out)
  end
  N = SL ? 4 : 5
  return chan_out[:,1:N]
end

function parse_sl(CC::Array{String,2})
  L = size(CC,1)
  S = CC[:,2].*(" "^L).*CC[:,1]
    P = [(i = isempty(i) ? "??" : i) for i in CC[:,3]] .*
        [(i = isempty(i) ? "???" : i) for i in CC[:,4]] .* fill(".", L) .*
        [(i = isempty(i) ? "D" : i) for i in CC[:,5]]
  return S,P
end

# FDSNWS
# http://service.iris.edu/fdsnws/dataselect/1/query?net=IU&sta=ANMO&loc=00&cha=BHZ&start=2010-02-27T06:30:00.000&end=2010-02-27T10:30:00.000
# IRISWS
# http://service.iris.edu/irisws/timeseries/1/query?net=IU&sta=ANMO&loc=00&cha=BHZ&starttime=2005-01-01T00:00:00&endtime=2005-01-02T00:00:00
# IRISWS
function build_stream_query(C::Array{String,1}, d0::String, d1::String; estr=""::String)
  net_str = isempty(C[1]) ? estr : "net="*C[1]
  sta_str = isempty(C[2]) ? estr : "&sta="*C[2]
  loc_str = isempty(C[3]) ? estr : "&loc="*C[3]
  cha_str = isempty(C[4]) ? estr : "&cha="*C[4]
  return net_str * sta_str * loc_str * cha_str * string("&start=", d0, "&end=", d1)
end

# ============================================================================
# Purpose: Create the most compact set of requests, one per row
# """
#     minreq!(S::Array{String,2})
#
# Reduce `S` to the most compact possible set of SeedLink request query strings
# that completely cover its string requests.
# """
function minreq!(S::Array{String,2})
  d = ','
  (M,N) = size(S)
  K = Array{Int64, 1}(undef, N)
  T = Array{String, 2}(undef, M, N)
  for n = 1:N
    for m = 1:M
      T[m,n] = join(S[m,1:N.!=n],d)
    end
    K[n] = length(unique(T[:,n]))
  end
  (L,J) = findmin(K)
  if L != M
    V = T[:,J]
    U = unique(V)
    Q = Array{String,2}(undef, L, N)
    for i = 1:L
      j = findall(V.==U[i])
      Q[i,1:N.!=J] = split(V[j[1]],d)
      Q[i,J] = join(S[j,J],d)
    end
    S = Q
  end
  S = map(String, S)
  return nothing
end
minreq(S::Array{String,2}) = (T = deepcopy(S); minreq!(T); return T)


text:
"""
    rouge_n(references::Array{T}, candidate::AbstractString, n; avg::Bool, lang::Language) where T<: AbstractString

Compute n-gram recall between `candidate` and the `references` summaries.

See [Rouge: A package for automatic evaluation of summaries](http://www.aclweb.org/anthology/W04-1013)

See also: [`rouge_l_sentence`](@ref), [`rouge_l_summary`](@ref)
"""
function rouge_n(references, candidate, n; avg = true, lang = Languages.English())
    ng_candidate = ngramize(lang, candidate, n)
    ng_refs = [ngramize(lang, ref, n) for ref in references]

    rouge_recall = Array{Float64,1}()
    for ref in ng_refs
        push!(rouge_recall, rouge_match_score(keys(ref), ng_candidate) / sum(values(ref)) )
    end

    avg == true && return jackknife_avg(rouge_recall)
    return rouge_recall
end

function rouge_match_score(ref, candidate::Dict)
    matches = 0
    for p in keys(candidate)
        p ∉ ref && continue
        matches += candidate[p]
    end
    return matches
end

"""
    rouge_l_sentence(references, candidate, β, average)

Calculate the ROUGE-L score between `references` and `candidate` at sentence level.

See [Rouge: A package for automatic evaluation of summaries](http://www.aclweb.org/anthology/W04-1013)

See also: [`rouge_n`](@ref), [`rouge_l_summary`](@ref)
"""
function rouge_l_sentence(references, candidate, β=8, average = true)
    ngram_cand = tokenize(Languages.English(), candidate)
    rouge_l_list = []

    for ref in references
        ngram_ref = tokenize(Languages.English(), ref)
        r_lcs = weighted_lcs(ngram_ref, ngram_cand, true, false, sqrt) / length(ngram_ref)
        p_lcs = weighted_lcs(ngram_ref, ngram_cand, true, false, sqrt) / length(ngram_cand)
        score = fmeasure_lcs(r_lcs, p_lcs, β)
        push!(rouge_l_list, score)
    end

    if average == true
        rouge_l_list = jackknife_avg(rouge_l_list)
    end
    return rouge_l_list
end

"""
    rouge_l_summary(references, candidate, β, average)

Calculate the ROUGE-L score between `references` and `candidate` at summary level.

See [Rouge: A package for automatic evaluation of summaries](http://www.aclweb.org/anthology/W04-1013)

See also: [`rouge_l_sentence()`](@ref), [`rouge_l_summary`](@ref)
"""
function rouge_l_summary(references, candidate, β, averaging=true)
    rouge_l_list = []
    cand_sent_list = split_sentences(candidate)

    for ref in references
        ref_sent_list = split_sentences(ref)
        sum_value = 0

        for ref_sent in ref_sent_list
            l_ = []
            arg1 = tokenize(Languages.English(), ref)
            for cand_sent in cand_sent_list
                arg2 = tokenize(Languages.English(), cand_sent)
                d = tokenize(Languages.English(), weighted_lcs(arg1, arg2, false, true, sqrt))
                append!(l_,d)
            end
            sum_value += length(unique(l_))
        end

        r_lcs = sum_value / length(tokenize(Languages.English(), ref))
        p_lcs = sum_value / length(tokenize(Languages.English(), candidate))
        score = fmeasure_lcs(r_lcs, p_lcs, β)
        push!(rouge_l_list,score)
    end

    averaging == true && return jackknife_avg(rouge_l_list)
    return rouge_l_list
end


text:
import RestrictedBoltzmannMachines as RBMs
using Test: @test, @testset
using Statistics: mean

@testset "initialization Binary" begin
    data = rand(2, 3, 10^6) .≤ 3/4
    rbm = RBMs.RBM(RBMs.Binary(2,3), RBMs.Binary(0), randn(2,3,0))
    RBMs.initialize!(rbm, data)
    @test RBMs.mean_from_inputs(rbm.visible) ≈ reshape(mean(data; dims=3), size(rbm.visible)) rtol=0.01
end

@testset "initialization Spin" begin
    data = Int8.(sign.(rand(2, 3, 10^6) .- 1/4))
    rbm = RBMs.RBM(RBMs.Spin(2,3), RBMs.Binary(0), randn(2,3,0))
    RBMs.initialize!(rbm, data)
    @test RBMs.mean_from_inputs(rbm.visible) ≈ reshape(mean(data; dims=3), size(rbm.visible)) rtol=0.01
end

@testset "initialization Potts" begin
    rbm = RBMs.RBM(RBMs.Potts(3,2,3), RBMs.Binary(0), randn(3,2,3,0))
    data = RBMs.onehot_encode(RBMs.onehot_decode(rand(3,2,3,10^6) .- rand(3,2,3,1)), 1:3)
    @assert all(sum(data; dims=1) .== 1)
    RBMs.initialize!(rbm, data)
    @test RBMs.mean_from_inputs(rbm.visible) ≈ reshape(mean(data; dims=4), size(rbm.visible)) rtol=0.01
end