text:
using DimensionalData, Test, Unitful, Combinatorics, Dates
using DimensionalData.LookupArrays, DimensionalData.Dimensions
using .LookupArrays: between, at, near, contains

a = [1 2  3  4
     5 6  7  8
     9 10 11 12]

dims_ = X(Sampled(10:10:20; sampling=Intervals())),
        Y(Sampled(5:7; sampling=Intervals()))
A = DimArray([1 2 3; 4 5 6], dims_)

@testset "selector primitives" begin

    @testset "Explicit Intervals" begin
        fwdindex = 11.0:30.0
        revindex = 30.0:-1.0:11.0

        args = Intervals(Start()), NoMetadata()
        fwdmatrix = vcat(fwdindex', (fwdindex .+ 1)')
        revmatrix = vcat(revindex', (revindex .+ 1)')
        startfwd = Sampled(fwdindex, ForwardOrdered(), Explicit(fwdmatrix), args...)
        startrev = Sampled(revindex, ReverseOrdered(), Explicit(revmatrix), args...)

        args = Intervals(Center()), NoMetadata()
        fwdmatrix = vcat((fwdindex .- 0.5)', (fwdindex .+ 0.5)')
        revmatrix = vcat((revindex .- 0.5)', (revindex .+ 0.5)')
        centerfwd = Sampled(fwdindex, ForwardOrdered(), Explicit(fwdmatrix), args...)
        centerrev = Sampled(revindex, ReverseOrdered(), Explicit(revmatrix), args...)

        args = Intervals(End()), NoMetadata()
        fwdmatrix = vcat((fwdindex .- 1)', fwdindex')
        revmatrix = vcat((revindex .- 1)', revindex')
        endfwd = Sampled(fwdindex, ForwardOrdered(), Explicit(fwdmatrix), args...)
        endrev = Sampled(revindex, ReverseOrdered(), Explicit(revmatrix), args...)

        @testset "Any at" begin
            @test at(startfwd, At(30)) == 20
            @test at(startrev, At(30)) == 1
            @test at(startfwd, At(29.9; atol=0.2)) == 20
            @test at(startrev, At(29.9; atol=0.2)) == 1
            @test at(startfwd, At(30.1; atol=0.2)) == 20
            @test at(startrev, At(30.1; atol=0.2)) == 1
        end

        @testset "Start between" begin
            @test between(startfwd, Between(0, 11.9)) === 1:0
            @test between(startfwd, Between(0, 12)) === 1:1
            @test between(startfwd, Between(11, 14)) === 1:3
            @test between(startfwd, Between(30, 50)) === 20:20
            @test between(startfwd, Between(31, 50)) === 21:20
            @test between(startfwd, Between(11.1, 13.9)) === 2:2
            @test between(startrev, Between(0, 11.9)) === 21:20
            @test between(startrev, Between(0, 12)) === 20:20
            @test between(startrev, Between(11, 14)) === 18:20
            @test between(startrev, Between(11.1, 13.9)) === 19:19
            @test between(startrev, Between(30, 50)) === 1:1
            @test between(startrev, Between(30.1, 50)) === 1:0
            @test between(startfwd, Between(0, 40)) === 1:20
            @test between(startrev, Between(0, 40)) === 1:20
            # Input order doesn't matter
            @test between(startfwd, Between(14, 11)) === 1:3
        end

        @testset "Center between" begin
            @test between(centerfwd, Between(0, 11.4)) === 1:0
            @test between(centerfwd, Between(0, 11.5)) === 1:1
            @test between(centerfwd, Between(10.5, 14.6)) === 1:4
            @test between(centerfwd, Between(10.6, 14.4)) === 2:3
            @test between(centerfwd, Between(29.5, 50.0)) === 20:20
            @test between(centerfwd, Between(29.6, 50.0)) === 21:20
            @test between(centerrev, Between(10.5, 14.6)) === 17:20
            @test between(centerrev, Between(10.6, 14.4)) === 18:19
            @test between(centerrev, Between(0, 11.4)) === 21:20
            @test between(centerrev, Between(0, 11.5)) === 20:20
            @test between(centerrev, Between(29.5, 50.0)) === 1:1
            @test between(centerrev, Between(29.6, 50.0)) === 1:0
            @test between(centerfwd, Between(0, 40)) === 1:20
            @test between(centerrev, Between(0, 40)) === 1:20
            # Input order doesn't matter
            @test between(centerfwd, Between(15, 10)) === 1:4
        end

        @testset "End between" begin
            @test between(endfwd, Between(0.0, 10.9)) === 1:0
            @test between(endfwd, Between(0, 11)) === 1:1
            @test between(endfwd, Between(10.1, 14.9)) === 2:4
            @test between(endfwd, Between(10, 15)) === 1:5
            @test between(endfwd, Between(29.0, 30.0)) === 20:20
            @test between(endfwd, Between(29.1, 50.0)) === 21:20
            @test between(endrev, Between(0.0, 10.9)) === 21:20
            @test between(endrev, Between(0.0, 11.0)) === 20:20
            @test between(endrev, Between(10.1, 14.9)) === 17:19
            @test between(endrev, Between(10, 15)) === 16:20
            @test between(endrev, Between(29.0, 50.0)) === 1:1
            @test between(endrev, Between(29.1, 50.0)) === 1:0
            @test between(endfwd, Between(0, 40)) === 1:20
            @test between(endrev, Between(0, 40)) === 1:20
            # Input order doesn't matter
            @test between(endfwd, Between(15, 10)) === 1:5
        end

        @testset "Start contains" begin
            @test_throws BoundsError contains(startfwd, Contains(10.9))
            @test_throws BoundsError contains(startfwd, Contains(31))
            @test_throws BoundsError contains(startrev, Contains(31))
            @test_throws BoundsError contains(startrev, Contains(10.9))
            @test contains(startfwd, Contains(11)) == 1
            @test contains(startfwd, Contains(11.9)) == 1
            @test contains(startfwd, Contains(12.0)) == 2
            @test contains(startfwd, Contains(30.0)) == 20
            @test contains(startfwd, Contains(30.9)) == 20
            @test contains(startfwd, Contains(29.9)) == 19
            @test contains(startrev, Contains(11.9)) == 20
            @test contains(startrev, Contains(12.0)) == 19
            @test contains(startrev, Contains(30.0)) == 1
            @test contains(startrev, Contains(30.9)) == 1
            @test contains(startrev, Contains(29.0)) == 2
        end

        @testset "Center contains" begin
            @test_throws BoundsError contains(centerfwd, Contains(10.4))
            @test_throws BoundsError contains(centerfwd, Contains(30.5))
            @test_throws BoundsError contains(centerrev, Contains(10.4))
            @test_throws BoundsError contains(centerrev, Contains(30.5))
            @test contains(centerfwd, Contains(10.5)) == 1
            @test contains(centerfwd, Contains(30.4)) == 20
            @test contains(centerfwd, Contains(29.5)) == 20
            @test contains(centerfwd, Contains(29.4)) == 19
            @test contains(centerrev, Contains(10.5)) == 20
            @test contains(centerrev, Contains(30.4)) == 1
            @test contains(centerrev, Contains(29.5)) == 1
            @test contains(centerrev, Contains(29.4)) == 2
        end

        @testset "End contains" begin
            @test_throws BoundsError contains(endfwd, Contains(10))
            @test_throws BoundsError contains(endfwd, Contains(30.1))
            @test_throws BoundsError contains(endrev, Contains(10))
            @test_throws BoundsError contains(endrev, Contains(30.1))
            @test contains(endfwd, Contains(10.1)) == 1
            @test contains(endfwd, Contains(11.0)) == 1
            @test contains(endfwd, Contains(11.1)) == 2
            @test contains(endfwd, Contains(29.0)) == 19
            @test contains(endfwd, Contains(29.1)) == 20
            @test contains(endfwd, Contains(30.0)) == 20
            @test contains(endrev, Contains(10.1)) == 20
            @test contains(endrev, Contains(11.0)) == 20
            @test contains(endrev, Contains(11.1)) == 19
            @test contains(endrev, Contains(29.0)) == 2
            @test contains(endrev, Contains(29.1)) == 1
            @test contains(endrev, Contains(30.0)) == 1
            @test contains(endrev, Contains(11.0)) == 20
            @test contains(endrev, Contains(11.1)) == 19
        end

        @testset "Center near" begin
            @test near(centerfwd, Near(10.4)) == 1
            @test near(centerfwd, Near(30.5)) == 20
            @test near(centerrev, Near(10.4)) == 20
            @test near(centerrev, Near(30.5)) == 1
            @test near(centerfwd, Near(10.5)) == 1
            @test near(centerfwd, Near(30.4)) == 20
            @test near(centerfwd, Near(29.5)) == 20
            @test near(centerfwd, Near(29.4)) == 19
            @test near(centerrev, Near(10.5)) == 20
            @test near(centerrev, Near(30.4)) == 1
            @test near(centerrev, Near(29.5)) == 1
            @test near(centerrev, Near(29.4)) == 2
        end

    end

    @testset "Regular Intervals with range" begin
        args = Intervals(Start()), NoMetadata()
        startfwd = Sampled(11.0:30.0,      ForwardOrdered(), Regular(1), args...)
        startrev = Sampled(30.0:-1.0:11.0, ReverseOrdered(), Regular(-1), args...)

        args = Intervals(Center()), NoMetadata()
        centerfwd = Sampled(11.0:30.0,      ForwardOrdered(), Regular(1), args...)
        centerrev = Sampled(30.0:-1.0:11.0, ReverseOrdered(), Regular(-1), args...)

        args = Intervals(End()), NoMetadata()
        endfwd = Sampled(11.0:30.0,      ForwardOrdered(), Regular(1), args...)
        endrev = Sampled(30.0:-1.0:11.0, ReverseOrdered(), Regular(-1), args...)

        @testset "Any at" begin
            @test at(startfwd, At(30)) == 20
            @test at(startrev, At(30)) == 1
            @test at(startfwd, At(29.9; atol=0.2)) == 20
            @test at(startrev, At(29.9; atol=0.2)) == 1
            @test at(startfwd, At(30.1; atol=0.2)) == 20
            @test at(startrev, At(30.1; atol=0.2)) == 1
        end

        @testset "Start between" begin
            @test between(startfwd, Between(0, 11.9)) === 1:0
            @test between(startfwd, Between(0, 12)) === 1:1
            @test between(startfwd, Between(11, 14)) === 1:3
            @test between(startfwd, Between(30, 50)) === 20:20
            @test between(startfwd, Between(31, 50)) === 21:20
            @test between(startfwd, Between(11.1, 13.9)) === 2:2
            @test between(startrev, Between(0, 11.9)) === 21:20
            @test between(startrev, Between(0, 12)) === 20:20
            @test between(startrev, Between(11, 14)) === 18:20
            @test between(startrev, Between(11.1, 13.9)) === 19:19
            @test between(startrev, Between(30, 50)) === 1:1
            @test between(startrev, Between(30.1, 50)) === 1:0
            @test between(startfwd, Between(0, 40)) === 1:20
            @test between(startrev, Between(0, 40)) === 1:20
            # Input order doesn't matter
            @test between(startfwd, Between(14, 11)) === 1:3
        end

        @testset "Center between" begin
            @test between(centerfwd, Between(0, 11.4)) === 1:0
            @test between(centerfwd, Between(0, 11.5)) === 1:1
            @test between(centerfwd, Between(10.5, 14.6)) === 1:4
            @test between(centerfwd, Between(10.6, 14.4)) === 2:3
            @test between(centerfwd, Between(29.5, 50.0)) === 20:20
            @test between(centerfwd, Between(29.6, 50.0)) === 21:20
            @test between(centerrev, Between(10.5, 14.6)) === 17:20
            @test between(centerrev, Between(10.6, 14.4)) === 18:19
            @test between(centerrev, Between(0, 11.4)) === 21:20
            @test between(centerrev, Between(0, 11.5)) === 20:20
            @test between(centerrev, Between(29.5, 50.0)) === 1:1
            @test between(centerrev, Between(29.6, 50.0)) === 1:0
            @test between(centerfwd, Between(0, 40)) === 1:20
            @test between(centerrev, Between(0, 40)) === 1:20
            # Input order doesn't matter
            @test between(centerfwd, Between(15, 10)) === 1:4
        end

        @testset "End between" begin
            @test between(endfwd, Between(0.0, 10.9)) === 1:0
            @test between(endfwd, Between(0, 11)) === 1:1
            @test between(endfwd, Between(10.1, 14.9)) === 2:4
            @test between(endfwd, Between(10, 15)) === 1:5
            @test between(endfwd, Between(29.0, 30.0)) === 20:20
            @test between(endfwd, Between(29.1, 50.0)) === 21:20
            @test between(endrev, Between(0.0, 10.9)) === 21:20
            @test between(endrev, Between(0.0, 11.0)) === 20:20
            @test between(endrev, Between(10.1, 14.9)) === 17:19
            @test between(endrev, Between(10, 15)) === 16:20
            @test between(endrev, Between(29.0, 50.0)) === 1:1
            @test between(endrev, Between(29.1, 50.0)) === 1:0
            @test between(endfwd, Between(0, 40)) === 1:20
            @test between(endrev, Between(0, 40)) === 1:20
            # Input order doesn't matter
            @test between(endfwd, Between(15, 10)) === 1:5
        end

        @testset "Start contains" begin
            @test_throws BoundsError contains(startfwd, Contains(10.9))
            @test_throws BoundsError contains(startrev, Contains(31))
            @test contains(startfwd, Contains(11)) == 1
            @test contains(startfwd, Contains(11.9)) == 1
            @test contains(startfwd, Contains(12.0)) == 2
            @test contains(startfwd, Contains(30.0)) == 20
            @test contains(startfwd, Contains(29.9)) == 19
            @test contains(startrev, Contains(11.9)) == 20
            @test contains(startrev, Contains(12.0)) == 19
            @test contains(startrev, Contains(30.9)) == 1
            @test contains(startrev, Contains(30.0)) == 1
            @test contains(startrev, Contains(29.0)) == 2
        end

        @testset "Center contains" begin
            @test_throws BoundsError contains(centerfwd, Contains(10.4))
            @test_throws BoundsError contains(centerfwd, Contains(30.5))
            @test_throws BoundsError contains(centerrev, Contains(10.4))
            @test_throws BoundsError contains(centerrev, Contains(30.5))
            @test contains(centerfwd, Contains(10.5)) == 1
            @test contains(centerfwd, Contains(30.4)) == 20
            @test contains(centerfwd, Contains(29.5)) == 20
            @test contains(centerfwd, Contains(29.4)) == 19
            @test contains(centerrev, Contains(10.5)) == 20
            @test contains(centerrev, Contains(30.4)) == 1
            @test contains(centerrev, Contains(29.5)) == 1
            @test contains(centerrev, Contains(29.4)) == 2
        end

        @testset "End contains" begin
            @test_throws BoundsError contains(endfwd, Contains(10))
            @test_throws BoundsError contains(endfwd, Contains(30.1))
            @test_throws BoundsError contains(endrev, Contains(10))
            @test_throws BoundsError contains(endrev, Contains(30.1))
            @test contains(endfwd, Contains(10.1)) == 1
            @test contains(endfwd, Contains(11.0)) == 1
            @test contains(endfwd, Contains(11.1)) == 2
            @test contains(endfwd, Contains(29.0)) == 19
            @test contains(endfwd, Contains(29.1)) == 20
            @test contains(endfwd, Contains(30.0)) == 20
            @test contains(endrev, Contains(10.1)) == 20
            @test contains(endrev, Contains(11.0)) == 20
            @test contains(endrev, Contains(11.1)) == 19
            @test contains(endrev, Contains(29.0)) == 2
            @test contains(endrev, Contains(29.1)) == 1
            @test contains(endrev, Contains(30.0)) == 1
            @test contains(endrev, Contains(11.0)) == 20
            @test contains(endrev, Contains(11.1)) == 19
        end

        @testset "Start near" begin
            @test bounds(startfwd) == bounds(startrev)
            @test near(startfwd, Near(-100)) == 1
            @test near(startfwd, Near(11.9)) == 1
            @test near(startfwd, Near(12.0)) == 2
            @test near(startfwd, Near(30.0)) == 20
            @test near(startfwd, Near(29.9)) == 19
            @test near(startrev, Near(11.9)) == 20
            @test near(startrev, Near(12.0)) == 19
            @test near(startrev, Near(29.0)) == 2
            @test near(startrev, Near(30.0)) == 1
            @test near(startfwd, Near(100)) == 20
        end

        @testset "Center near" begin
            @test near(centerfwd, Near(10.4)) == 1
            @test near(centerfwd, Near(30.5)) == 20
            @test near(centerrev, Near(10.4)) == 20
            @test near(centerrev, Near(30.5)) == 1
            @test near(centerfwd, Near(10.5)) == 1
            @test near(centerfwd, Near(30.4)) == 20
            @test near(centerfwd, Near(29.5)) == 20
            @test near(centerfwd, Near(29.4)) == 19
            @test near(centerrev, Near(10.5)) == 20
            @test near(centerrev, Near(30.4)) == 1
            @test near(centerrev, Near(29.5)) == 1
            @test near(centerrev, Near(29.4)) == 2
        end

        @testset "End near" begin
            @test near(endfwd, Near(10)) == 1
            @test near(endfwd, Near(11.0)) == 1
            @test near(endfwd, Near(11.1)) == 2
            @test near(endfwd, Near(29.0)) == 19
            @test near(endfwd, Near(29.1)) == 20
            @test near(endfwd, Near(30.0)) == 20
            @test near(endfwd, Near(31.1)) == 20
            @test near(endrev, Near(10)) == 20
            @test near(endrev, Near(11.0)) == 20
            @test near(endrev, Near(11.1)) == 19
            @test near(endrev, Near(29.0)) == 2
            @test near(endrev, Near(29.1)) == 1
            @test near(endrev, Near(30.0)) == 1
            @test near(endrev, Near(31.1)) == 1
            @test near(endrev, Near(11.0)) == 20
            @test near(endrev, Near(11.1)) == 19
            @test near(endrev, Near(29.0)) == 2
            @test near(endrev, Near(29.1)) == 1
            @test near(endfwd, Near(-100)) == 1
            @test near(endfwd, Near(100)) == 20
        end

    end

    @testset "Regular Intervals with array" begin
        startfwd = Sampled([1, 3, 4, 5], ForwardOrdered(), Regular(1), Intervals(Start()), NoMetadata())
        startrev = Sampled([5, 4, 3, 1], ReverseOrdered(), Regular(-1), Intervals(Start()), NoMetadata())
        @test_throws BoundsError contains(startfwd, Contains(0.9))
        @test contains(startfwd, Contains(1.0)) == 1
        @test contains(startfwd, Contains(1.9)) == 1
        @test_throws ArgumentError contains(startfwd, Contains(2))
        @test_throws ArgumentError contains(startfwd, Contains(2.9))
        @test contains(startfwd, Contains(3)) == 2
        @test contains(startfwd, Contains(5.9)) == 4
        @test_throws BoundsError contains(startfwd, Contains(6))
        @test_throws BoundsError contains(startrev, Contains(0.9))
        @test contains(startrev, Contains(1.0)) == 4
        @test contains(startrev, Contains(1.9)) == 4
        @test_throws ArgumentError contains(startrev, Contains(2))
        @test_throws ArgumentError contains(startrev, Contains(2.9))
        @test contains(startrev, Contains(3)) == 3
        @test contains(startrev, Contains(5.9)) == 1
        @test_throws BoundsError contains(startrev, Contains(6))
    end

    @testset "Irregular Intervals with array" begin
        # Order: index, array, relation (array order is irrelevent here, it's just for plotting)
        # Varnames: locusindexorderrelation
        args = Irregular(1.0, 121.0), Intervals(Start()), NoMetadata()
        startfwd = Sampled((1:10).^2,    ForwardOrdered(), args...)
        startrev = Sampled((10:-1:1).^2, ReverseOrdered(), args...)

        args = Irregular(0.5, 111.5), Intervals(Center()), NoMetadata()
        centerfwd = Sampled((1.0:10.0).^2,      ForwardOrdered(), args...)
        centerrev = Sampled((10.0:-1.0:1.0).^2, ReverseOrdered(), args...)

        args = Irregular(0.0, 100.0), Intervals(End()), NoMetadata()
        endfwd = Sampled((1.0:10.0).^2,      ForwardOrdered(), args...)
        endrev = Sampled((10.0:-1.0:1.0).^2, ReverseOrdered(), args...)

        @testset "Any at" begin
            @test at(startfwd, At(25)) == 5
            @test at(startrev, At(25)) == 6
            @test at(startfwd, At(100)) == 10
            @test at(startrev, At(100)) == 1
        end

        @testset "Start between" begin
            @test between(startfwd, Between(9, 36)) === 3:5
            @test between(startfwd, Between(9.1, 35.0)) === 4:4
            @test between(startrev, Between(9, 36)) === 6:8
            @test between(startrev, Between(9.1, 35.9)) === 7:7
            # Input order doesn't matter
            @test between(startfwd, Between(36, 9)) === 3:5
            # Handle searchorted overflow
            @test between(startfwd, Between(0, 0)) === 1:0
            @test between(startfwd, Between(130, 130)) === 11:10
            @test between(startrev, Between(0, 0)) === 11:10
            @test between(startrev, Between(130, 130)) === 1:0
            @test between(startfwd, Between(-100, 9)) === 1:2
            @test between(startfwd, Between(80, 150)) === 9:10
            @test between(startrev, Between(-100, 9)) === 9:10
            @test between(startrev, Between(80, 150)) === 1:2
            @test between(startfwd, Between(-200, 200)) === 1:10
            @test between(startrev, Between(-200, 200)) === 1:10
        end

        @testset "Center between" begin
            @test between(centerfwd, Between(6.5, 30.5)) === 3:5
            @test between(centerfwd, Between(6.6, 30.4)) === 4:4
            @test between(centerrev, Between(6.5, 30.5)) === 6:8
            @test between(centerrev, Between(6.6, 30.4)) === 7:7
            # Input order doesn't matter
            @test between(centerfwd, Between(30.5, 6.5)) === 3:5
            # Handle searchorted overflow
            @test between(centerfwd, Between(0, 0)) === 1:0
            @test between(centerfwd, Between(-100, 9)) === 1:2
            @test between(centerfwd, Between(70, 150)) === 9:10
            @test between(centerfwd, Between(130, 130)) === 11:10
            @test between(centerrev, Between(0, 0)) === 11:10
            @test between(centerrev, Between(-100, 9)) === 9:10
            @test between(centerrev, Between(70, 150)) === 1:2
            @test between(centerrev, Between(130, 130)) === 1:0
            @test between(centerfwd, Between(-200, 200)) === 1:10
            @test between(centerrev, Between(-200, 200)) === 1:10
        end

        @testset "End between" begin
            @test between(endfwd, Between(4, 25)) === 3:5
            @test between(endrev, Between(4, 25)) === 6:8
            @test between(endfwd, Between(4.1, 24.9)) === 4:4
            @test between(endrev, Between(4.1, 24.9)) === 7:7
            # Input order doesn't matter
            @test between(endfwd, Between(25, 4)) === 3:5
            # Handle searchorted overflow
            @test between(endfwd, Between(0, 0)) === 1:0
            @test between(endfwd, Between(130, 130)) === 11:10
            @test between(endrev, Between(0, 0)) === 11:10
            @test between(endrev, Between(130, 130)) === 1:0
            @test between(endfwd, Between(-100, 4)) === 1:2
            @test between(endfwd, Between(64, 150)) === 9:10
            @test between(endrev, Between(-100, 4)) === 9:10
            @test between(endrev, Between(64, 150)) === 1:2
            @test between(endfwd, Between(-200, 200)) === 1:10
            @test between(endrev, Between(-200, 200)) === 1:10
        end

        @testset "Start contains" begin
            @test_throws BoundsError contains(startfwd, Contains(0.9))
            @test_throws BoundsError contains(startfwd, Contains(121.1))
            @test contains(startfwd, Contains(1)) == 1
            @test contains(startfwd, Contains(3.9)) == 1
            @test contains(startfwd, Contains(4.0)) == 2
            @test contains(startfwd, Contains(100.0)) == 10
            @test contains(startfwd, Contains(99.9)) == 9
            @test_throws BoundsError contains(startrev, Contains(0.9))
            @test_throws BoundsError contains(startrev, Contains(121.1))
            @test contains(startrev, Contains(3.9)) == 10
            @test contains(startrev, Contains(4.0)) == 9
            @test contains(startrev, Contains(120.9)) == 1
            @test contains(startrev, Contains(100.0)) == 1
            @test contains(startrev, Contains(99.0)) == 2
        end

        @testset "Center contains" begin
            @test_throws BoundsError contains(centerfwd, Contains(0.4))
            @test_throws BoundsError contains(centerfwd, Contains(111.5))
            @test contains(centerfwd, Contains(0.5)) == 1
            @test contains(centerfwd, Contains(111.4)) == 10
            @test contains(centerfwd, Contains(90.5)) == 10
            @test contains(centerfwd, Contains(90.4)) == 9
            @test_throws BoundsError contains(centerrev, Contains(0.4))
            @test_throws BoundsError contains(centerrev, Contains(111.5))
            @test contains(centerrev, Contains(72.5)) == 2
            @test contains(centerrev, Contains(72.4)) == 3
            @test contains(centerrev, Contains(0.5)) == 10
            @test contains(centerrev, Contains(111.4)) == 1
            @test contains(centerrev, Contains(90.5)) == 1
            @test contains(centerrev, Contains(90.4)) == 2
        end

        @testset "End contains" begin
            @test_throws BoundsError contains(endfwd, Contains(-0.1))
            @test_throws BoundsError contains(endfwd, Contains(100.1))
            @test_throws BoundsError contains(endrev, Contains(-0.1))
            @test_throws BoundsError contains(endrev, Contains(100.1))
            @test contains(endfwd, Contains(0.1)) == 1
            @test contains(endfwd, Contains(1.0)) == 1
            @test contains(endfwd, Contains(1.1)) == 2
            @test contains(endfwd, Contains(81.0)) == 9
            @test contains(endfwd, Contains(81.1)) == 10
            @test contains(endfwd, Contains(100.0)) == 10
            @test contains(endrev, Contains(0.1)) == 10
            @test contains(endrev, Contains(1.0)) == 10
            @test contains(endrev, Contains(1.1)) == 9
            @test contains(endrev, Contains(81.0)) == 2
            @test contains(endrev, Contains(81.1)) == 1
            @test contains(endrev, Contains(100.0)) == 1
        end

    end

    @testset "Points lookup" begin
        fwd = Sampled(5.0:30.0;      order=ForwardOrdered(), sampling=Points())
        rev = Sampled(30.0:-1.0:5.0; order=ReverseOrdered(), sampling=Points())

        @testset "between" begin
            @test between(fwd, Between(0, 4.9)) === 1:0
            @test between(fwd, Between(0, 5)) === 1:1
            @test between(fwd, Between(30, 31)) === 27:26
            @test between(fwd, Between(10, 14.9)) === 6:10
            @test between(fwd, Between(10, 15)) === 6:11
            @test between(rev, Between(10, 14.9)) === 17:21
            @test between(rev, Between(10, 15)) === 16:21
            @test between(rev, Between(0, 4.9)) === 27:26
            @test between(rev, Between(0, 5)) === 26:26
            @test between(rev, Between(30, 31)) === 1:0
            # Input order doesn't matter
            @test between(fwd, Between(15, 10)) === 6:11
        end

        @testset "at" begin
            @test at(fwd, At(30)) == 26
            @test at(rev, At(30)) == 1
        end

        @testset "near" begin
            @test near(fwd, Near(50))   == 26
            @test near(fwd, Near(0))    == 1
            @test near(fwd, Near(29.4)) == 25
            @test near(fwd, Near(29.5)) == 26
            @test near(rev, Near(29.4)) == 2
            @test near(rev, Near(30.1)) == 1
            @test_throws ArgumentError near(Sampled((5.0:30.0); order=Unordered(), sampling=Points()), Near(30.1))
        end

        @testset "contains" begin
            @test_throws ArgumentError contains(fwd, Contains(50))
        end

    end

end

@testset "Selectors on Sampled Points" begin
    da = DimArray(a, (Y(Sampled(10:10:30)), Ti(Sampled((1:4)u"s"))))

    @test At(10.0) == At(10.0, nothing, nothing)
    @test At(10.0; atol=0.0, rtol=Base.rtoldefault(Float64)) ==
          At(10.0, 0.0, Base.rtoldefault(Float64))
    Near([10, 20])

    @test Between(10, 20) == Between((10, 20))

    @testset "selectors with dim wrappers" begin
        @test @inferred da[Y(At([10, 30])), Ti(At([1u"s", 4u"s"]))] == [1 4; 9 12]
        @test_throws ArgumentError da[Y(At([9, 30])), Ti(At([1u"s", 4u"s"]))]
        @test @inferred view(da, Y(At(20)), Ti(At((3:4)u"s"))) == [7, 8]
        @test @inferred view(da, Y(Near(17)), Ti(Near([1.5u"s", 3.1u"s"]))) == [6, 7]
        @test @inferred view(da, Y(Between(9, 21)), Ti(At((3:4)u"s"))) == [3 4; 7 8]
    end

    @testset "selectors without dim wrappers" begin
        @test @inferred da[At(20:10:30), At(1u"s")] == [5, 9]
        @test @inferred view(da, Between(9, 31), Near((3:4)u"s")) == [3 4; 7 8; 11 12]
        @test @inferred view(da, Near(22), At([3.0u"s", 4.0u"s"])) == [7, 8]
        @test @inferred view(da, At(20), At((2:3)u"s")) == [6, 7]
        @test @inferred view(da, Near(13), Near([1.3u"s", 3.3u"s"])) == [1, 3]
        # Near works with a tuple input
        @test @inferred view(da, Near([13]), Near([1.3u"s", 3.3u"s"])) == [1 3]
        @test @inferred view(da, Between(11, 20), At((2:3)u"s")) == [6 7]
        # Between also accepts a tuple input
        @test @inferred view(da, Between((11, 20)), Between((2u"s", 3u"s"))) == [6 7]
    end

    @testset "mixed selectors and standard" begin
        selectors = [
            (Between(9, 31), Near((3:4)u"s")),
            (Near(22), At([3.0u"s", 4.0u"s"])),
            (At(20), At((2:3)u"s")),
            (Near(13), Near([1.3u"s", 3.3u"s"])),
            (Near([13]), Near([1.3u"s", 3.3u"s"])),
            (Between(11, 20), At((2:3)u"s"))
        ]
        positions =  [
            (1:3, [3, 4]),
            (2, [3, 4]),
            (2, [2, 3]),
            (1, [1, 3]),
            ([1], [1, 3]),
            (2:2, [2, 3])
        ]
        for (selector, pos) in zip(selectors, positions)
            pairs = collect(zip(selector, pos))
            cases = [(i, j) for i in pairs[1], j in pairs[2]]
            for (case1, case2) in combinations(cases, 2)
                @test da[case1...] == da[case2...]
                @test view(da, case1...) == view(da, case2...)
                dac1, dac2 = copy(da), copy(da)
                sample = da[case1...]
                replacement  = sample isa Integer ? 100 : rand(Int, size(sample))
                # Test return value
                @test setindex!(dac1, replacement, case1...) == setindex!(dac2, replacement, case2...)
                # Test mutation
                @test dac1 == dac2
            end
        end
    end

    @testset "single-arity standard index" begin
        indices = [
            1:3,
            [1, 2, 4],
            4:-2:1,
        ]
        for idx in indices
            from2d = da[idx]
            @test from2d == parent(da)[idx]
            @test !(from2d isa AbstractDimArray)
            from1d = da[Y(At(10))][idx]
            @test from1d == parent(da)[1, :][idx]
            @test from1d isa AbstractDimArray
        end
    end

    @testset "single-arity views" begin
        indices = [
            3,
            1:3,
            [1, 2, 4],
            4:-2:1,
        ]
        for idx in indices
            idx = 3
            from2d = view(da, idx)
            @test from2d == view(parent(da), idx)
            @test from2d isa SubArray
            from1d = view(da[Y(At(10))], idx)
            @test from1d == view(parent(da)[1, :], idx)
            @test from1d isa SubArray
        end
    end

    @testset "single-arity setindex!" begin
        indices = [
            3,
            1:3,
            [1, 2, 4],
            4:-2:1,
        ]
        for idx in indices
            idx = indices[1]
            # 2D case
            da2d = copy(da)
            a2d = copy(parent(da2d))
            replacement = zero(a2d[idx])
            @test setindex!(da2d, replacement, idx) == setindex!(a2d, replacement, idx)
            @test da2d == a2d
            # 1D array
            da1d = da[Y(At(10))]
            a1d = copy(parent(da1d))
            @test setindex!(da1d, replacement, idx) == setindex!(a1d, replacement, idx)
            @test da1d == a1d
        end
    end

    @testset "more Unitful dims" begin
        dimz = Ti(1.0u"s":1.0u"s":3.0u"s"), Y(1u"km":1u"km":4u"km")
        db = DimArray(a, dimz)
        @test db[Y(Between(2u"km", 3.9u"km")), Ti(At(3.0u"s"))] == [10, 11]
    end

    @testset "selectors work in reverse orders" begin
        a = [1 2  3  4
             5 6  7  8
             9 10 11 12]

        @testset "forward index" begin
            da_fff = DimArray(a, (Y(Sampled(10:10:30; order=ForwardOrdered())),
                                  Ti(Sampled((1:4)u"s"; order=ForwardOrdered()))))
            @test da_fff[Y(At(20)), Ti(At((3.0:4.0)u"s"))] == [7, 8]
            @test da_fff[Y(At([20, 30])), Ti(At((3.0:4.0)u"s"))] == [7 8; 11 12]
            @test da_fff[Y(Near(22)), Ti(Near([3.3u"s", 4.3u"s"]))] == [7, 8]
            @test da_fff[Y(Near([22, 42])), Ti(Near([3.3u"s", 4.3u"s"]))] == [7 8; 11 12]
            @test da_fff[Y(Between(20, 30)), Ti(Between(3.0u"s", 4.0u"s"))] == [7 8; 11 12]
        end

        @testset "reverse index" begin
            da_rff = DimArray(a, (Y(Sampled(30:-10:10; order=ReverseOrdered())),
                                  Ti(Sampled((4:-1:1)u"s"; order=ReverseOrdered()))))
            @test da_rff[Y(At(20)), Ti(At((3.0:4.0)u"s"))] == [6, 5]
            @test da_rff[Y(At([20, 30])), Ti(At((3.0:4.0)u"s"))] == [6 5; 2 1]
            @test da_rff[Y(Near(22)), Ti(Near([3.3u"s", 4.3u"s"]))] == [6, 5]
            @test da_rff[Y(Near([22, 42])), Ti(Near([3.3u"s", 4.3u"s"]))] == [6 5; 2 1]
            # Between hasn't reverse the index order
            @test da_rff[Y(Between(20, 30)), Ti(Between(3.0u"s", 4.0u"s"))] == [1 2; 5 6]
        end

    end

    @testset "setindex! with selectors" begin
        c = deepcopy(a)
        dc = DimArray(c, (Y(10:10:30), Ti((1:4)u"s")))
        dc[Near(11), At(3u"s")] = 100
        @test c[1, 3] == 100
        @inferred setindex!(dc, [200, 201, 202], Ti(Near(2.2u"s")), Y(Between(10, 30)))
        @test c[1:3, 2] == [200, 201, 202]
    end

    @testset "Where " begin
        @test val(Where(identity)) == identity
        dimz = Ti((1:1:3)u"s"), Y(10:10:40)
        da = DimArray(a, dimz)
        wda = da[Y(Where(x -> x >= 30)), Ti(Where(x -> x in([2u"s", 3u"s"])))]
        @test parent(wda) == [7 8; 11 12]
        @test index(wda) == ([2u"s", 3u"s"], [30, 40])
    end

end

@testset "Selectors on Sampled Intervals" begin
    da = DimArray(a, (Y(Sampled(10:10:30; sampling=Intervals())),
                      Ti(Sampled((1:4)u"s"; sampling=Intervals()))))

    @testset "with dim wrappers" begin
        @test @inferred da[Y(At([10, 30])), Ti(At([1u"s", 4u"s"]))] == [1 4; 9 12]
        @test_throws ArgumentError da[Y(At([9, 30])), Ti(At([1u"s", 4u"s"]))]
        @test @inferred view(da, Y(At(20)), Ti(At((3:4)u"s"))) == [7, 8]
        @test @inferred view(da, Y(Contains(17)), Ti(Contains([1.9u"s", 3.1u"s"]))) == [5, 7]
        @test @inferred view(da, Y(Between(4, 26)), Ti(At((3:4)u"s"))) == [3 4; 7 8]
    end

    @testset "without dim wrappers" begin
        @test @inferred da[At(20:10:30), At(1u"s")] == [5, 9]
        @test @inferred view(da, Between(4, 36), Near((3:4)u"s")) == [3 4; 7 8; 11 12]
        @test @inferred view(da, Near(22), At([3.0u"s", 4.0u"s"])) == [7, 8]
        @test @inferred view(da, At(20), At((2:3)u"s")) == [6, 7]
        @test @inferred view(da, Near(13), Near([1.3u"s", 3.3u"s"])) == [1, 3]
        @test @inferred view(da, Near([13]), Near([1.3u"s", 3.3u"s"])) == [1 3]
        @test @inferred view(da, Between(11, 26), At((2:3)u"s")) == [6 7]
        # Between also accepts a tuple input
        @test @inferred view(da, Between((11, 26)), Between((2u"s", 4u"s"))) == [6 7]
    end

    @testset "out of bounds" begin
        @test size(view(da, Between(0, 4), At((2:3)u"s"))) == (0, 2)
        @test view(da, Between(0, 4), At((2:3)u"s")) isa DimArray{Int64,2}
        @test size(view(da, Between(40, 45), At((2:3)u"s"))) == (0, 2)
        @test view(da, Between(40, 45), At((2:3)u"s")) isa DimArray{Int64,2}
    end

    @testset "with DateTime index" begin
        @testset "Start locus" begin
            timedim = Ti(Sampled(DateTime(2001):Month(1):DateTime(2001, 12);
                span=Regular(Month(1)), sampling=Intervals(Start())
            ))
            da = DimArray(1:12, timedim)
            @test @inferred da[Ti(At(DateTime(2001, 3)))] == 3
            @test @inferred da[Near(DateTime(2001, 4, 7))] == 4
            @test @inferred da[Between(DateTime(2001, 4, 7), DateTime(2001, 8, 30))] == [5, 6, 7]
        end
        @testset "End locus" begin
            timedim = Ti(Sampled(DateTime(2001):Month(1):DateTime(2001, 12);
                span=Regular(Month(1)), sampling=Intervals(End()))
            )
            da = DimArray(1:12, timedim)
            @test @inferred da[Ti(At(DateTime(2001, 3)))] == 3
            @test @inferred da[Near(DateTime(2001, 4, 7))] == 5
            @test @inferred da[Between(DateTime(2001, 4, 7), DateTime(2001, 8, 30))] == [6, 7, 8]
        end
    end

    # @testset "with Val index" begin
    #     valdimz = Ti(Sampled(Val(2.0:2.0:6.0))),
    #               Y(Sampled(Val(10.0:10.0:40)))
    #     da = DimArray(a, valdimz)
    #     @test @inferred da[Ti=Val(4.0), Y=Val(40.0)] == 8
    #     @test @inferred da[At(2.0), At(20.0)] == 2
    #     @test @inferred da[Near(Val{3.2}()), At(Val{20.0}())] == 6
    #     @test @inferred da[Near(3.2), At(20.0)] == 6
    # end

end

@testset "Selectors on Sampled Explicit Intervals" begin
    span_y = Explicit(vcat((5.0:10.0:25.0)', (15.0:10.0:35.0)'))
    span_ti = Explicit(vcat(((0.5:1.0:3.5)u"s")', ((1.5:1.0:4.5)u"s")'))
    da = DimArray(a, (Y(Sampled(10:10:30; span=span_y, sampling=Intervals())),
                      Ti(Sampled((1:4)u"s"; span=span_ti, sampling=Intervals()))))

    @testset "with dim wrappers" begin
        @test @inferred da[Y(At([10, 30])), Ti(At([1u"s", 4u"s"]))] == [1 4; 9 12]
        @test_throws ArgumentError da[Y(At([9, 30])), Ti(At([1u"s", 4u"s"]))]
        @test @inferred view(da, Y(At(20)), Ti(At((3:4)u"s"))) == [7, 8]
        @test @inferred view(da, Y(Contains(17)), Ti(Contains([1.4u"s", 3.1u"s"]))) == [5, 7]
    end

    @testset "without dim wrappers" begin
        @test @inferred da[At(20:10:30), At(1u"s")] == [5, 9]
        @test @inferred view(da, Between(4, 36), Contains((3:4)u"s")) == [3 4; 7 8; 11 12]
        @test @inferred view(da, Contains(22), At([3.0u"s", 4.0u"s"])) == [7, 8]
        @test @inferred view(da, At(20), At((2:3)u"s")) == [6, 7]
        @test @inferred view(da, Contains(13), Contains([1.3u"s", 3.3u"s"])) == [1, 3]
        @test @inferred view(da, Contains([13]), Contains([1.3u"s", 3.3u"s"])) == [1 3]
        @test @inferred view(da, Between(11, 26), At((2:3)u"s")) == [6 7]
        # Between also accepts a tuple input
        @test @inferred view(da, Between((11, 26)), Between((1.4u"s", 4u"s"))) == [6 7]
    end

    @testset "with DateTime index" begin
        span_ti = Explicit(vcat(
            reshape((DateTime(2001, 1):Month(1):DateTime(2001, 12)), 1, 12),
            reshape((DateTime(2001, 2):Month(1):DateTime(2002, 1)), 1, 12)
        ))
        timedim = Ti(Sampled(DateTime(2001, 1, 15):Month(1):DateTime(2001, 12, 15);
            span=span_ti, sampling=Intervals(Center())
        ))
        da = DimArray(1:12, timedim)
        @test @inferred da[Ti(At(DateTime(2001, 3, 15)))] == 3
        @test @inferred da[Contains(DateTime(2001, 4, 7))] == 4
        @test @inferred da[Between(DateTime(2001, 4, 7), DateTime(2001, 8, 30))] == [5, 6, 7]
    end

end

@testset "Selectors on Categorical" begin
    a = [1 2  3  4
         5 6  7  8
         9 10 11 12]

    dimz = Ti(Categorical([:one, :two, :three]; order=ForwardOrdered())),
           Y(Categorical([:a, :b, :c, :d]; order=ForwardOrdered()))
    da = DimArray(a, dimz)
    @test @inferred da[Ti(At([:one, :two])), Y(Contains(:b))] == [2, 6]
    @test @inferred da[At(:two), Between(:b, :d)] == [6, 7, 8]
    @test @inferred da[At(:two), At(:b)] == 6
    # Near and contains are just At
    @test @inferred da[Contains([:one, :three]), Near([:b, :c, :d])] == [2 3 4; 10 11 12]

    dimz = Ti(Categorical([:one, :two, :three]; order=Unordered())),
           Y(Categorical([:a, :b, :c, :d]; order=Unordered()))
    da = DimArray(a, dimz)
    @test_throws ArgumentError da[At(:two), Between(:b, :d)] == [6, 7, 8]

    unordered_dimz = Ti(Categorical([:one, :two, :three]; order=Unordered())),
                     Y(Categorical([:a, :b, :c, :d]; order=Unordered()))
    unordered_da = DimArray(a, unordered_dimz)
    unordered_da[Near(:two), Near(:d)]

    a = [1 2  3  4
         5 6  7  8
         9 10 11 12]
end

@testset "Selectors on NoLookup" begin
    dimz = Ti(), Y()
    da = DimArray(a, dimz)
    @test @inferred da[Ti(At([1, 2])), Y(Contains(2))] == [2, 6]
    @test @inferred da[Near(2), Between(2, 4)] == [6, 7, 8]
    @test @inferred da[Contains([1, 3]), Near([2, 3, 4])] == [2 3 4; 10 11 12]
end

@testset "Selectors on TranformedIndex" begin
    using CoordinateTransformations

    m = LinearMap([0.5 0.0; 0.0 0.5])

    dimz = Dim{:trans1}(Transformed(m, X())),
           Dim{:trans2}(Transformed(m, Y())), Z()
    @test DimensionalData.dimsmatch(dimz[1], X())
    @test DimensionalData.dimsmatch(dimz[2], Y())

    @testset "permutedims works on lookup dimensions" begin
        # @test @inferred sortdims((Y(), Z(), X()), dimz) == (X(), Y(), Z())
    end

    da = DimArray(reshape(a, 3, 4, 1), dimz)

    @testset "Indexing with array dims indexes the array as usual" begin
        da2 = da[1:3, 1:1, 1:1]
        @test @inferred da2[Dim{:trans1}(3), Dim{:trans2}(1), Z(1)] == 9
        # Using selectors works the same as indexing with lookup
        # dims - it applies the transform function.
        # It's not clear this should be allowed or makes sense,
        # but it works anyway because the permutation is correct either way.
        @test @inferred da[Dim{:trans1}(At(6)), Dim{:trans2}(At(2)), Z(1)] == 9
    end

    @testset "Indexing with lookup dims uses the transformation" begin
        @test @inferred da[X(Near(6.1)), Y(Near(8.5)), Z(1)] == 12
        @test @inferred da[X(At(4.0)), Y(At(2.0)), Z(1)] == 5
        @test_throws ArgumentError da[trans1=At(4.0)]
        @test_throws InexactError da[X(At(6.1)), Y(At(8)), Z(1)]
        # Indexing directly with lookup dims also just works, but maybe shouldn't?
        @test @inferred da[X(2), Y(2), Z(1)] == 6
    end

end

@testset "selectindices" begin
    @test selectindices(A[X(1)], Contains(7)) == (3,)
    @test selectindices(dims_, ()) == ()
end

@testset "errors" begin
    @test_throws ArgumentError DimensionalData.selectindices(X(Sampled(1:4, sampling=Points())), Contains(1))
end

@testset "hasselection" begin
    @test hasselection(A, X(At(20)))
    @test hasselection(dims(A, X), X(At(20)))
    @test hasselection(dims(A, X), At(19; atol=2))
    @test hasselection(A, (Y(At(7)),))
    @test hasselection(A, (X(At(10)), Y(At(7))))

    args = Intervals(Start()), NoMetadata()
    startfwd = Ti(Sampled(11.0:30.0,      ForwardOrdered(), Regular(1), args...))
    startrev = Ti(Sampled(30.0:-1.0:11.0, ReverseOrdered(), Regular(-1), args...))

    cases = (startfwd, startrev)
    @test all(map(d -> hasselection(d, At(30.0)), cases))
    @test all(map(d -> !hasselection(d, At(31.0)), cases))
    @test all(map(d -> hasselection(d, Contains(12.8)), cases))
    @test all(map(d -> !hasselection(d, Contains(400.0)), cases))
    @test all(map(d -> hasselection(d, Near(0.0)), cases))
end


text:
using ACEds

using StaticArrays
using ACEds.equivRotations3D: _mrange

import LinearAlgebra: dot

cg=ClebschGordan()

cg(5, -4, 5, 4, 2, 0)

ll = @SVector [1,1,3,2]
mm = @SVector [1,-1,1,-1]
kk = @SVector [1,-1,1,-1]
irotc = ACEds.equivRotations3D.IRot3DCoeffs(Float64)
irotc(ll,mm,kk)

#a=3
#a::Int64

#dicttype(::Val{N}) where {N} =
#   Dict{Tuple{SVector{N,Int}, SVector{N,Int}, SVector{N,Int}}, Float64}

#SMatrix{3,3,Float64}
#

#SMatrix{3,3,Float64,9}
erotc = ACEds.equivRotations3D.ERot3DCoeffs(Float64)
erotc(ll,mm,kk)

"""
Playground: some simple tests
"""

length(_mrange(@SVector([1,2,3,1])))

for (im, mm) in enumerate(_mrange(@SVector([2,2,2])))
   print((im,mm),", \n")
end

struct Squares
   count::Int
end

Base.iterate(S::Squares, state=1) = state > S.count ? nothing : (state*state, state+1)
for i in Squares(7)
   println(i)
end

ll = @SVector([1])
Iterators.Stateful(CartesianIndices(ntuple(i -> -ll[i]:ll[i], length(ll)-1)))

for g in Iterators.Stateful(CartesianIndices(ntuple(i -> -ll[i]:ll[i], length(ll)-1)))
   print(g,"\n")
end

g= filter((x) -> abs(sum(x))<= 1, Tuple.(CartesianIndices(ntuple(i -> -ll[i]:ll[i], length(ll)))))
gi = Iterators.Stateful(g)
length(gi)

for (i,val) in enumerate(gi)
   print(i,val)
end
for i in Tuple(CartesianIndex(1, 2, -3, -1))
  println(i)
end

###################

"""
Given an l-vector `ll` iterate over all combinations of `mm` vectors  of
the same length such that `sum(mm) == 1`
"""
struct MRange2{N, T2}
   ll::SVector{N, Int}
   cartrg::T2
end

Base.length(mr::MRange2) = sum(_->1, _mrange2(mr.ll))

_mrange2(ll) = MRange2(ll, Iterators.Stateful(
               filter((x) -> abs(sum(x))<= 1, Tuple.(CartesianIndices(ntuple(i -> -ll[i]:ll[i], length(ll)))))
                     ))

function Base.iterate(mr::MRange2, args...)
   while true
      if isempty(mr.cartrg)
         return nothing
      end
      mpre = popfirst!(mr.cartrg)
      return SVector(mpre), nothing
   end
end

"""
Compare behaviour of iterate(mr::MRange) with iterate(mr::MRange2)
"""
ll = @SVector([1,2])
mr = ACEds.equivRotations3D._mrange(ll)
for (i,val) in enumerate(mr)
   print(i,val)
end

#mr2 = _mrange2(ll)
mr2 = ACEds.equivRotations3D._mrange1(ll)
for (i,val) in enumerate(mr2)
   print(i,",", val,"\n")
end

collect(mr2)
a = collect(Iterators.product(1:3, mr2))

#for (i,val) in Iterators.product(1:3, Squares2(7))
#   print(i,",",val,"\n")
#end

#collect(Iterators.product(1:2, Squares(7)))

#struct Squares2
#   count::Int
#end

#Base.iterate(S::Squares2, state=1) = state > S.count ? nothing : ( string("blabla ", state*state), state+1)

for i in mr2
   println(i)
end

mr2 = ACEds.equivRotations3D._mrange1(ll)
for (counter, (i,k)) in enumerate(Iterators.product(1:3, mr2))
   print(counter,",",i,",",k,"\n")
end

a=SMatrix{3, 3, ComplexF64, 9}(1/6, 1im/6, 0, -1im/6, 1/6, 0, 0, 0, 0)

function erot_dot(jmu,mu, jm,m)
   E =
   return dot
end
print(dot(a[:,1],a[:,1]))

print("\n")
print("----------------")
print("\n")
print(a)
print("\n")
print(a[2,1])
print("\n")
print(conj(conj(a[:,1])))

#compute_gl(A::ERot3DCoeffs{T}, ll::SVector, ::Val{false})

ext: jl


text:
using CoupledDipoles, Revise, LinearAlgebra, Folds
using Random

#=
This file contain steps to compute improve the computation of
    Light scatterd on the far field approximation:

I(θ, φ) = ∑ₙ∑ₘ conj(βₙ)βₘ exp( im*k₀*n⋅(rₙ - rₘ) ) + ∑ₙ (1 + zₙ)/2

∑ₙ∑ₘ = n=1:N; m=1:N; n≠m

βⱼ = atom expected value of σ⁻ (j ∈ [n,m])
zₙ = atom expected value of σᶻ
k₀ = unit of length of the system (k₀ = 1)
n = versor of the sensor ( n = r_sensor./norm(r_sensor) )
rⱼ = atom position (j ∈ [n,m])

=#

"""
    Following the definition of I(θ, φ) without any optimization
"""
function scattering_version_naive(β,z, n_hat, r; k₀=1)
    N = length(β)
    intensity = ComplexF64(0)

    for n=1:N
        for m=1:N
            if n≠m
                intensity += conj(β[n])*β[m]*exp( im*k₀*dot(n_hat, r[:, n]-r[:, m])   )
            end
        end
    end
    for n=1:N
        intensity += ( 1 + z[n] )/2
    end
    return intensity
end

"""
   - I compute `conj(β[n])*β[m]` once
   - avoid accessing position rₙ inside loop m (r_n = @view(r[:, n]))
   - put @view to access all collumns
"""
function scattering_optimization_1(β, z, n_hat, r; k₀=1)
    N = length(β)

    βₙₘ = zeros(ComplexF64, N, N)
    for n=1:N
        for m=1:N
            if n≠m
                βₙₘ[n,m] = conj(β[n])*β[m]
            end
        end
    end

    intensity = ComplexF64(0)
    for n=1:N
        r_n = @view(r[:, n])
        for m=1:N
            if n≠m
                intensity += βₙₘ[n,m]*exp( im*k₀*dot(n_hat, r_n - @view(r[:,m])  )   )
            end
        end
    end
    for n=1:N
        intensity += (1+z[n])/2
    end
    return intensity
end

"""
   - `zeros(N,N)` --> `Array{ComplexF64}(undef, N,N)`
   - expand the dot product
   - change `exp(im*x)` by `cis(x)`
   - removed the loop for `∑ₙ (1 + zₙ)/2`, and used a `sum` function
"""
function scattering_optimization_2(β, z, n_hat, r; k₀=1)
    N = length(β)

    βₙₘ = Array{ComplexF64}(undef, N,N)
    for n=1:N
        for m=1:N
            if n≠m
                βₙₘ[n,m] = conj(β[n])*β[m]
            end
        end
    end

    intensity = ComplexF64(0)
    r_nm = zeros(3)
    for n=1:N
        r_n = @view(r[:, n])
        for m=1:N
            if n≠m
                r_nm .= r_n - @view(r[:,m])
                dot_n_r = n_hat[1]*r_nm[1] + n_hat[2]*r_nm[2] + n_hat[3]*r_nm[3]
                intensity += βₙₘ[n,m]*cis( k₀*dot_n_r  )
            end
        end
    end

    intensity += sum( (1 .+ z)./2 )

    return intensity
end

"""
   - stored in a matriz the vectors of `rₙ - rₘ`
"""
function scattering_optimization_3(β, z, n_hat, r; k₀=1)
    N = length(β)

    βₙₘ = Array{ComplexF64}(undef, N,N)
    for n=1:N
        for m=1:N
            βₙₘ[n,m] = conj(β[n])*β[m]
        end
    end

    rₙₘ = Array{ComplexF64}(undef, 3, N^2)
    cont = 1
    for n=1:N
        r_n = @view(r[:,n])
        for m=1:N
            rₙₘ[:,cont] = r_n - @view(r[:,m])
            cont += 1
        end
    end

    intensity = ComplexF64(0)
    cont = 1
    for n=1:N
        for m=1:N
            if n≠m
                dot_n_r = n_hat[1]*rₙₘ[1, cont] + n_hat[2]*rₙₘ[2, cont] + n_hat[3]*rₙₘ[3, cont]
                intensity += βₙₘ[n,m]*cis( k₀*dot_n_r  )
            end
            cont += 1
        end
    end

    intensity += sum( (1 .+ z)./2 )

    return intensity
end

"""
    - using complex number relation: 2real(z) = z + conj(z)
        - Then I return: 2real(intensity)
    - I don't need to compute `conj(z)` terms, which are the lower diagonal
        of βₙₘ and rₙₘ.
    - Also, to avoid the `if n≠m`, I begin with upper diagonal  `m = (n+1):N`
    for n=1:N
        for m = (n+1):N  # <--- here is the change
            ...
        end
    end
    - I cannot forget to divide the diagonal by 2:
        `sum( (1 .+ z)./2 )/2` --> `sum( (1 .+ z)./4 )`
    - Isolate the main computation inside the function `aux4`
"""
function scattering_optimization_4(β, z, n_hat, r; k₀=1)
    N = length(β)

    number_configurations = ((N^2)÷2 - N÷2)

    βₙₘ = Array{ComplexF64}(undef, number_configurations)
    cont = 1
    for n=1:N
        for m=(n+1):N
            βₙₘ[cont] = conj(β[n])*β[m]
            cont += 1
        end
    end

    rₙₘ = Array{ComplexF64}(undef, 3, number_configurations)
    cont = 1
    for n=1:N
        r_n = @view(r[:,n])
        for m=(n+1):N
            rₙₘ[1,cont] = r_n[1] - r[1,m]
            rₙₘ[2,cont] = r_n[2] - r[2,m]
            rₙₘ[3,cont] = r_n[3] - r[3,m]
            cont += 1
        end
    end

    intensity = ComplexF64(0)
    cont = 1
    for n=1:N
        for m=(n+1):N
            intensity += aux4(n_hat, rₙₘ, βₙₘ, cont)
            cont += 1
        end
    end

    intensity +=  sum( (1 .+ z)./4 )
    return 2real(intensity)
end
function aux4(n_hat, rₙₘ, βₙₘ, cont; k₀=1)
    dot_n_r = n_hat[1]*rₙₘ[1, cont] + n_hat[2]*rₙₘ[2, cont] + n_hat[3]*rₙₘ[3, cont]
    return βₙₘ[cont]*cis( k₀*dot_n_r  )
end

"""
    - create views of βₙₘ and rₙₘ.
    - do not call for external function to compute the intensity
"""
function scattering_optimization_5(β, z, n_hat, r; k₀=1)
    N = length(β)

    number_configurations = ((N^2)÷2 - N÷2)

    βₙₘ = Array{ComplexF64}(undef, number_configurations)
    cont = 1
    for n=1:N
        for m=(n+1):N
            βₙₘ[cont] = conj(β[n])*β[m]
            cont += 1
        end
    end
    vβₙₘ = view(βₙₘ, :)

    rₙₘ = Array{ComplexF64}(undef, 3, number_configurations)
    cont = 1
    for n=1:N
        r_n = @view(r[:,n])
        for m=(n+1):N
            rₙₘ[1,cont] = r_n[1] - r[1,m]
            rₙₘ[2,cont] = r_n[2] - r[2,m]
            rₙₘ[3,cont] = r_n[3] - r[3,m]
            cont += 1
        end
    end
    vrₙₘ = view(rₙₘ,:, :)

    intensity = ComplexF64(0)
    for cont = 1:number_configurations
        dot_n_r = n_hat[1]*vrₙₘ[1, cont] + n_hat[2]*vrₙₘ[2, cont] + n_hat[3]*vrₙₘ[3, cont]
        intensity += vβₙₘ[cont]*cis( k₀*dot_n_r  )
    end

    intensity +=  sum( (1 .+ z)./4 )
    return 2real(intensity)
end

"""
    - use multi-threading of `Folds.mapreduce`
    - I still need to use another @view when slicing vrₙₘ
"""
function scattering_optimization_6(β, z, n_hat, r; k₀=1)
    N = length(β)

    number_configurations = ((N^2)÷2 - N÷2)

    βₙₘ = Array{ComplexF64}(undef, number_configurations)
    cont = 1
    for n=1:N
        for m=(n+1):N
            βₙₘ[cont] = conj(β[n])*β[m]
            cont += 1
        end
    end
    vβₙₘ = view(βₙₘ, :)

    rₙₘ = Array{ComplexF64}(undef, 3, number_configurations)
    cont = 1
    for n=1:N
        r_n = @view(r[:,n])
        for m=(n+1):N
            rₙₘ[1,cont] = r_n[1] - r[1,m]
            rₙₘ[2,cont] = r_n[2] - r[2,m]
            rₙₘ[3,cont] = r_n[3] - r[3,m]
            cont += 1
        end
    end
    vrₙₘ = view(rₙₘ,:, :)

    intensity = ComplexF64(0)
    intensity = Folds.mapreduce(+, 1:number_configurations) do k
        (  aux_opt6(n_hat, @view(vrₙₘ[:, k]), vβₙₘ[k]) )
    end
    intensity +=  sum( (1 .+ z)./4 )
    return 2real(intensity)
end
function aux_opt6(n_hat, vrₙₘ, vβₘₙ; k₀=1)
    dot_n_r = cis( k₀*(n_hat[1]*vrₙₘ[1]
                    + n_hat[2]*vrₙₘ[2]
                    + n_hat[3]*vrₙₘ[3]))
   return dot_n_r*vβₘₙ
end

"""
    - do not call for external functions inside mapreduce
    - use @inbounds
"""
function scattering_optimization_7(β, z, n_hat, r; k₀=1)
    N = length(β)

    number_configurations = ((N^2)÷2 - N÷2)

    βₙₘ = Array{ComplexF64}(undef, number_configurations)
    cont = 1
    for n=1:N
        for m=(n+1):N
            βₙₘ[cont] = conj(β[n])*β[m]
            cont += 1
        end
    end
    vβₙₘ = view(βₙₘ, :)

    rₙₘ = Array{ComplexF64}(undef, 3, number_configurations)
    cont = 1
    for n=1:N
        r_n = @view(r[:,n])
        for m=(n+1):N
            rₙₘ[1,cont] = r_n[1] - r[1,m]
            rₙₘ[2,cont] = r_n[2] - r[2,m]
            rₙₘ[3,cont] = r_n[3] - r[3,m]
            cont += 1
        end
    end
    vrₙₘ = view(rₙₘ,:, :)

    intensity = ComplexF64(0)
    intensity = Folds.mapreduce(+, 1:number_configurations) do k
        (
            begin
                @inbounds dot_n_r = n_hat[1]*vrₙₘ[1, k] + n_hat[2]*vrₙₘ[2, k] + n_hat[3]*vrₙₘ[3, k];
                @inbounds vβₙₘ[k]*cis( k₀*dot_n_r)
            end
        )
    end
    intensity +=  sum( (1 .+ z)./4 )
    return 2real(intensity)
end

# initialize
N = 400
kL  = 10
s = 1e-6
Δ = 0.0
sensors = get_sensors_ring(;num_pts=90, kR=1.5kL, θ=5π/12)

Random.seed!(1301)
atoms = Atom(Cube(), N, kL)
laser = Laser(Gaussian3D(kL/8), s, Δ)
simulation = NonLinearOptics(MeanField(), atoms, laser)

β = rand(ComplexF64, N)
z = 2conj(β).*β .- 1
r = atoms.r
n = sensors[:,1]./norm(sensors[:,1])

print("naive:");@time  I_naive = scattering_version_naive(β,z, n, r);
print("opt1: "); @time I_opt_1 = scattering_optimization_1(β,z, n, r); @assert I_naive ≈ I_opt_1
print("opt2: "); @time I_opt_2 = scattering_optimization_2(β,z, n, r); @assert I_naive ≈ I_opt_2
print("opt3: "); @time I_opt_3 = scattering_optimization_3(β,z, n, r); @assert I_naive ≈ I_opt_3
print("opt4: "); @time I_opt_4 = scattering_optimization_4(β,z, n, r); @assert I_naive ≈ I_opt_4
print("opt5: "); @time I_opt_5 = scattering_optimization_5(β,z, n, r); @assert I_naive ≈ I_opt_5
print("opt6: "); @time I_opt_6 = scattering_optimization_6(β,z, n, r); @assert I_naive ≈ I_opt_6
print("opt7: "); @time I_opt_7 = scattering_optimization_7(β,z, n, r); @assert I_naive ≈ I_opt_7

# benchmarks
N = 6000

Random.seed!(1301)
atoms = Atom(Cube(), N, kL)
laser = Laser(Gaussian3D(kL/8), s, Δ)
simulation = NonLinearOptics(MeanField(), atoms, laser)

β = rand(ComplexF64, N)
z = 2conj(β).*β .- 1
r = atoms.r

print("naive:");@time  I_naive = scattering_version_naive(β,z, n, r);
print("opt1: "); @time I_opt_1 = scattering_optimization_1(β,z, n, r); @assert I_naive ≈ I_opt_1
print("opt2: "); @time I_opt_2 = scattering_optimization_2(β,z, n, r); @assert I_naive ≈ I_opt_2
print("opt3: "); @time I_opt_3 = scattering_optimization_3(β,z, n, r); @assert I_naive ≈ I_opt_3
print("opt4: "); @time I_opt_4 = scattering_optimization_4(β,z, n, r); @assert I_naive ≈ I_opt_4
print("opt5: "); @time I_opt_5 = scattering_optimization_5(β,z, n, r); @assert I_naive ≈ I_opt_5
print("opt6: "); @time I_opt_6 = scattering_optimization_6(β,z, n, r); @assert I_naive ≈ I_opt_6
print("opt7: "); @time I_opt_7 = scattering_optimization_7(β,z, n, r); @assert I_naive ≈ I_opt_7

"""
    Code for production is ideal up to 90 sensors,
        because the function takes advantage of CPU frequency burst
"""
function scattering_current_best(β, z, r, sensors; k₀=1)
    N = length(β)

    number_configurations = ((N^2)÷2 - N÷2)

    βₙₘ = Array{ComplexF64}(undef, number_configurations)
    cont = 1
    for n=1:N
        for m=(n+1):N
            βₙₘ[cont] = conj(β[n])*β[m]
            cont += 1
        end
    end
    vβₙₘ = view(βₙₘ, :)

    rₙₘ = Array{ComplexF64}(undef, 3, number_configurations)
    cont = 1
    for n=1:N
        r_n = @view(r[:,n])
        for m=(n+1):N
            rₙₘ[1,cont] = r_n[1] - r[1,m]
            rₙₘ[2,cont] = r_n[2] - r[2,m]
            rₙₘ[3,cont] = r_n[3] - r[3,m]
            cont += 1
        end
    end
    vrₙₘ = view(rₙₘ,:, :)

    n_sensors = size(sensors, 2)
    intensities = Float64[]

    if n_sensors==1
        n_hat = sensors./norm(sensors)
        intensity = _oneSensor_MeanField_scattering(n_hat, vβₙₘ, vrₙₘ, number_configurations)
        intensity +=  sum( (1 .+ z)./4 )

        push!(intensities, 2real(intensity)   )
    else
        const_sum = sum( (1 .+ z)./4 )
        for oneSensor in eachcol(sensors)
            n_hat =  oneSensor./norm(oneSensor)

            intensity = _manySensors_MeanField_scattering(n_hat, vβₙₘ, vrₙₘ, number_configurations)
            intensity +=  const_sum

            push!(intensities, 2real(intensity)   )
        end
    end
    return intensities
end
function _oneSensor_MeanField_scattering(n_hat, vβₙₘ, vrₙₘ, number_configurations; k₀=1)
    intensity = ComplexF64(0)
    for cont = 1:number_configurations
        dot_n_r = n_hat[1]*vrₙₘ[1, cont] + n_hat[2]*vrₙₘ[2, cont] + n_hat[3]*vrₙₘ[3, cont]
        intensity += vβₙₘ[cont]*cis( k₀*dot_n_r  )
    end
    return intensity
end

function _manySensors_MeanField_scattering(n_hat, vβₙₘ, vrₙₘ, number_configurations; k₀=1)
    intensity = ComplexF64(0)
    intensity = Folds.mapreduce(+, 1:number_configurations) do k
        (
            begin
                @inbounds vβₙₘ[k]*cis( k₀*(n_hat[1]*vrₙₘ[1, k] + n_hat[2]*vrₙₘ[2, k] + n_hat[3]*vrₙₘ[3, k]))
            end
        )
    end
    return intensity
end
@profview scattering_current_best(β, z, atoms.r, sensors)

for i=1:10
    @time scattering_current_best(β, z, atoms.r, sensors)
end


text:
struct DummyFilm <: GRAYTR.Film
    resX::Int64
    resY::Int64
end
DummyFilm() = DummyFilm(10,10)

@testset "Orthographic camera" begin
    @testset "Orthographic function" begin
        T = GRAYTR.orthographic(0.0, 1.0)
        @test T.M ≈ Matrix{Float64}(I, 4, 4)
        @test T.MInv ≈ Matrix{Float64}(I, 4, 4)
        T = GRAYTR.orthographic(1.0, 3.0)
        @test diag(T.M) ≈ [1.0, 1.0, 0.5, 1.0]
        @test T.M[:,4] ≈ [0.0, 0.0, -0.5, 1.0]
    end
    @testset "Creations" begin
        T = Transformation()
        window = [-1.0, 1.0, -1.0, 1.0]
        C = OrthographicCamera(T, window, 0.0, 0.0, DummyFilm())
        sample = GRAYTR.CameraSample(5.0, 5.0, 0.5, 0.5)
        weight, R = GRAYTR.generate_ray(C, sample)
        @test R.origin ≈ Point3(0.0, 0.0, 0.0)
        @test R.direction ≈ Vector3(0.0, 0.0, 1.0)
        @test R.tmin == 0.0
        @test R.tmax == Inf
        @test R.depth == 1
    end
    @testset "Transformations" begin
        T = Transformation()
        window = [-1.0, 1.0, -1.0, 1.0]
        C = OrthographicCamera(T, window, 0.0, 0.0, DummyFilm())
        @test C.camera_to_world.M ≈ Matrix{Float64}(I, 4, 4)
        @test C.raster_to_camera(Point3(0)) ≈ Point3(-1, -1, 0)
        @test C.screen_to_camera(Point3(0)) ≈ Point3(0)
    end
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "OpenCL_Headers"
version = v"2020.03.13"

# Collection of sources required to build this package
sources = [
    ArchiveSource("https://github.com/KhronosGroup/OpenCL-Headers/archive/v2020.03.13.tar.gz",
                  "664bbe587e5a0a00aac267f645b7c413586e7bc56dca9ff3b00037050d06f476"),
]

# Bash recipe for building across all platforms
script = raw"""
cd OpenCL-Headers-*
install_license LICENSE

mkdir ${prefix}/include
mv CL ${prefix}/include/
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
products = [
    FileProduct("include/CL/cl.h", :cl_h)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
### Diagonal ops

# we use Base.Diagonal(f) for diagonal operators so very little specific code is
# actually needed here.

simulate(rng::AbstractRNG, D::DiagOp{F}) where {F<:Field} = sqrt(D) * white_noise(rng, F)
global_rng_for(::Type{<:DiagOp{F}}) where {F} = global_rng_for(F)

# automatic basis conversion (and NaN-zeroing)
(*)(D::DiagOp{<:Field{B}}, f::Field) where {B} = D.diag .* B(f)
(\)(D::DiagOp{<:Field{B}}, f::Field) where {B} = nan2zero.(D.diag .\ B(f))

# broadcasting
struct DiagOpStyle{FS} <: AbstractArrayStyle{2} end
BroadcastStyle(::Type{D}) where {F<:Field,D<:DiagOp{F}} = DiagOpStyle{typeof(BroadcastStyle(F))}()
BroadcastStyle(::DiagOpStyle{FS1}, ::DiagOpStyle{FS2}) where {FS1,FS2} = DiagOpStyle{typeof(result_style(FS1(),FS2()))}()
BroadcastStyle(S::DiagOpStyle, ::DefaultArrayStyle{0}) = S
similar(bc::Broadcasted{DiagOpStyle{FS}}, ::Type{T}) where {FS,T} = Diagonal(similar(FS,T))
instantiate(bc::Broadcasted{<:DiagOpStyle}) = bc
diag_data(D::Diagonal) = D.diag
diag_data(x) = x
function copyto!(dest::DiagOp, bc::Broadcasted{Nothing})
    copyto!(dest.diag, map_bc_args(diag_data, bc))
    dest
end

# the generic versions of these kind of suck so we need these specializized
# versions:
(*)(x::Adjoint{<:Any,<:Field}, D::Diagonal) = (D*parent(x))'
(*)(x::Adjoint{<:Any,<:Field}, D::Diagonal, y::Field) = x*(D*y)
diag(L::DiagOp) = L.diag

# use getproperty here to ensure no basis conversion is done
getindex(D::DiagOp, s::Symbol) = Diagonal(getproperty(D.diag,s))

# for printing
size(::DiagOp{<:ImplicitField}) = ()
axes(::DiagOp{<:ImplicitField}, I) = OneTo(0)

# the generic version of this is prohibitively slow so we need this
hash(D::DiagOp, h::UInt64) = hash(D.diag, h)

### Derivative ops

# These ops just store the coordinate with respect to which a derivative is
# being taken, and each Field type F implements how to actually take the
# derivative.
#
# ∇diag{coord, covariance, prefactor} is a Field which represents the diagonal
# entries of the Fourier representation of the gradient operation, with `coord`
# specifiying the coordinate for the gradient direction (1 or 2), `covariance`
# being :covariant or :contravariant, and multiplied by a `prefactor` (generally
# 1 or -1 to keep track of if an adjoint has been taken). ∇[i] returns one of
# these wrapped in a Diagonal so that multipliying by these does the necessary
# basis conversion. ∇ (along with ∇ⁱ, and ∇ᵢ) are StaticVectors which can also
# be used with FieldVector algebra, e.g. ∇*f.
#
# Note: we are cheating a bit, because while ∇ is diagonal when a derivative is
# taken in Fourier space via FFT's, its tri-diagonal if taking a map-space
# derivative, so conceptually these should not be wrapped in a Diagonal. Its a
# bit ugly and I'd like to fix it, but since the auto-basis conversion mechanism
# is based on Diagonal, keep things this way for now. Individual field types
# just need to intercept the broadcast machinery at copyto! to implement the
# map-space derivative, e.g. see: flat_generic.jl.
#
# Also note: We define the components of vectors, including ∇, to be with respect to
# the _unnormalized_ covariant or contravariant basis vectors, hence ∇ⁱ = d/dxᵢ
# and ∇ᵢ = d/dxⁱ. This is different than with respect to the _normalized)
# covariant basis vectors, which, e.g., in spherical coordinates, gives the more
# familiar ∇ = (d/dθ, 1/sinθ d/dϕ), (but whose components are neither covariant
# nor contravariant).

struct ∇diag{coord, covariance, prefactor} <: ImplicitField{DerivBasis,Spin,Pix} end
struct ∇²diag <: ImplicitField{DerivBasis,Spin,Pix} end

# adjoint(D::Diagonal{<:∇diag}) calls conj(D.diag) and here lazily we keep track
# of that a conjugate was taken
conj(::∇diag{coord,covariance,prefactor}) where {coord,covariance,prefactor} = ∇diag{coord,covariance,-prefactor}()
-(::DiagOp{∇diag{coord,covariance,prefactor}}) where {coord,covariance,prefactor} = Diagonal(∇diag{coord,covariance,-prefactor}())

# Gradient vector which can be used with FieldVector algebra.
struct ∇Op{covariance,prefactor} <: StaticVector{2,Diagonal{Float32,∇diag{<:Any,covariance,prefactor}}} end
getindex(::∇Op{covariance,prefactor}, i::Int) where {covariance,prefactor} = Diagonal(∇diag{i,covariance,prefactor}())
-(::∇Op{covariance,prefactor}) where {covariance,prefactor} = ∇Op{covariance,-prefactor}()
const ∇ⁱ = ∇Op{:contravariant,1}()
const ∇ᵢ = ∇Op{:covariant,1}()
const ∇ = ∇ⁱ # ∇ is contravariant by default if not specified
const ∇² = Diagonal(∇²diag())

@doc doc"""
    gradhess(f)

Compute the gradient $g^i = \nabla^i f$, and the hessian, $H_j^{\,i} = \nabla_j \nabla^i f$.
"""
function gradhess(f)
    g = ∇ⁱ*f
    H = @SMatrix[∇ᵢ[1]*g[1] ∇ᵢ[2]*g[1]; ∇ᵢ[1]*g[2] ∇ᵢ[2]*g[2]]
    @namedtuple(g,H)
end

# this is not strictly true (∇[1] is generically a gradient w.r.t. the first
# coordinate, e.g. ∂θ), but this is useful shorthand to have for the flat-sky:
const ∂x = ∇[1]
const ∂y = ∇[2]

### FuncOp

# An Op which applies some arbitrary function to its argument.
# Transpose and/or inverse operations which are not specified will return an error.
@kwdef struct FuncOp <: ImplicitOp{Basis,Spin,Pix}
    op   = nothing
    opᴴ  = nothing
    op⁻¹ = nothing
    op⁻ᴴ = nothing
end
SymmetricFuncOp(;op=nothing, op⁻¹=nothing) = FuncOp(op,op,op⁻¹,op⁻¹)
FuncOp(op::Function) = FuncOp(op=op)
SymmetricFuncOp(op::Function) = SymmetricFuncOp(op=op)
*(L::FuncOp, f::Field) =
    L.op   != nothing ? L.op(f)   : error("op*f not implemented")
\(L::FuncOp, f::Field) =
    L.op⁻¹ != nothing ? L.op⁻¹(f) : error("op\\f not implemented")
*(f::Adjoint{<:Any,Field}, L::FuncOp) =
    L.opᴴ  != nothing ? L.opᴴ(f)  : error("opᴴ*f not implemented")
adjoint(L::FuncOp) = FuncOp(L.opᴴ,L.op,L.op⁻ᴴ,L.op⁻¹)
inv(L::FuncOp) = FuncOp(L.op⁻¹,L.op⁻ᴴ,L.op,L.opᴴ)
adapt_structure(to, L::FuncOp) = FuncOp(adapt(to, fieldvalues(L))...)

### BandPassOp

# An op which applies some bandpass, like a high or low-pass filter. This object
# stores the bandpass weights, Wℓ, and each Field type F should implement
# broadcast_data(::Type{F}, ::BandPassOp) to describe how this is actually
# applied.

abstract type HarmonicBasis <: Basislike end

struct BandPass{W<:InterpolatedCℓs} <: ImplicitField{HarmonicBasis,Spin,Pix}
    Wℓ::W
end
BandPassOp(ℓ,Wℓ) = Diagonal(BandPass(InterpolatedCℓs(promote(collect(ℓ),collect(Wℓ))...)))
BandPassOp(Wℓ::InterpolatedCℓs) = Diagonal(BandPass(Wℓ))
cos_ramp_up(length) = @. (cos($range(π,0,length=length))+1)/2
cos_ramp_down(length) = 1 .- cos_ramp_up(length)
HighPass(ℓ; Δℓ=50) = BandPassOp(ℓ:20000, [cos_ramp_up(Δℓ); ones(20000-ℓ-Δℓ+1)])
LowPass(ℓ; Δℓ=50) = BandPassOp(0:ℓ, [ones(ℓ-Δℓ+1); cos_ramp_down(Δℓ)])
MidPass(ℓmin, ℓmax; Δℓ=50) = BandPassOp(ℓmin:ℓmax, [cos_ramp_up(Δℓ); ones(ℓmax-ℓmin-2Δℓ+1); cos_ramp_down(Δℓ)])
MidPasses(ℓedges; Δℓ=10) = [MidPass(ℓmin-Δℓ÷2,ℓmax+Δℓ÷2; Δℓ=Δℓ) for (ℓmin,ℓmax) in zip(ℓedges[1:end-1],ℓedges[2:end])]

### ParamDependentOp

# A LinOp which depends on some parameters, θ.
# L(;θ...) recomputes the operator at a given set of parameters, but the
# operator can also be used as-is in which case it is evaluated at a fiducial θ
# (which is stored inside the operator when it is first constructed).

@doc doc"""
    ParamDependentOp(recompute_function::Function)

Creates an operator which depends on some parameters $\theta$ and can be
evaluated at various values of these parameters.

`recompute_function` should be a function which accepts keyword arguments for
$\theta$ and returns the operator. Each keyword must have a default value; the
operator will act as if evaluated at these defaults unless it is explicitly
evaluated at other parameters.

Example:

```julia
Cϕ₀ = Diagonal(...) # some fixed Diagonal operator
Cϕ = ParamDependentOp((;Aϕ=1)->Aϕ*Cϕ₀) # create ParamDependentOp

Cϕ(Aϕ=1.1) * ϕ   # Cϕ(Aϕ=1.1) is equal to 1.1*Cϕ₀
Cϕ * ϕ           # Cϕ alone will act like Cϕ(Aϕ=1) because that was the default above
```

Note: if you are doing parallel work, global variables referred to in the
`recompute_function` need to be distributed to all workers. A more robust
solution is to avoid globals entirely and instead ensure all variables are
"closed" over (and hence will automatically get distributed). This will happen
by default if defining the `ParamDependentOp` inside any function, or can be
forced at the global scope by wrapping everything in a `let`-block, e.g.:

```julia
Cϕ = let Cϕ₀=Cϕ₀
    ParamDependentOp((;Aϕ=1)->Aϕ*Cϕ₀)
end
```

After executing the code above, `Cϕ` is now ready to be (auto-)shipped to any workers
and will work regardless of what global variables are defined on these workers.
"""
struct ParamDependentOp{B, S, P, L<:LinOp{B,S,P}, F<:Function} <: ImplicitOp{B,S,P}
    op::L
    recompute_function::F
    parameters::Vector{Symbol}
end
function ParamDependentOp(recompute_function::Function)
    # invokelatest here allows creating a ParamDependent op which calls a
    # BinRescaledOp (eg this is the case for the mixing matrix G which depends
    # on Cϕ) from inside function. this would otherwise fail due to
    # BinRescaledOp eval'ed function being too new
    kwarg_names = get_kwarg_names(recompute_function)
    if endswith(string(kwarg_names[end]), "...") && !startswith(string(kwarg_names[end]),"_")
        kwarg_decl = empty!(kwarg_names) # to indicate it depends on anything
    end
    ParamDependentOp(Base.invokelatest(recompute_function), recompute_function, kwarg_names)
end
function (L::ParamDependentOp)(θ::NamedTuple)
    if depends_on(L,θ)
        # filtering out non-dependent parameters disabled until I can find a fix to:
        # https://discourse.julialang.org/t/can-zygote-do-derivatives-w-r-t-keyword-arguments-which-get-captured-in-kwargs/34553/8
        # dependent_θ = filter(((k,_),)->k in L.parameters, pairs(θ))
        L.recompute_function(;θ...)
    else
        L.op
    end
end
(L::ParamDependentOp)(;θ...) = L((;θ...))
*(L::ParamDependentOp, f::Field) = L.op * f
\(L::ParamDependentOp, f::Field) = L.op \ f
for F in (:inv, :pinv, :sqrt, :adjoint, :Diagonal, :diag, :simulate, :zero, :one, :logdet, :global_rng_for)
    @eval $F(L::ParamDependentOp) = $F(L.op)
end
getindex(L::ParamDependentOp, x) = getindex(L.op, x)
simulate(rng::AbstractRNG, L::ParamDependentOp) = simulate(rng, L.op)
depends_on(L::ParamDependentOp, θ) = depends_on(L, keys(θ))
depends_on(L::ParamDependentOp, θ::Tuple) = isempty(L.parameters) || any(L.parameters .∈ Ref(θ))
depends_on(L,                   θ) = false

adapt_structure(to, L::ParamDependentOp) =
    ParamDependentOp(adapt(to, L.op), adapt(to, L.recompute_function), L.parameters)

### LazyBinaryOp

# we use LazyBinaryOps to create new operators composed from other operators
# which don't actually evaluate anything until they've been multiplied by a
# field
struct LazyBinaryOp{F,A<:Union{LinOrAdjOp,Scalar},B<:Union{LinOrAdjOp,Scalar}} <: ImplicitOp{Basis,Spin,Pix}
    a::A
    b::B
    LazyBinaryOp(op,a::A,b::B) where {A,B} = new{op,A,B}(a,b)
end
# creating LazyBinaryOps
for op in (:+, :-, :*)
    @eval ($op)(a::ImplicitOrAdjOp,          b::ImplicitOrAdjOp)          = LazyBinaryOp($op,a,b)
    @eval ($op)(a::Union{LinOrAdjOp,Scalar}, b::ImplicitOrAdjOp)          = LazyBinaryOp($op,a,b)
    @eval ($op)(a::ImplicitOrAdjOp,          b::Union{LinOrAdjOp,Scalar}) = LazyBinaryOp($op,a,b)
    # explicit vs. lazy binary operations on Diagonals:
    @eval ($op)(D1::DiagOp{<:Field{B}},  D2::DiagOp{<:Field{B}})  where {B}     = Diagonal(broadcast($op,D1.diag,D2.diag))
    @eval ($op)(D1::DiagOp{<:Field{B1}}, D2::DiagOp{<:Field{B2}}) where {B1,B2} = LazyBinaryOp($op,D1,D2)
end
/(op::ImplicitOrAdjOp, n::Real) = LazyBinaryOp(/,op,n)
^(op::Union{ImplicitOrAdjOp,DiagOp{<:ImplicitField}}, n::Int) = LazyBinaryOp(^,op,n)
inv(op::Union{ImplicitOrAdjOp,DiagOp{<:ImplicitField}}) = LazyBinaryOp(^,op,-1)
-(op::ImplicitOrAdjOp) = -1 * op
pinv(op::LazyBinaryOp{*}) = pinv(op.b) * pinv(op.a)
# evaluating LazyBinaryOps
for op in (:+, :-)
    @eval *(lz::LazyBinaryOp{$op}, f::Field) = ($op)(lz.a * f, lz.b * f)
    @eval diag(lz::LazyBinaryOp{$op}) = ($op)(diag(lz.a), diag(lz.b))
    @eval adjoint(lz::LazyBinaryOp{$op}) = LazyBinaryOp(($op),adjoint(lz.a),adjoint(lz.b))
end
*(lz::LazyBinaryOp{/}, f::Field) = (lz.a * f) / lz.b
*(lz::LazyBinaryOp{*}, f::Field) = lz.a * (lz.b * f)
\(lz::LazyBinaryOp{*}, f::Field) = lz.b \ (lz.a \ f)
*(lz::LazyBinaryOp{^}, f::Field) = foldr((lz.b>0 ? (*) : (\)), fill(lz.a,abs(lz.b)), init=f)
adjoint(lz::LazyBinaryOp{*}) = LazyBinaryOp(*,adjoint(lz.b),adjoint(lz.a))
ud_grade(lz::LazyBinaryOp{op}, args...; kwargs...) where {op} = LazyBinaryOp(op,ud_grade(lz.a,args...;kwargs...),ud_grade(lz.b,args...;kwargs...))
adapt_structure(to, lz::LazyBinaryOp{op}) where {op} = LazyBinaryOp(op, adapt(to,lz.a), adapt(to,lz.b))
function diag(lz::LazyBinaryOp{*})
    _diag(x) = diag(x)
    _diag(x::Int) = x
    da, db = _diag(lz.a), _diag(lz.b)
    # if basis(da)!=basis(db)
    #     error("Can't take diag(A*B) where A::$(typeof(lz.a)) and B::$(typeof(lz.b)).")
    # end
    da .* db
end

### OuterProdOp

# an operator L which represents L = V*W'
# this could also be represented by a LazyBinaryOp, but this allows us to
# define a few extra functions like simulate or diag, which get used in various places
struct OuterProdOp{TV,TW} <: ImplicitOp{Basis,Spin,Pix}
    V::TV
    W::TW
end
OuterProdOp(V) = OuterProdOp(V,V)
_check_sym(L::OuterProdOp) = L.V === L.W ? L : error("Can't do this operation on non-symmetric OuterProdOp.")
pinv(L::OuterProdOp{<:LazyBinaryOp{*}}) = (_check_sym(L); OuterProdOp(pinv(L.V.a)' * pinv(L.V.b)'))
*(L::OuterProdOp, f::Field) = L.V * (L.W' * f)
\(L::OuterProdOp{<:LinOp,<:LinOp}, f::Field) = L.W' \ (L.V \ f)
adjoint(L::OuterProdOp) = OuterProdOp(L.W,L.V)
adapt_structure(to, L::OuterProdOp) = OuterProdOp((V′=adapt(to,L.V);), (L.V === L.W ? V′ : adapt(to,L.W)))
diag(L::OuterProdOp{<:Field{B},<:Field}) where {B} = L.V .* conj.(B(L.W))
*(D::DiagOp{<:Field{B}}, L::OuterProdOp{<:Field{B},<:Field{B}}) where {B} = OuterProdOp(diag(D) .* L.V, L.W)
*(L::OuterProdOp{<:Field{B},<:Field{B}}, D::DiagOp{<:Field{B}}) where {B} = OuterProdOp(L.V, L.W .* diag(D))
tr(L::OuterProdOp{<:Field{B},<:Field{B}}) where {B} = dot(L.V, L.W)

### BinRescaledOp

"""
    BinRescaledOp(C₀, Cbins, θname::Symbol)

Create a [`ParamDependentOp`](@ref) which has a parameter named `θname` which is
an array that controls the amplitude of bandpowers in bins given by `Cbins`.

For example, `BinRescaledOp(C₀, [Cbin1, Cbin2], :A)` creates the operator:

    ParamDependentOp( (;A=[1,1], _...) -> C₀ + (A[1]-1) * Cbin1 + (A[2]-1) * Cbin2 )

where `C₀`, `Cbin1`, and `Cbin2` should be some `LinOp`s. Note `Cbins` are
directly the power which is added, rather than a mask.

The resulting operator is differentiable in `θname`.

"""
function BinRescaledOp(C₀, Cbins, θname::Symbol)
    # for some reason, if I eval this into CMBLensing as opposed to Main, it
    # doesn't get shipped to workers correctly.
    # see also: https://discourse.julialang.org/t/closure-not-shipping-to-remote-workers-except-from-main
    @eval Main begin
        # ensure Cbins is a tuple and not Array so that `adapt` works recursively through it
        let C₀=$C₀, T=$(real(eltype(C₀))), Cbins=$(tuple(Cbins...))
            $ParamDependentOp(function (;($θname)=$(ones(Int,length(Cbins))), _...)
                $(Expr(:call, :+, :C₀, [:(T($θname[$i] - 1) * Cbins[$i]) for i=1:length(Cbins)]...))
            end)
        end
    end
end


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

spy(randn((10, 10)), Coord.cartesian(fixed=true))


text:
export CircleRepresentationDrawing, RainbowDrawing

"""
`CircleRepresentationDrawing(list, fontsize=8)` draws the circle representation
implicit in `list`. Requires `PyPlot`. Not exporting this for that
reason.
"""
function CircleRepresentationDrawing(list::Array{T,1}, font_size = 8) where {T}
    font_size = 8
    positions = SimpleGraphRepresentations.list2locs(list)
    nn = length(list)
    verts = collect(keys(positions))
    rad = 2
    theta = collect(0:360) * pi / 180
    x = rad * map(cos, theta)
    y = rad * map(sin, theta)
    plot(x, y, linestyle = :dot, color = "black")

    for v in verts
        p = positions[v]
        a = (p[1] - 1) * 2 * pi / nn
        b = (p[2] - 1) * 2 * pi / nn
        x = rad * [cos(a), cos(b)]
        y = rad * [sin(a), sin(b)]
        plot!(x, y, color = "black")

        if font_size > 0
            mu = 1.1
            for k = 1:2
                # text(mu*x[k],mu*y[k],string(v))
                annotate!((mu * x[k], mu * y[k], string(v), font_size))
            end
        end

    end
    finish()
end

"""
`CircleRepresentationDrawing(str::String)` uses `str` as the
list of symbols.
"""
function CircleRepresentationDrawing(str::String, labels::Bool = true)
    list = [ch for ch in str]
    CircleRepresentationDrawing(list, labels)
end

function semicircle(x1::Number, x2::Number, diam = 2)
    ctr = (x1 + x2) / 2
    rad = abs(x1 - x2) / 2

    angles = collect(0:180) * (pi / 180)
    x = rad * [cos(t) for t in angles] .+ ctr
    y = rad * [sin(t) for t in angles]
    plot!(x, y, color = "black", linestyle = :dot)
    plot!([x[1]], [y[1]], marker = diam, markercolor = "black")
    plot!([x[end]], [y[end]], marker = diam, markercolor = "black")
end

"""
`RainbowDrawing(lst,font_size=8)` creates a rainbow representation
of a circle graph from a list of the positions.
"""
function RainbowDrawing(list::Array{T,1}, font_size::Int = 8) where {T}
    positions = SimpleGraphRepresentations.list2locs(list)
    nn = length(list)
    verts = collect(keys(positions))

    x = [0, nn + 1]
    y = [0, 0]
    plot(x, y, color = "black")

    gap = -0.5
    for v in verts
        p = positions[v]
        semicircle(p[1], p[2])
        if font_size > 0
            annotate!(p[1], gap, string(v), font_size)
            annotate!(p[2], gap, string(v), font_size)
        end
    end

    finish()
end


text:
using PackageCompiler

const EXT = if Sys.isapple()
    ".dylib"
elseif Sys.iswindows()
    ".dll"
elseif Sys.islinux()
    ".so"
else
    @warn "This operating system may not be supported."
    ".so"
end

PackageCompiler.create_sysimage(
    [:JuliaFormatter, :JSON];
    precompile_execution_file=joinpath(@__DIR__, "precompile.jl"),
    sysimage_path=joinpath(@__DIR__, "juliaformatter$EXT")
)


text:
mutable struct Polynomial <: AbstractComponent
    coeff::Vector{Parameter}

    function Polynomial(args...)
        coeff = [Parameter(arg) for arg in args]
        return new(coeff)
    end
end

function evaluate!(buffer, comp::Polynomial, x::AbstractDomain{1},
                   coeffs...)
    buffer .= coeffs[1]
    for deg in 1:length(coeffs)-1
        buffer .+= x.^deg .* coeffs[deg+1]
    end
end


text:
using Derivatives, InplaceRealFFT
using LESFilter: rfftfreq, fftfreq
using Base.Test

#=
Test function: f = sin(3x)*cos(3y)*sin(4z)
=#
for (lx,ly,lz) in ((1.,1.,1.),(2.,2.,2.),(3.,3.,3.))

# lx = 1.0
# ly = 1.0
# lz = 1.0

nx = 32
ny = 32
nz = 32

x = reshape(linspace(0,lx*2π*(1-1/nx),nx),(nx,1,1))
y = reshape(linspace(0,ly*2π*(1-1/ny),ny),(1,ny,1))
z = reshape(linspace(0,lz*2π*(1-1/nz),nz),(1,1,nz))

field = @. sin(3*x) * cos(3*y) * sin(4*z)

correct(f::typeof(Derivatives.dx)) = @. 3*cos(3*x) * cos(3*y) * sin(4*z)
correct(f::typeof(Derivatives.dy)) = @. sin(3*x) * -3*sin(3*y) * sin(4*z)
correct(f::typeof(Derivatives.dz)) = @. sin(3*x) * cos(3*y) * 4*cos(4*z)

correct2(f::typeof(Derivatives.dx)) = @. -9*sin(3*x) * cos(3*y) * sin(4*z)
correct2(f::typeof(Derivatives.dy)) = @. sin(3*x) * -9*cos(3*y) * sin(4*z)
correct2(f::typeof(Derivatives.dz)) = @. sin(3*x) * cos(3*y) * -16*sin(4*z)

correct3(f::typeof(Derivatives.dx)) = @. -27*cos(3*x) * cos(3*y) * sin(4*z)
correct3(f::typeof(Derivatives.dy)) = @. sin(3*x) * 27*sin(3*y) * sin(4*z)
correct3(f::typeof(Derivatives.dz)) = @. sin(3*x) * cos(3*y) * -64*cos(4*z)

for (c,n) in zip((correct,correct2,correct3),(1,2,3))
  for (f,l) in zip((Derivatives.dx,Derivatives.dy,Derivatives.dz),(lx,ly,lz))
    @test f(field,l,n) ≈ c(f)
  end
end

field2 = PaddedArray(field)

for (c,n) in zip((correct,correct2,correct3),(1,2,3))
  for (f,l) in zip((Derivatives.dx,Derivatives.dy,Derivatives.dz),(lx,ly,lz))
    @test real(f(field,l,n)) ≈ c(f)
  end
end

rfft!(field2)

@inferred Derivatives.loopdx!(complex(field2),rfftfreq(nx,lx))
@inferred Derivatives.loopdy!(complex(field2),fftfreq(ny,ly))
@inferred Derivatives.loopdz!(complex(field2),fftfreq(nz,lz))

end


text:
using Distributions, Econometrics
function main()
# sample is from exponential, prior is lognormal, proposal is random walk lognormal
y = rand(Exponential(3.0),30)
# set prior, likelihood and proposal
Prior = θ -> pdf(LogNormal(1.0,1.0), θ)
lnL = θ -> sum(logpdf.(Ref(Exponential(θ)),y))
tuning = 0.5
Proposal = θ -> rand(LogNormal(log(θ),tuning))
# get the chain, plot posterior, and descriptive stats
chain = mcmc(1.0, 100000, 10000, Prior, lnL, Proposal) # start value, chain length, and burnin
p = npdensity(chain[:,1]) # nonparametric plot of posterior density
plot!(p, title="posterior density: true value = 3.0") # add a title
display(p)
dstats(chain)
return
end
main()

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
################################################################################
#
#  Lifting again
#
################################################################################

mutable struct Emat
  R
  n
  i
  j
  a
end

function matrix(e::Emat)
  z = identity_matrix(e.R, e.n)
  z[e.i, e.j] = e.a
  return z
end

inv(e::Emat) = Emat(e.R, e.n, e.i, e.j, -e.a)

function elementary_matrix(R, n, i, j, a)
  M = identity_matrix(R, n)
  M[i, j] = a
  return M
end

function _lift2(MM)
  M = deepcopy(MM)
  R = base_ring(M)
  S = _base_ring(R)
  @assert det(M) == 1
  k = 1
  n = nrows(M)
  left = []
  left2 = []
  while k < ncols(M)
    @assert det(M) == 1
    # scan the first column for a unit
    v = Int[euclid(M[i, k]) for i in k:n]
    #@show v
    if isone(minimum(abs.(v)))
      l = findfirst(isone, v) + (k - 1)
      #@show k, M
      #@show isunit(M[l, k])
      #@show M[l, k]
      if l != k
        @assert l isa Int
        #@show l
        b = M[l, k]
        @assert det(M) == 1
        E1 = elementary_matrix(R, n, k, l, one(R))
        E2 = elementary_matrix(R, n, l, k, -one(R))
        E3 = elementary_matrix(R, n, k, l, one(R))
        M = (E1 * E2 * E3) * M
        @assert det(M) == 1
        # push first E3
        push!(left, elementary_matrix(S, n, k, l, one(S)))
        @assert map_entries(R, left[end]) == E3
        push!(left2, Emat(S, n, k, l, one(S)))
        push!(left, elementary_matrix(S, n, l, k, -one(S)))
        push!(left2, Emat(S, n, l, k, -one(S)))
        @assert map_entries(R, left[end]) == E2
        push!(left, elementary_matrix(S, n, k, l, one(S)))
        push!(left2, Emat(S, n, k, l, one(S)))
        @assert map_entries(R, left[end]) == E1
      end
      @assert det(M) == 1
      @assert isunit(M[k, k])
      for j in (k+1):n
        if !iszero(M[j, k])
          N = deepcopy(M)
          q = divexact(M[j, k], M[k, k])
          E = elementary_matrix(R, n, j, k, -q)
          push!(left, elementary_matrix(S, n, j, k, lift(-q)))
          push!(left2, Emat(S, n, j, k, lift(-q)))
          @assert map_entries(R, left[end]) == E
          for o in 1:n
            M[j, o] = M[j, o] - q * M[k, o]
          end
          @assert E * N == M
        end
      end
      k += 1
      #@show M
    else # no unit there # I could do this in one go by putting a unit in the first position
      _, l = findmin(abs.(v))
      l = l + (k - 1)
      #@show l
      #@show euclid(M[l, k])
      #@show v
      for j in k:n
        if j == l
          continue
        end
        fl, _ = divides(M[j, k], M[l, k])
        if !fl
          #@show M[j, k], M[l, k]
          N = deepcopy(M)
          #@show euclid(M[j, k])
          q, r = divrem(M[j, k], M[l, k])
          #@show euclid(r)
          E = elementary_matrix(R, n, j, l, -q)
          for o in 1:n
            M[j, o] = M[j, o] - q * M[l, o]
          end
          push!(left, elementary_matrix(S, n, j, l, lift(-q)))
          push!(left2, Emat(S, n, j, l, lift(-q)))
          @assert map_entries(R, left[end]) == E
          @assert E * N == M
          break
        end
      end
    end
    @assert det(M) == 1
  end
  #println("M should now be lower triangular")
  #@show M
  @assert det(M) == 1
  # Now M is lower triangular with units on the diagonal
  #@show "here"
  for i in n:-1:2
    for j in (i - 1):-1:1
      # set jth row to jth row - M[k, j]/M[k, k] * jth row
      #@show isunit(M[i, i])
      #@show M[i, i]
      q = -divexact(M[j, i], M[i, i])
      @assert (-q) * M[i, i] == M[j, i]
      E = elementary_matrix(R, n, j, i, q)
      push!(left, elementary_matrix(S, n, j, i, lift(q)))
      push!(left2, Emat(S, n, j, i, lift(q)))
      @assert map_entries(R, left[end]) == E
      N = deepcopy(M)
      M[j, i] = M[j, i] + q * M[i, i]
      @assert E * N == M
    end
  end
  @assert det(M) == 1
  # Now sort the diagonal
  # We use this neat trick of Rosenberg, p. 60
  for i in 1:(n - 1)
    a = inv(M[i, i])
    ainv = M[i, i]
    # we multiply with (1,...1,a,a^-1,1,...1)
    E1 = elementary_matrix(R, n, i, i + 1, a)
    E2 = elementary_matrix(R, n, i + 1, i, -ainv)
    E3 = E1
    #E4 = identity_matrix(R, n)
    #E4[i, i] = 0
    #E4[i, i + 1] = -1
    #E4[i + 1, i + 1] = 0
    #E4[i + 1, i] = 1
    E5 = elementary_matrix(R, n, i, i + 1, -one(R))
    E6 = elementary_matrix(R, n, i + 1, i, one(R))
    E4 = E6

    #E4S = identity_matrix(S, n)
    #E4S[i, i] = 0
    #E4S[i, i + 1] = -1
    #E4S[i + 1, i + 1] = 0
    #E4S[i + 1, i] = 1
    E1S = elementary_matrix(S, n, i, i + 1, lift(a))
    E1SS = Emat(S, n, i, i + 1, lift(a))
    E2S = elementary_matrix(S, n, i + 1, i, lift(-ainv))
    E2SS = Emat(S, n, i + 1, i, lift(-ainv))
    E3S = E1S
    E3SS = E1SS
    E5S = elementary_matrix(S, n, i, i + 1, -one(S))
    E5SS = Emat(S, n, i, i + 1, -one(S))
    E6S = elementary_matrix(S, n, i + 1, i, one(S))
    E6SS = Emat(S, n, i + 1, i, one(S))
    E4S = E6S
    E4SS = E6SS
    push!(left, E6S)
    push!(left, E5S)
    push!(left, E4S)
    push!(left, E3S)
    push!(left, E2S)
    push!(left, E1S)
    push!(left2, E6SS)
    push!(left2, E5SS)
    push!(left2, E4SS)
    push!(left2, E3SS)
    push!(left2, E2SS)
    push!(left2, E1SS)

    M = (E1 * E2 * E3 * E4 * E5 * E6) * M
    @assert det(M) == 1
  end

  @assert isone(M)

  for i in 1:length(left)
    @assert matrix(left2[i]) == left[i]
  end

  @assert map(R, prod(reverse(left))) * MM == 1

  return prod(matrix.(inv.(left2)))

  return M, left
end

euclid(n::nmod) = gcd(n.data, modulus(parent(n)))
euclid(n::Nemo.fmpz_mod) = gcd(n.data, modulus(parent(n)))

function Base.divrem(n::T, m::T) where T <: Union{nmod,Nemo.fmpz_mod}
  @assert !iszero(m)
  R = parent(n)
  e = euclid(m)
  if isone(e)
    fl, q = divides(n, m)
    @assert fl
    return q, zero(R)
  end

  cp = coprime_base(fmpz[n.data, m.data, modulus(R)])::Array{fmpz, 1}

  q = Array{Tuple{fmpz, fmpz}, 1}()
  for i=1:length(cp)
    v = valuation(modulus(R), cp[i])::Int
    if v != 0
      pk = cp[i]^v
      nv = iszero(n.data % pk) ? inf : valuation(n.data % pk, cp[i])
      mv = iszero(m.data % pk) ? inf : valuation(m.data % pk, cp[i])
      if nv < mv
        push!(q, (pk, 0))
      else
        if nv === inf
          push!(q, (pk, 1))
        else
          push!(q, (pk, divexact(n.data % pk, cp[i]^nv)))
        end
      end
    end
  end
  qq = R(crt([x[2] for x = q], [x[1] for x = q])::fmpz)::T
  rr = n-qq*m
  @assert n == qq*m+rr
  @assert rr == 0 || euclid(rr) < e
  return (qq,rr)::Tuple{T, T}
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@reexport module LShaped

# Standard library
using LinearAlgebra
using SparseArrays
using Distributed
using Printf

# External libraries
using Parameters
using JuMP
using MathOptInterface
using StochasticPrograms
using StochasticPrograms: AcceptableTermination
using StochasticPrograms: UnspecifiedInstantiation, StageDecompositionStructure, AbstractScenarioProblems, ScenarioProblems, DistributedScenarioProblems, DecisionChannel
using StochasticPrograms: AbstractExecution, Serial, Synchronous, Asynchronous
using StochasticPrograms: AbstractStructuredOptimizer, set_master_optimizer!, set_subproblem_optimizer!
using StochasticPrograms: DecisionMap, all_decisions, set_decision!, set_stage!, SingleDecisionSet, NoSpecifiedConstraint
using StochasticPrograms: all_known_decisions, update_known_decisions!, KnownValuesChange, Known, KnownDecision
using StochasticPrograms: add_subscript
using StochasticPrograms: AbstractPenaltyTerm, Quadratic, InfNorm, ManhattanNorm, initialize_penaltyterm!, update_penaltyterm!, remove_penalty!, remove_penalty_variables!, remove_penalty_constraints!
using ProgressMeter
using Clustering

import Base: show, put!, wait, isready, take!, fetch, zero, +, length, size
import StochasticPrograms: supports_structure, num_iterations, default_structure, check_loadable, load_structure!, restore_structure!, optimize!, optimizer_name, master_optimizer, subproblem_optimizer, num_subproblems, remove_penalty_variables!, remove_penalty_constraints!, relax_decision_integrality

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const CI = MOI.ConstraintIndex
const CutConstraint = CI{AffineDecisionFunction{Float64}, MOI.GreaterThan{Float64}}

export
    AbstractLShapedAttribute,
    LShapedAlgorithm,
    num_cuts,
    num_iterations,
    FeasibilityStrategy,
    IgnoreFeasibility,
    FeasibilityCuts,
    Regularizer,
    RegularizationParameter,
    RawRegularizationParameter,
    set_regularization_attribute,
    set_regularization_attributes,
    DontRegularize,
    NoRegularization,
    RegularizedDecomposition,
    WithRegularizedDecomposition,
    RD,
    WithRD,
    TrustRegion,
    WithTrustRegion,
    TR,
    WithTR,
    LevelSet,
    WithLevelSets,
    LV,
    WithLV,
    IntegerStrategy,
    IgnoreIntegers,
    CombinatorialCuts,
    Convexification,
    Gomory,
    LiftAndProject,
    CuttingPlaneTree,
    Consolidator,
    DontConsolidate,
    Consolidate,
    Consolidation,
    ConsolidationParameter,
    RawConsolidationParameter,
    set_consolidation_attribute,
    set_consolidation_attributes,
    at_tolerance,
    num_iters,
    tolerance_reached,
    Aggregator,
    AggregationParameter,
    RawAggregationParameter,
    set_aggregation_attribute,
    set_aggregation_attributes,
    DontAggregate,
    NoAggregation,
    PartialAggregate,
    PartialAggregation,
    Aggregate,
    FullAggregation,
    DynamicAggregate,
    DynamicAggregation,
    ClusterAggregate,
    ClusterAggregation,
    HybridAggregate,
    HybridAggregation,
    GranulatedAggregate,
    GranulatedAggregation,
    SelectUniform,
    SelectDecaying,
    SelectRandom,
    SelectClosest,
    SortByReference,
    StaticCluster,
    ClusterByReference,
    Kmedoids,
    Hierarchical,
    absolute_distance,
    angular_distance,
    spatioangular_distance

# Include files
include("types/types.jl")
include("integer/integer.jl")
include("consolidators/consolidation.jl")
include("aggregators/aggregation.jl")
include("regularizers/regularization.jl")
include("execution/execution.jl")
include("solver.jl")
include("MOI_wrapper.jl")

end # module


text:
using SimplexSplitting

@testset "Triangulation" begin
    @testset "From embedding" begin
        npts = 15; cov = 0.4; emb = gaussian_embedding(npts, cov)
        t = triang_from_embedding(emb)

        # Verify that triangulation is created
        @test typeof(t) == Triangulation

        # Check functions to extract parts of the triangulation
        n_simplices = size(t.simplex_inds, 1)
        @test length(get_simplices(t)) == n_simplices
        @test length(get_imagesimplices(t)) == n_simplices
        @test typeof(find_simplex(t, 1)) == Simplex
        @test typeof(find_imsimplex(t, 1)) == Simplex
    end
end

max_stars_repo_name: kahaaga/SimplexSplitting

max_issues_repo_name: kahaaga/SimplexSplitting

max_forks_repo_name: kahaaga/SimplexSplitting


text:
@testset "Null Time Change" begin

    testEvents = sort(rand(100)*100)
    null = HawkesProcesses.time_change_null(testEvents, 100)

    @test length(null) == 2
    inferTests && @inferred HawkesProcesses.time_change_null(testEvents, 100)

end

@testset "Hawkes Time Change" begin
    bg = 0.5
    kappa = 0.5
    kernD = Distributions.Exponential(1/0.5)
    events = sort(rand(10))

    res = HawkesProcesses.time_change_hawkes(events, bg, kappa, kernD)
    @test length(res) == 2
    inferTests && @inferred HawkesProcesses.time_change_hawkes(events, bg, kappa, kernD)

end

@testset "Internal Function" begin

    res = HawkesProcesses._time_change(sort(rand(100)))

    @test length(res) == 2
    inferTests && @inferred HawkesProcesses._time_change(sort(rand(100)))

end


text:
export IsothermalEnergyDistribution

"Initial isothermal energy distribution."
struct IsothermalEnergyDistribution{N} <: InitialEnergyDistribution{N}
    "Temperature."
    T::Float

    IsothermalEnergyDistribution{N}(T::Number) where {N} =
        @check(T > 0) && new(T)
end

Base.show(io::IO, ::Type{IsothermalEnergyDistribution{N}}) where {N} =
    print(io, "IsothermalEnergyDistribution")

function (distribution::IsothermalEnergyDistribution{N})(
    positions::AbstractVector{SVectorF{N}},
    mass_densities::AbstractVector{Float},
    eos::EquationOfState,
) where {N}
    @checkeq(length(positions), length(mass_densities))
    fill(specific_energy_from(eos, distribution.T), length(positions))
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
scene_spec =
[
    (:fill, Any, nothing),
    (:fillOpacity, Any, nothing),
    (:stroke, Any, nothing),
    (:strokeOpacity, Any, nothing),
    (:strokeWidth, Any, nothing),
    (:strokeDash, Any, nothing),
    (:strokeDashOffset, Any, nothing),
]

primitivefactory(:VegaScene, scene_spec)

max_stars_repo_name: jpthor/vega