text:
using DimensionalData, Dates, Test, BenchmarkTools

using DimensionalData: val, basetypeof, slicedims, dims2indices, mode, dimsmatch,
      @dim, reducedims, XDim, YDim, ZDim, commondims, dim2key, key2dim, dimstride,
      _call, AlwaysTuple, MaybeFirst, _wraparg, _reducedims, combinedims

@testset "dimsmatch" begin
    @test (@inferred dimsmatch(Y(), Y())) == true
    @test (@inferred dimsmatch(X(), Y())) == false
    @test (@inferred dimsmatch(Y(), Y)) == true
    @test (@inferred dimsmatch(X(), Y)) == false
    @test (@inferred dimsmatch(Y, Y())) == true
    @test (@inferred dimsmatch(X, Y())) == false
    @test (@inferred dimsmatch(X, XDim)) == true
    @test (@inferred dimsmatch(Y, ZDim)) == false
    @test (@inferred dimsmatch(ZDim, Dimension)) == true
    @test (@inferred dimsmatch((Z(), ZDim), (ZDim, Dimension))) == true
    @test (@inferred dimsmatch((Z(), ZDim), (ZDim, XDim))) == false

    @test (@inferred dimsmatch(XDim, nothing)) == false
    @test (@inferred dimsmatch(X(), nothing)) == false
    @test (@inferred dimsmatch(nothing, ZDim)) == false
    @test (@inferred dimsmatch(nothing, Z())) == false
    @test (@inferred dimsmatch(nothing, nothing)) == false

    @test (@ballocated dimsmatch(ZDim, Dimension)) == 0
    @test (@ballocated dimsmatch((Z(), ZDim), (ZDim, XDim))) == 0
end

@dim Tst

@testset "key2dim" begin
    @test key2dim(:test) == Dim{:test}()
    @test key2dim(:X) == X()
    @test key2dim(:x) == Dim{:x}()
    @test key2dim(:Ti) == Ti()
    @test key2dim(:ti) == Dim{:ti}()
    @test key2dim(:Tst) == Tst()
    @test key2dim(Ti) == Ti
    @test key2dim(Ti()) == Ti()
    @test key2dim(Val{TimeDim}()) == Val{TimeDim}()
    @test key2dim((:test, Ti, Ti())) == (Dim{:test}(), Ti, Ti())
end

@testset "dim2key" begin
    @test dim2key(X()) == :X
    @test dim2key(Dim{:x}) == :x
    @test dim2key(Tst()) == :Tst
    @test dim2key(Dim{:test}()) == :test
    @test (@ballocated dim2key(Dim{:test}())) == 0
    @test dim2key(Val{TimeDim}()) == :TimeDim
end

@testset "_wraparg" begin
    @test _wraparg(X()) == (X(),)
    @test _wraparg((X, Y,), X,) == ((Val(X), Val(Y),), Val(X),)
    @test _wraparg((X, X(), :x), Y, Y(), :y) ==
        ((Val(X), X(), Dim{:x}()), Val(Y), Y(), Dim{:y}())
    f1 = () -> _wraparg((X(), Y()), X(), Y())
    f2 = () -> _wraparg((:x, :y), :x, :y)
    f3 = () -> _wraparg((X(), Y()), X, Y)
    f4 = () -> _wraparg((X, Y), X, Y)
    f5 = () -> _wraparg((X, X(), :x), Y, Y(), :y)
    f6 = () -> _wraparg((:x, :y, :z), :x, :y)
    f7 = () -> _wraparg((:x, :y), :x, :y, :z)
    f8 = () -> _wraparg((:x, :y, :z, :a, :b), :x, :y, :z, :a, :b)
    f9 = () -> _wraparg((:x, :y, :z, :a, :b, :c, :d, :e, :f, :g, :h), :x, :y, :z, :a, :b, :c, :d, :e, :f, :g, :h)
    @inferred f1()
    @inferred f2()
    @inferred f3()
    @inferred f4()
    @inferred f5()
    @inferred f6()
    @inferred f7()
    @inferred f8()
    @test (@inferred f9()) == (DimensionalData.key2dim.(((:x, :y, :z, :a, :b, :c, :d, :e, :f, :g, :h))),
                        DimensionalData.key2dim.(((:x, :y, :z, :a, :b, :c, :d, :e, :f, :g, :h)))...)
end

@testset "_call" begin
    @test _call((f, args...) -> args, AlwaysTuple(), (Z, :a, :b), Ti, XDim, :x) ==
          _call((f, args...) -> args, MaybeFirst(), (Z, :a, :b), Ti, XDim, :x) ==
          ((Val(Z), Dim{:a}(), Dim{:b}()), (Val(Ti), Val(XDim), Dim{:x}()))
    @test _call((f, args...) -> args, MaybeFirst(), (Z, :a, :b), Ti) ==
        (Val(Z), Dim{:a}(), Dim{:b}())
    @testset "_call" begin
        f1 = t -> _call((f, args...) -> args, t, (X, :a, :b), (TimeDim, X(), :a, :b, Ti))
        @inferred f1(AlwaysTuple())
        @test (@ballocated $f1(AlwaysTuple())) == 0
    end
end

@testset "sortdims" begin
    dimz = (X(), Y(), Z(), Ti())
    @test (@inferred sortdims((Y(1:2), X(1)), dimz)) == (X(1), Y(1:2), nothing, nothing)
    @test (@inferred sortdims((Ti(1),), dimz)) == (nothing, nothing, nothing, Ti(1))
    @test (@inferred sortdims((Z(), Y(), X()), dimz)) == (X(), Y(), Z(), nothing)
    @test (@inferred sortdims(dimz, (Y(), Z()))) == (Y(), Z())
    @test (@inferred sortdims((Y(), X(), Z(), Ti()), dimz)) == (X(), Y(), Z(), Ti())
    @test (@ballocated sortdims((Y(), X(), Z(), Ti()), $dimz)) == 0
    f1 = (dimz) -> sortdims((Y, X, Z, Ti), dimz)
    @test (@inferred f1(dimz)) == (X, Y, Z, Ti)
    @ballocated $f1($dimz)
    @test (@ballocated $f1($dimz)) == 0
    f2 = (dimz) -> sortdims(dimz, (Y, X, Z, Ti))
    @test (@inferred f2(dimz)) == (Y(), X(), Z(), Ti())
    @test (@ballocated $f2($dimz)) == 0
    # Val
    @inferred sortdims(dimz, (Val{Y}(), Val{Ti}(), Val{Z}(), Val{X}()))
    @test (@ballocated sortdims($dimz, (Val{Y}(), Val{Ti}(), Val{Z}(), Val{X}()))) == 0
    # Transformed
    @test (@inferred sortdims((Y(1:2; mode=Transformed(identity, Z())), X(1)), (X(), Z()))) ==
                             (X(1), Y(1:2; mode=Transformed(identity, Z())))
    # Abstract
    @test sortdims((Z(), Y(), X()), (XDim, TimeDim)) == (X(), nothing)
    # Repeating
    @test sortdims((Z(:a), Y(), Z(:b)), (YDim, Z(), ZDim)) == (Y(), Z(:a), Z(:b))
end

a = [1 2 3; 4 5 6]
da = DimArray(a, (X((143, 145)), Y((-38, -36))))
dimz = dims(da)

@testset "dims" begin
    dimz = dims(da)
    @test dims(X()) == X()
    @test (@inferred dims(dimz, X())) isa X
    @test (@ballocated dims($dimz, X())) == 0
    @test (@inferred dims(dimz, Y)) isa Y
    @test (@ballocated dims($dimz, Y)) == 0
    @test dims(dimz, (X(), Y())) isa Tuple{<:X,<:Y}
    @test (@ballocated dims($dimz, (X(), Y()))) == 0
    @test dims(dimz, (XDim, YDim)) isa Tuple{<:X,<:Y}
    # It's hard to make this infer - Types in a tuple are UnionAll
    f1 = (da) -> dims(da, (XDim, YDim))
    @test (@inferred f1(da)) isa Tuple{<:X,<:Y}
    @test (@ballocated $f1($dimz)) == 0

    @test dims(da, X()) isa X
    @test dims(da, Z()) isa Nothing
    @test (@inferred dims(da, XDim, YDim)) isa Tuple{<:X,<:Y}
    @test (@ballocated dims($da, XDim, YDim)) == 0
    @test dims(da, ()) == ()
    @test dims(dimz, 1) isa X
    @test dims(dimz, (2, 1)) isa Tuple{<:Y,<:X}
    # Mixed Int/dim indexing is no longer supported
    @test_broken dims(dimz, (2, Y)) isa Tuple{<:Y,<:Y}

    @testset "with Dim{X} and symbols" begin
        A = DimArray(zeros(4, 5), (:one, :two))
        @test dims(A) ==
            (Dim{:one}(Base.OneTo(4), NoIndex(), NoMetadata()),
             Dim{:two}(Base.OneTo(5), NoIndex(), NoMetadata()))
        @test dims(A, :two) == Dim{:two}(Base.OneTo(5), NoIndex(), NoMetadata())
    end

    # @test_throws ArgumentError dims(da, Ti)
    # @test_throws ArgumentError dims(dimz, Ti)
    @test_throws ArgumentError dims(nothing, X)

    @test dims(dimz) === dimz
    @test dims(dimz, X) === dimz[1]
    @test dims(dimz, Y) === dimz[2]
    @test typeof(dims(da)) ==
        Tuple{X{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},
              Y{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}}
end

@testset "commondims" begin
    @test commondims(da, X) == (dims(da, X),)
    # Dims are always in the base order
    @test (@inferred commondims(da, (Y(), X()))) == dims(da, (X, Y))
    @test (@ballocated commondims($da, (Y(), X()))) == 0
    f1 = (da) -> commondims(da, (X, Y))
    @test (@inferred f1(da)) == dims(da, (X, Y))
    @test (@ballocated $f1($da)) == 0
    @test (@inferred commondims(da, X, Y)) == dims(da, (X, Y))
    @test (@ballocated commondims($da, X, Y)) == 0
    @test basetypeof(commondims(da, DimArray(zeros(5), Y()))[1]) <: Y

    @testset "with Dim{X} and symbols" begin
        dimz = Dim{:a}(), Y(), Dim{:b}()
        f1 = (dimz) -> commondims((Dim{:a}(), Dim{:b}()), (:a, :c))
        @test f1(dimz) == (Dim{:a}(),)
        @test (@ballocated $f1(dimz)) == 0
        f2 = (dimz) -> commondims(dimz, (Y, Dim{:b}))
        @test f2(dimz) == (Y(), Dim{:b}())
        f3 = (dimz) -> commondims(dimz, Y, :b)
        @test (@inferred f3(dimz)) == (Y(), Dim{:b}())
        @test (@ballocated $f3($dimz)) == 0
        f4 = (dimz) -> commondims(dimz, Y, Z, Dim{:b})
        @test (@inferred f4(dimz)) == (Y(), Dim{:b}())
        @test (@ballocated $f4($dimz)) == 0
    end

    @testset "with abstract types" begin
        @test commondims((Z(), Y(), Ti()), (ZDim, Dimension)) == (Z(), Y())
        @test commondims(>:, (TimeDim, YDim), (Z(), Ti(), Dim{:a})) == (TimeDim,)
    end
end

@testset "dimnum" begin
    @test dimnum(da, Y()) == dimnum(da, 2) == 2
    @test (@ballocated dimnum($da, Y())) == 0
    @test dimnum(da, X) == 1
    @test (@ballocated dimnum($da, X)) == 0
    @test dimnum(da, (Y(), X())) == (2, 1)
    @ballocated dimnum($da, (Y(), X()))
    @test (@ballocated dimnum($da, (Y(), X()))) == 0
    f1 = (da) -> dimnum(da, (Y, X()))
    @test (@inferred f1(da)) == (2, 1)
    @ballocated $f1($da)
    @test (@ballocated $f1($da)) == 0
    @test dimnum(da, Y, X) == (2, 1)
    @ballocated dimnum($da, Y, X)
    @test (@ballocated dimnum($da, Y, X)) == 0
    @testset "with Dim{X} and symbols" begin
        @test dimnum((Dim{:a}(), Dim{:b}()), :a) == 1
        @test dimnum((Dim{:a}(), Y(), Dim{:b}()), (:b, :a, Y)) == (3, 1, 2)
        dimz = (Dim{:a}(), Y(), Dim{:b}())
        @test (@ballocated dimnum($dimz, :b, :a, Y)) == 0
    end
end

@testset "hasdim" begin
    @test hasdim(da, X()) == true
    @test (@ballocated hasdim($da, X())) == 0
    @test hasdim(da, Ti) == false
    @test (@ballocated hasdim($da, Ti)) == 0
    @test hasdim(dims(da), Y) == true
    @ballocated hasdim(dims($da), Y)
    @test (@ballocated hasdim(dims($da), Y)) == 0
    @test hasdim(dims(da), (X, Y)) == (true, true)
    f1 = (da) -> hasdim(dims(da), (X, Ti, Y, Z))
    @test @inferred f1(da) == (true, false, true, false)
    @test (@ballocated $f1($da)) == 0
    f2 = (da) -> hasdim(dims(da), X, Ti, Y, Z)
    @test (@ballocated $f2($da)) == 0

    @testset "hasdim for Abstract types" begin
        @test hasdim(dims(da), (XDim, YDim)) == (true, true)
        @test hasdim(dims(da), (XDim, XDim)) == (true, false)
        @test hasdim(dims(da), (ZDim, YDim)) == (false, true)
        @test hasdim(dims(da), (ZDim, ZDim)) == (false, false)
        @test (@ballocated hasdim($da, YDim)) == 0
        @test (@ballocated hasdim($da, (ZDim, YDim))) == 0
        @test (@ballocated hasdim($da, ZDim, YDim)) == 0
    end

    @testset "with Dim{X} and symbols" begin
        @test hasdim((Dim{:a}(), Dim{:b}()), (:a, :c)) == (true, false)
        @test hasdim((Dim{:a}(), Dim{:b}()), (:b, :a, :c, :d, :e)) == (true, true, false, false, false)
        @test hasdim((Dim{:a}(), Dim{:b}()), (:b, :a, :c, :d, :e)) == (true, true, false, false, false)
        @ballocated hasdim((Dim{:a}(), Dim{:b}()), (:b, :a, :c, :d, :e))
        @test (@ballocated hasdim((Dim{:a}(), Dim{:b}()), (:b, :a, :c, :d, :e))) == 0
    end
    @test_throws ArgumentError hasdim(nothing, X)
end

@testset "otherdims" begin
    A = DimArray(ones(5, 10, 15), (X, Y, Z));
    @test otherdims(A, X()) == dims(A, (Y, Z))
    @test (@ballocated otherdims($A, X())) == 0
    @test otherdims(A, Y) == dims(A, (X, Z))
    @test otherdims(A, Z) == dims(A, (X, Y))
    @test otherdims(A, (X, Z)) == dims(A, (Y,))
    f1 = A -> otherdims(A, (X, Z))
    @ballocated $f1($A)
    @test (@ballocated $f1($A)) == 0
    f2 = (A) -> otherdims(A, Ti)
    @test f2(A) == dims(A, (X, Y, Z))
    @test (@ballocated $f2($A)) == 0
    @testset "with Dim{X} and symbols" begin
        dimz = (Z(), Dim{:a}(), Y(), Dim{:b}())
        f3 = (dimz) -> otherdims(dimz, (:b, :a, Y))
        @test (@inferred f3(dimz)) == (Z(),)
        @test (@ballocated $f3($dimz)) == 0
        @test otherdims((Dim{:a}(), Dim{:b}(), Ti()), (:a, :c)) == (Dim{:b}(), Ti())
    end
    @test_throws ArgumentError otherdims(nothing, X)
end

@testset "combinedims" begin
    @test combinedims((X(1:10), Y(1:5)), (X(1:10), Z(3:10))) == (X(1:10), Y(1:5), Z(3:10))
    @test_throws DimensionMismatch combinedims((X(1:2), Y(1:5)), (X(1:10), Z(3:10)))
end

@testset "setdims" begin
    A = setdims(da, X(LinRange(150,152,2)))
    @test index(dims(dims(A), X())) == LinRange(150,152,2)
    @test dims(dims(A)) isa Tuple{<:X,<:Y}
    A = setdims(da, Y(10:12), X(LinRange(150,152,2)))
    @test index(dims(dims(A), Y())) == 10:12
    @test dims(dims(A)) isa Tuple{<:X,<:Y}
    @testset "set an empty tuple" begin
        A = setdims(da, ())
        @test dims(A) === dims(da)
    end
end

@testset "swapdims" begin
    @testset "swap type wrappers" begin
        A = swapdims(da, (Z, Dim{:test1}))
        @test dims(A) isa Tuple{<:Z,<:Dim{:test1}}
        @test map(val, dims(A)) == map(val, dims(da))
        @test map(mode, dims(A)) == map(mode, dims(da))
    end
    @testset "swap whole dim instances" begin
        A = swapdims(da, Z(2:2:4), Dim{:test2}(3:5))
        @test dims(A) isa Tuple{<:Z,<:Dim{:test2}}
        @test map(val, dims(A)) == (2:2:4, 3:5)
        @test map(mode, dims(A)) ==
            (Sampled(Ordered(), Regular(2), Points()),
             Sampled(Ordered(), Regular(1), Points()))
    end
    @testset "passing `nothing` keeps the original dim" begin
        A = swapdims(da, (Z(2:2:4), nothing))
        dims(A) isa Tuple{<:Z,<:Y}
        @test map(val, dims(A)) == (2:2:4, val(dims(da, 2)))
        A = swapdims(da, (nothing, Dim{:test3}))
        @test dims(A) isa Tuple{<:X,<:Dim{:test3}}
    end
    @testset "new instances are checked against the array" begin
        @test_throws DimensionMismatch swapdims(da, (Z(2:2:4), Dim{:test4}(3:6)))
    end
end

@testset "slicedims" begin
    @testset "Regular" begin
        @test slicedims(dimz, (1:2, 3)) == slicedims(dimz, 1:2, 3) == slicedims(dimz, (), (1:2, 3)) ==
            ((X(LinRange(143,145,2), Sampled(Ordered(), Regular(2.0), Points()), NoMetadata()),),
             (Y(-36.0, Sampled(Ordered(), Regular(1.0), Points()), NoMetadata()),))
        @test slicedims(dimz, (Z(),), (1:2, 3)) == slicedims(dimz, (Z(),), 1:2, 3) ==
            ((X(LinRange(143,145,2), Sampled(Ordered(), Regular(2.0), Points()), NoMetadata()),),
             (Z(), Y(-36.0, Sampled(Ordered(), Regular(1.0), Points()), NoMetadata()),))
        @test slicedims(dimz, 2:2, :) == slicedims(dimz, (), 2:2, :) ==
            ((X(LinRange(145,145,1), Sampled(Ordered(), Regular(2.0), Points()), NoMetadata()),
              Y(LinRange(-38.0,-36.0, 3), Sampled(Ordered(), Regular(1.0), Points()), NoMetadata())), ())
        @test slicedims((), (1:2, 3)) == slicedims((), (), (1:2, 3)) ==
              slicedims((), 1:2, 3) == slicedims((), (), 1:2, 3) == ((), ())
        @test slicedims(dimz, CartesianIndex(2, 3)) ==
            ((), (X(145.0, Sampled(Ordered(), Regular(2.0), Points()), NoMetadata()),
                         Y(-36.0, Sampled(Ordered(), Regular(1.0), Points()), NoMetadata())),)
    end

    @testset "Irregular" begin
        irreg = DimArray(a, (X([140.0, 142.0]; mode=Sampled(Ordered(), Irregular(140.0, 144.0), Intervals(Start()))),
                             Y([10.0, 20.0, 40.0]; mode=Sampled(Ordered(), Irregular(0.0, 60.0), Intervals(Center()))), ))
        irreg_dimz = dims(irreg)
        @test slicedims(irreg, (1:2, 3)) == slicedims(irreg, 1:2, 3) ==
            ((X([140.0, 142.0], Sampled(Ordered(), Irregular(140.0, 144.0), Intervals(Start())), NoMetadata()),),
                    (Y(40.0, Sampled(Ordered(), Irregular(30.0, 60.0), Intervals(Center())), NoMetadata()),))
        @test slicedims(irreg, (2:2, 1:2)) == slicedims(irreg, 2:2, 1:2) ==
            ((X([142.0], Sampled(Ordered(), Irregular(142.0, 144.0), Intervals(Start())), NoMetadata()),
              Y([10.0, 20.0], Sampled(Ordered(), Irregular(0.0, 30.0), Intervals(Center())), NoMetadata())), ())
        @test slicedims((), (1:2, 3)) == slicedims((), (), (1:2, 3)) == ((), ())
    end

    @testset "Val index" begin
        da = DimArray(a, (X(Val((143, 145))), Y(Val((:x, :y, :z)))))
        dimz = dims(da)
        @test slicedims(dimz, (1:2, 3)) == slicedims(dimz, (), (1:2, 3)) ==
            ((X(Val((143,145)), Categorical(), NoMetadata()),),
             (Y(:z, Categorical(), NoMetadata()),))
        @test slicedims(dimz, (2:2, :)) ==
            ((X(Val((145,)), Categorical(), NoMetadata()),
              Y(Val((:x, :y, :z)), Categorical(), NoMetadata())), ())
    end

    @testset "NoIndex" begin
        da = DimArray(a, (X(), Y()))
        dimz = dims(da)
        @test slicedims(dimz, (1:2, 3)) == ((X(1:2, NoIndex()),), (Y(3, NoIndex()),))
        @test slicedims(dimz, (2:2, :)) == ((X(2:2, NoIndex()), Y(Base.OneTo(3), NoIndex())), ())
    end
    @testset "No slicing" begin
        da = DimArray(a, (X([143, 145]), Y([:x, :y, :z])))
        dimz = dims(da)
        @test slicedims(dimz, ()) == (dimz, ())
    end
end

@testset "reducedims" begin
    @test _reducedims((X(3:4; mode=Sampled(Ordered(), Regular(1), Points())),
                      Y(1:5; mode=Sampled(Ordered(), Regular(1), Points()))), (X, Y)) ==
                     (X([4], Sampled(Ordered(), Regular(2), Points()), NoMetadata()),
                      Y([3], Sampled(Ordered(), Regular(5), Points()), NoMetadata()))
    @test _reducedims((X(3:4; mode=Sampled(Ordered(), Regular(1), Intervals(Start()))),
                      Y(1:5; mode=Sampled(Ordered(), Regular(1), Intervals(End())))), (X, Y)) ==
        (X([3], Sampled(Ordered(), Regular(2), Intervals(Start())), NoMetadata()),
         Y([5], Sampled(Ordered(), Regular(5), Intervals(End())), NoMetadata()))

    @test _reducedims((X(3:4; mode=Sampled(Ordered(), Irregular(2.5, 4.5), Intervals(Center()))),
                      Y(1:5; mode=Sampled(Ordered(), Irregular(0.5, 5.5), Intervals(Center())))), (X, Y))[1] ==
                     (X([4], Sampled(Ordered(), Irregular(2.5, 4.5), Intervals(Center())), NoMetadata()),
                      Y([3], Sampled(Ordered(), Irregular(0.5, 5.5), Intervals(Center())), NoMetadata()))[1]
    @test _reducedims((X(3:4; mode=Sampled(Ordered(), Irregular(3, 5), Intervals(Start()))),
                      Y(1:5; mode=Sampled(Ordered(), Irregular(0, 5), Intervals(End()  )))), (X, Y))[1] ==
                     (X([3], Sampled(Ordered(), Irregular(3, 5), Intervals(Start())), NoMetadata()),
                      Y([5], Sampled(Ordered(), Irregular(0, 5), Intervals(End()  )), NoMetadata()))[1]

    @test _reducedims((X(3:4; mode=Sampled(Ordered(), Irregular(), Points())),
                      Y(1:5; mode=Sampled(Ordered(), Irregular(), Points()))), (X, Y)) ==
        (X([4], Sampled(Ordered(), Irregular(), Points()), NoMetadata()),
         Y([3], Sampled(Ordered(), Irregular(), Points()), NoMetadata()))
    @test _reducedims((X(3:4; mode=Sampled(Ordered(), Regular(1), Points())),
                      Y(1:5; mode=Sampled(Ordered(), Regular(1), Points()))), (X, Y)) ==
                     (X([4], Sampled(Ordered(), Regular(2), Points()), NoMetadata()),
                      Y([3], Sampled(Ordered(), Regular(5), Points()), NoMetadata()))

    @test _reducedims((X([:a,:b]; mode=Categorical()),
                      Y(["1","2","3","4","5"]; mode=Categorical())), (X, Y)) ==
                     (X([:combined]; mode=Categorical()),
                      Y(["combined"]; mode=Categorical()))

    @test _reducedims((X(Base.OneTo(10); mode=NoIndex()),
                       Y(Base.OneTo(10); mode=NoIndex())), (X(), Y())) ==
        (X(Base.OneTo(1); mode=NoIndex()), Y(Base.OneTo(1); mode=NoIndex()))

    @testset "Special case CompoundPeriod" begin
        step_ = Dates.CompoundPeriod([Month(1), Day(3)])
        timespan = [DateTime(2001, 1), DateTime(2001, 1, 3)]
        teststep = Dates.CompoundPeriod([Month(2), Day(6)])
        testdim = Ti(timespan[2:2], Sampled(Ordered(), Regular(teststep), Points()), NoMetadata())
        reduceddim = _reducedims((Ti(timespan; mode=Sampled(Ordered(), Regular(step_), Points())),), (Ti,))[1]
        @test typeof(testdim) == typeof(reduceddim)
        @test val(testdim) == val(reduceddim)
        @test step(testdim) == step(reduceddim)
    end
end

@testset "dimstride" begin
    dimz = (X(), Y(), Dim{:test}())
    da = DimArray(ones(3, 2, 3), dimz; name=:data)
    @test dimstride(da, X()) == 1
    @test dimstride(da, Y()) == 3
    @test dimstride(da, Dim{:test}()) == 6
    @test_throws ArgumentError dimstride(nothing, X())
end


text:
########################
### TESTING ROUTINES ###
########################

# TEST THE EXCHANGE ALGORITHM ON RANDOM DATA
#
# This function simulates random uncorrelated Gaussian data to test the exchange algorithm.
#
# Arguments:
# -- n is the desired number of samples.
# -- p is the desired number of predictors.
# -- r is the desired number of true nonzero components, or the size of the true model.
# -- extra is the amount of extra model sizes to compute beyond r. Thus, the length of the computed
#    regularization path is r + extra.
#
# Optional Arguments:
# -- max_iter caps the number of outer iterations in the exchange algorithm. Defaults to 100.
# -- quiet is a Bool to control algorithm output to the console. Defaults to false (print output).
# -- window is the number of predictors to consider for swapping. Defaults to r + extra, which tells the algorithm
#    to consider all active predictors for swapping.
# -- seed sets the random seed for the data generation. Defaults to 2015.
# -- noise indicates the magnitude of the random noise added to the response y. Defaults to 0.5, which indicates
#    noise drawn from the N(0,0.25) distribution
#
# coded by Kevin L. Keys (2015)
# klkeys@g.ucla.edu
function test_exleastsq(n::Int, p::Int, r::Int, extra::Int; max_iter::Int = 100, quiet::Bool = false, window::Int = r+extra, seed::Int = 2015, noise::Float64 = 0.5)
	tic()

	# testing data
	x         = randn(n,p)
	y         = randn(n)
	b         = zeros(p)
	perm      = collect(1:p)
	inner     = Dict{Int,DenseArray{Float64,1}}()

	# declare all temporary arrays
	res        = zeros(n)	# Y - Xbeta
	df         = zeros(p)	# X'(Y - Xbeta)
	tempn      = zeros(n)	# temporary array of length n
	tempn2     = zeros(n)	# temporary array of length n
	tempp      = zeros(p)	# temporary array of length p
	dotprods   = zeros(p)	# hold in memory the dot products for current index

	# precompute sum of squares for each column of x
	const nrmsq = vec(sumabs2(x,1))

	# declare return array
	betas = zeros(p,r+extra)

	# use a random permutation of indices for true model
	true_model = sample(perm, r)
	b[true_model] = randn(r)
	b_true = copy(b)
	BLAS.gemv!('N', 1.0, x, b, 1.0, y)	# y = xb + noise
	fill!(b, 0.0)
	iter = 0

	# print time spent on parameter initialization
	quiet || println("Parameter initialization took ", toq(), " seconds.")

	# reset timer and compute path
	tic()
	@time begin
		for k = 1:(r+extra)
			b = exchange_leastsq!(b, x, y, perm, k, inner=inner, max_iter=max_iter, quiet=quiet, nrmsq=nrmsq, res=res, df=df, tempn=tempn, tempn2=tempn2, tempp=tempp, dotprods=dotprods, window=min(k,window))
			rss = sumabs2(res)

			# print results from this step of regularization path
			if !quiet
				println("#pred = ", k, ", rss = ", rss, ", #truepos = ", countnz(b[true_model]))
				print("\t")
				print("\n")
			end
			betas[:,k] = b
			update_col!(betas, b, k, n=p, p=r+extra, a=1.0)
		end
	end

	# print time spent on computation of regularization path
	quiet || println("Path from 1 to ", r+extra, " took ", toq(), " seconds to compute.")

	# return sparsified betas
	return b_true, sparse(betas)
end

# TEST THE EXCHANGE ALGORITHM FOR LEAST SQUARES REGRESSION
# This subroutine runs the IHT on a dataset of chromosome 1 data from the WTCCC.
function test_exchangeleastsq(x_path::String = "/Users/kkeys/Downloads/wtccc-n2000p32307.txt", y_path::String = "/Users/kkeys/Downloads/withnoiselevelsd0_1/Y.100.1", b_path::String = "/Users/kkeys/Downloads/withnoiselevelsd0_1/causal.100.1", r::Int = 100, tol::Float64 = 1e-6, max_iter::Int = 100, quiet::Bool = false, extra::Int = 100, window::Int = r)

    # repeat input options
    println("This function will test the exchange algorithm.")
    println("Given options:")
    println("\tPath to X  = ", x_path)
    println("\tPath to Y  = ", y_path)
    println("\tPath to B  = ", b_path)
    println("\tModel size = ", r)
    println("\ttolerance  = ", tol, "\n\tmax_iter   = ", max_iter)

    # precompile @time macro for later use
    println("\nCompiling @time macro...")
    @time 1+1;

    # now load data
    println("\nLoading data...")
    tic()

    # load design matrix
    x = readdlm(x_path)
    const (n,p) = size(x)

    # load response vector
    y = readdlm(y_path)
    y = vec(y)  # need y to be 1D

    # load model
    B = readdlm(b_path)

    # need indices for true model
    # add 1 since Julia uses 1-indexing
    bidx = convert(Array{Int,1}, B[:,2]) + 1

    # need components of true model
    be = convert(Array{Float64,1}, B[:,3])

    # how long did it take to load files?
	file_time = toq()
    quiet || println("Files took ", file_time, " seconds to load.")

	# load parameters
	println("Allocating temporary arrays...")
	tic()

	# declare all temporary arrays
    b          = zeros(p)	# parameter vector
	res        = zeros(n)	# Y - Xbeta
	df         = zeros(p)	# X'(Y - Xbeta)
	tempn      = zeros(n)	# temporary array of length n
	tempn2     = zeros(n)	# temporary array of length n
	tempp      = zeros(p)	# temporary array of length p
	dotprods   = zeros(p)	# hold in memory the dot products for current index
	perm       = collect(1:p) # hold permutation vector of b

	# declare associative array for storing inner products
	inner = Dict{Int,DenseArray{Float64,1}}()

	# how long did the allocation take?
	alloc_time = toq()
	quiet || println("Temporary arrays allocated in ", alloc_time, " seconds.")

	# precompute sum of squares for each column of x
	quiet || println("Precomputing squared Euclidean norms of matrix columns...")
	tic()
	const nrmsq = vec(sumabs2(x,1))
	norm_time = toq()
	quiet || println("Column norms took ", norm_time, " seconds to compute")

	# declare any return values
	iter = 0

	# print time spent on parameter initialization
	parameter_time = norm_time + alloc_time
	quiet || println("Parameter initialization took ", parameter_time, " seconds.")

    # run exchange
    println("\nRunning exchange algorithm...")

	# reset timer and compute path
	tic()
	@time begin
		for i = 1:(r+extra)
			b = exchange_leastsq!(b, x, y, perm, i, inner=inner, max_iter=max_iter, quiet=quiet, nrmsq=nrmsq, res=res, df=df, tempn=tempn, tempn2=tempn2, tempp=tempp, dotprods=dotprods, window=min(i,window))
			rss = sumabs2(res)

			# print results from this step of regularization path
			quiet || println("#pred = ", i, ", rss = ", rss, ", #truepos = ", countnz(b[bidx]))
		end
	end

	# print time spent on computation of regularization path
	quiet || println("Path from 1 to ", r+extra, " took ", toq(), " seconds to compute.")

    # recover vector from output
    bk = b[bidx]

    # evaluate model
    println("\nTrue positives: ", countnz(bk), "/", length(be), ".")
    println("Distances to true model:")
    println("Chebyshev (L-Inf) = ", chebyshev(bk, be), ".")
    println("Euclidean (L-2)   = ", euclidean(bk, be), ".")
end

# TEST THE EXCHANGE ALGORITHM FOR LEAST SQUARES REGRESSION OVER PLINK BINARY DATA
# This subroutine runs the exchange algorithm on a GWAS dataset from the WTCCC.
function test_exchangeleastsq_plink(x_path::String = "/Users/kkeys/Downloads/wtccc_full/wtccc-n2k_chr1_clean.bed", xt_path::String = "/Users/kkeys/Downloads/wtccc_full/wtccc-n2k_chr1_clean_t.bed", y_path::String = "/Users/kkeys/Downloads/withnoiselevelsd0_1/Y.100.1", b_path::String = "/Users/kkeys/Downloads/withnoiselevelsd0_1/causal.100.1", r::Int = 10, tol::Float64 = 1e-4, max_iter::Int = 100, quiet::Bool = false, extra::Int = 10, window::Int = 10)

    # repeat input options
	quiet || begin
		println("This function will test the IHT algorithm.")
		println("Given options:")
		println("\tPath to X  = ", x_path)
		println("\tPath to X' = ", xt_path)
		println("\tPath to Y  = ", y_path)
		println("\tPath to B  = ", b_path)
		println("\tModel size = ", r)
		println("\ttolerance  = ", tol, "\n\tmax_iter   = ", max_iter)
	end

    # precompile @time macro for later use
    quiet || println("\nCompiling @time macro...")
    @time 1+1;

    # now load data
    quiet || println("\nLoading data...")
    tic()

    # load design matrix
	x = BEDFile(x_path, xt_path)
	const (n,p) = size(x)

    # load response vector
    y = readdlm(y_path)
	y = vec(y) # need y to be 1D
	const y = convert(SharedArray{Float64,1}, y)

	const means   = PLINK.mean(Float64, x, shared=true)
	const invstds = PLINK.invstd(x, means, shared=true)

    # load model
    B = readdlm(b_path)

    # need indices for true model
    # add 1 since Julia uses 1-indexing
	const bidx = convert(Array{Int,1}, B[:,2]) + 1

    # need components of true model
	const be = convert(Array{Float64,1}, B[:,3])

    # how long did it take to load files?
	file_time = toq()
    quiet || println("\nFiles took ", file_time, " seconds to load.")

	# discard B to recover some memory
	quiet || println("Object B is no longer needed, recovering memory...")
	tic()
	B = false
	gc()
	B_time = toq()
	quiet || println("Discarding B took ", B_time, " seconds.")

	# load parameters
	quiet || println("\nAllocating temporary arrays...")
	tic()

	# declare all temporary arrays
	b        = SharedArray(Float64, p, init = S -> S[localindexes(S)] = 0.0)		# parameter vector
	Xb       = SharedArray(Float64, n, init = S -> S[localindexes(S)] = 0.0)		# Xbeta
	res      = SharedArray(Float64, n, init = S -> S[localindexes(S)] = 0.0)		# Y - Xbeta
	df       = SharedArray(Float64, p, init = S -> S[localindexes(S)] = 0.0)		# X'(Y - Xbeta)
	tempn    = SharedArray(Float64, n, init = S -> S[localindexes(S)] = 0.0)		# temporary array of length n
	tempn2   = SharedArray(Float64, n, init = S -> S[localindexes(S)] = 0.0)		# temporary array of length n
	tempp    = SharedArray(Float64, p, init = S -> S[localindexes(S)] = 0.0)		# temporary array of length p
	dotprods = SharedArray(Float64, p, init = S -> S[localindexes(S)] = 0.0)		# hold in memory the dot products for current index
	indices  = falses(p)

	# permutation vector of b has a more complicated initialization
	perm = SharedArray(Int, p, init = S -> S[localindexes(S)] = localindexes(S))

	# declare associative array for storing inner products
	inner = Dict{Int,SharedArray{Float64,1}}()

	# how long did the allocation take?
	alloc_time = toq()
	quiet || println("Temporary arrays allocated in ", alloc_time, " seconds.")

	# precompute sum of squares for each column of x
	quiet || println("\nPrecomputing squared Euclidean norms of matrix columns...")
	tic()
	const nrmsq = sumsq(Float64, x, shared=true, means=means, invstds=invstds)
	norm_time = toq()
	quiet || println("Column norms took ", norm_time, " seconds to compute")

	# print time spent on parameter initialization
	parameter_time = norm_time + alloc_time + file_time + B_time
	quiet || println("Parameter initialization took ", parameter_time, " seconds.")

    # run exchange
    quiet || println("\nRunning exchange algorithm...")

	# reset timer and compute path
	tic()
	@time begin
		for i = 2:(r+extra)
			exchange_leastsq!(b, x, y, perm, i, p=p, n=n, inner=inner, max_iter=max_iter, quiet=quiet, nrmsq=nrmsq, res=res, df=df, tempn=tempn, tempn2=tempn2, tempp=tempp, dotprods=dotprods, window=min(i,window), Xb=Xb, means=means, invstds=invstds, indices=indices)
			rss = sumabs2(res)

			# print results from this step of regularization path
			quiet || println("#pred = ", i, ", rss = ", rss, ", #truepos = ", countnz(b[bidx]))
		end
	end

	# print time spent on computation of regularization path
	quiet || println("Path from 1 to ", r+extra, " took ", toq(), " seconds to compute.")

    # recover vector from output
    bk = b[bidx]

    # evaluate model
	quiet || begin
		println("\nTrue positives: ", countnz(bk), "/", length(be), ".")
		println("Distances to true model:")
		println("Chebyshev (L-Inf) = ", chebyshev(bk, be), ".")
		println("Euclidean (L-2)   = ", euclidean(bk, be), ".")
	end

	return bk, bidx
end

test_exchangeleastsq_plink(r::Int, extra::Int, window::Int) = test_exchangeleastsq_plink("/Users/kkeys/Downloads/wtccc_full/wtccc-n2k_chr1_clean.bed", "/Users/kkeys/Downloads/wtccc_full/wtccc-n2k_chr1_clean_t.bed", "/Users/kkeys/Downloads/withnoiselevelsd0_1/Y.100.1", "/Users/kkeys/Downloads/withnoiselevelsd0_1/causal.100.1", r, 1e-4, 100, false, extra, window)

# CROSSVALIDATION ROUTINE FOR EXCHANGE ALGORITHM
function test_cv_exlstsq(kend,numfolds; tol::Float64 = 1e-4, max_iter::Int = 1000, quiet::Bool=true, compute_model::Bool = false)

	# notify of number of folds
	quiet || println("Testing ", numfolds, "-fold crossvalidation from k = 1 to k = ", kend, "...")

	# start timing
	tic()

	# file paths for data, response, model
	x_path = "/Users/kkeys/Downloads/wtccc-n2000p32307.txt"
	y_path = "/Users/kkeys/Downloads/withnoiselevelsd0_1/Y.100.1"
	b_path = "/Users/kkeys/Downloads/withnoiselevelsd0_1/causal.100.1"

	# load design matrix
	x = readdlm(x_path)
	const (n,p) = size(x)
	const x = convert(SharedArray, x)

	# load response vector
	y = readdlm(y_path)
	y = vec(y)	# need y to be 1D
	const y = convert(SharedArray,y)

	# load model
	B = readdlm(b_path)

	# need indices for true model
	# add 1 since Julia uses 1-indexing
	const bidx = convert(Array{Int,1}, B[:,2]) + 1

	# need components of true model
	const be = convert(Array{Float64,1}, B[:,3])

	parameter_time = toq()

	# inform console that the data are loaded
	quiet || println("Data loaded in ", parameter_time, " seconds.")

	# obtain a partition of the sample set
	myfolds = cv_get_folds(y,numfolds)

	# inform console that the folds are computed
	quiet || println("Folds computed.")

	# run cv_iht to crossvalidate dataset
	# be careful! mses will be a tuple if compute_model=true
	tic()
	if compute_model
		mses, b = cv_exlstsq(x,y,kend,numfolds, p=p, tol=tol, max_iter=max_iter, quiet=quiet, folds=myfolds, compute_model=compute_model)
	else
		mses = cv_exlstsq(x,y,kend,numfolds, p=p, tol=tol, max_iter=max_iter, quiet=quiet, folds=myfolds, compute_model=compute_model)
	end

	cv_time = toq()
	quiet || println("CV computed in ", cv_time, " seconds.")

	compute_model && return mses, b
	return mses
end

function compare_exlstsq(x::Matrix, kend, numfolds; tol::Float64 = 1e-4, max_iter::Int = 1000, quiet::Bool=true, compute_model::Bool = false, y_path::String = "/Users/kkeys/Desktop/Thesis/Exchange/Y_1_same_1.0_2.0.txt", b_path::String = "/Users/kkeys/Desktop/Thesis/Exchange/b_1_same_1.0_2.0.txt", bidx_path::String = "/Users/kkeys/Desktop/Thesis/Exchange/bidx_1_same_1.0_2.0.txt")

	# notify of number of folds
	quiet || println("Testing ", numfolds, "-fold crossvalidation from k = 1 to k = ", kend, "...")

	# start timing
	tic()

	# dimensions of design matrix
	const n,p = size(x)
	x = convert(SharedArray{Float64,2}, x)

	# load response vector
	y = readdlm(y_path)
	y = vec(y)	# need y to be 1D
#	y = convert(Array{Float64,1}, y)
	y = convert(SharedArray{Float64,1}, y)

	# load model
	be = readdlm(b_path)
	be = vec(be)
	const be = convert(Array{Float64,1}, be)

	# need indices for true model
	bidx = readdlm(bidx_path)
	bidx = vec(bidx)
	const bidx = convert(Array{Int,1}, bidx)

	parameter_time = toq()

	# inform console that the data are loaded
	quiet || println("Data loaded in ", parameter_time, " seconds.")

	# start with exchange algorithm
	quiet || println("Running crossvalidation with exchange algorithm...")

	# obtain a partition of the sample set
	myfolds = cv_get_folds(y,numfolds)

	# inform console that the folds are computed
	quiet || println("Folds computed.")

	# run cv_iht to crossvalidate dataset
	# be careful! mses will be a tuple if compute_model=true
	tic()
#	mses = cv_exlstsq(x,y,kend,numfolds, p=p, tol=tol, max_iter=max_iter, quiet=quiet, folds=myfolds, compute_model=compute_model)
	if compute_model
		mses, b, ex_model = cv_exlstsq(x,y,kend,numfolds, p=p, tol=tol, max_iter=max_iter, quiet=quiet, folds=myfolds, compute_model=compute_model)
	else
		mses = cv_exlstsq(x,y,kend,numfolds, p=p, tol=tol, max_iter=max_iter, quiet=quiet, folds=myfolds, compute_model=compute_model)
	end

	cv_time = toq()
	quiet || println("CV computed in ", cv_time, " seconds.")

	# now compute LASSO model
	quiet || println("\nComputing LASSO model...")

	# translate folds

	# run glmnet
	tic()
	cv = glmnetcv(sdata(x), sdata(y), dfmax=kend, nlambda=kend, nfolds=5, parallel=true, folds=my_folds, tol=tol)

	# extract best beta from lasso
	lasso_btemp = cv.path.betas[:,indmin(cv.meanloss)]

	# inform console that lasso computation is complete
	lasso_time = toq()
	quiet || println("LASSO model crossvalidated in ", lasso_time, " seconds.")

	if compute_model
#		lasso_bidx = lasso_btemp .!= 0.0
		lasso_bidx = find( function f(x) x .!= 0.0; end, lasso_btemp)
		x_lasso = x[:,lasso_bidx]
		lasso_b = (x_lasso'*x_lasso) \ (x_lasso'*y)
	end

	# return values based on compute_model
	compute_model && return mses, cv.meanloss, b, lasso_b, be, ex_model, lasso_bidx, bidx
	return mses, cv.meanloss
end

function compare_exlstsq(x::Matrix, y::Vector, be::Vector, bidx::Vector, kend::Int, numfolds::Int; tol::Float64 = 1e-4, max_iter::Int = 10000, quiet::Bool=true, compute_model::Bool = false, n::Int = length(y), p::Int = length(be), window::Int = 20, nrmsq::DenseArray{Float64,1} = vec(sumsq(x,1)))

	isequal((n,p),size(x)) || throw(DimensionMismatch("(n,p) = ($n,$p) != ($(size(x)))"))
	0 <= window <= p       || throw(error("For current X, must have window in (0,$p)"))
	0 <= kend <= p         || throw(error("Path legnth must be positive and cannot exceed the number of predictors"))

	# notify of number of folds
	quiet || println("Testing ", numfolds, "-fold crossvalidation from k = 1 to k = ", kend, "...")

	# start with exchange algorithm
	quiet || println("Running crossvalidation with exchange algorithm...")

	# obtain a partition of the sample set
	myfolds = cv_get_folds(y,numfolds)

	# inform console that the folds are computed
	quiet || println("Folds computed.")

	# run cv_iht to crossvalidate dataset
	tic()
	if compute_model
		mses, b, ex_model = cv_exlstsq(x,y,kend,numfolds, p=p, tol=tol, max_iter=max_iter, quiet=quiet, folds=myfolds, compute_model=compute_model, window=window)
	else
		mses = cv_exlstsq(x,y,kend,numfolds, p=p, tol=tol, max_iter=max_iter, quiet=quiet, folds=myfolds, compute_model=compute_model, window=window)
	end
	cv_time = toq()
	quiet || println("CV computed in ", cv_time, " seconds.")

	# now compute LASSO model
	quiet || println("\nComputing LASSO model...")

	# run glmnet
	tic()
#	cv = glmnetcv(x,y, dfmax=kend, nlambda=kend, nfolds=numfolds, parallel=true, folds=my_folds, tol=tol)

	# extract best beta from lasso
#	lasso_btemp = cv.path.betas[:,indmin(cv.meanloss)]
	lasso_btemp = zeros(p)
	lasso_btemp[1] = 1.0

	# inform console that lasso computation is complete
	lasso_time = toq()
	quiet || println("LASSO model crossvalidated in ", lasso_time, " seconds.")

	if compute_model
#		lasso_bidx = lasso_btemp .!= 0.0
		lasso_bidx = find( function f(x) x .!= 0.0; end, lasso_btemp)
		x_lasso = x[:,lasso_bidx]
		lasso_b = (x_lasso'*x_lasso) \ (x_lasso'*y)
	end

	# return values based on compute_model
	compute_model && return mses, zeros(kend), b, lasso_b, be, ex_model, lasso_bidx, bidx, lasso_time, cv_time
	return mses, cv.meanloss, cv_time, lasso_time
end

ext: jl


text:
using Markdown
using GithubMarkdown
using HTMLSanitizer
using Highlights
using Downloads

function build_git_docs(packagespec, buildpath, uri; src_prefix="", href_prefix="")
    pkgname = packagespec.name
    return mktempdir() do dir
        return cd(dir) do
            run(`git clone --depth=1 $(uri) $(pkgname)`)
            docsproject = joinpath(dir, pkgname)
            return cd(docsproject) do
                return build_local_docs(packagespec, buildpath, nothing, docsproject, gitdirdocs = true; src_prefix=src_prefix, href_prefix=href_prefix)
            end
        end
    end
end

function build_hosted_docs(packagespec, buildpath, uri)
    pkgname = packagespec.name
    # js redirect
    open(joinpath(buildpath, "index.html"), "w") do io
        println(io,
            """
            <!DOCTYPE html>

                    <meta http-equiv="Content-Type" content="0; $(uri)" />

                    Redirecting to $(uri).

            """
        )
    end

    # check for meta redirects and follow them
    uri = maybe_redirect(uri)

    # download search index
    try
        Downloads.download(string(uri, "/search_index.js"), joinpath(buildpath, "search_index.js"))
    catch err
        @error("Search index download failed for `$(uri)`.", exception = err)
    end
    return Dict(
        "doctype" => :hosted,
        "success" => true,
        "hosted_uri" => uri
    )
end

using Gumbo, AbstractTrees
function maybe_redirect(uri)
    docspage = Downloads.download(uri)
    docspage = read(docspage, String)

    try
        doc = Gumbo.parsehtml(docspage)
        for el in PreOrderDFS(doc.root)
            if el isa HTMLElement && Gumbo.tag(el) == :meta && "content" in keys(Gumbo.attrs(el))
                content = getattr(el, "content")
                m = match(r"\d+;\s*url\=(.*)$", content)

                if m !== nothing
                    @info("Found redirect from `$(uri)` to `$(m[1])`.")
                    return m[1]
                end
            end
        end
    catch err
        @error("Failed to determine if there's a redirect in `$(uri)`.", err=err)
    end
    return uri
end

function build_local_docs(packagespec, buildpath, uri, pkgroot = nothing; gitdirdocs = false, src_prefix="", href_prefix="")
    uri = something(uri, "docs")
    mktempdir() do envdir
        pkgname = packagespec.name
        installable = try_install_package(packagespec, envdir)

        documenter_errored = false

        if !installable
            @error("$(pkgname) is not installable. Stopped building docs.")
            return Dict(
                "installable" => false,
                "success" => false
            )
        end
        if pkgroot === nothing
            pkgfile = Base.find_package(pkgname)
            pkgroot = normpath(joinpath(pkgfile, "..", ".."))
        end
        mod = try_use_package(packagespec)

        # actual Documenter docs
        try
            for docdir in joinpath.(pkgroot, unique([uri, "docs", "doc"]))
                if isdir(docdir)
                    @info("Building vendored Documenter.jl documentation at $(docdir).")
                    output = build_documenter(packagespec, docdir)
                    @info("Documentation built at $(output).")
                    if output !== nothing
                        @info("Copying build documentation from $(output) to $(buildpath)")
                        cp(output, buildpath, force = true)
                        return Dict(
                            "doctype" => gitdirdocs ? "git-repo" : "documenter",
                            "documenter_errored" => documenter_errored,
                            "installable" => true,
                            "using_failed" => mod == nothing,
                            "success" => true
                        )
                    else
                        @error("Errored while trying to generate Documenter docs at $(docdir).")
                        documenter_errored = true
                    end
                end
            end
        catch err
            @error("Errored while trying to generate Documenter of $(pkgname)", exception = (err, catch_backtrace()))

            documenter_errored = true
        end

        documenter_errored && @warn("Building vendored Documenter docs failed.")
        @info("Building fallback documentation.")

        # fallback docs (readme & docstrings)
        return mktempdir() do docsdir
            output = build_readme_docs(pkgname, pkgroot, docsdir, mod, src_prefix, href_prefix)
            if output !== nothing
                cp(output, buildpath, force = true)
                return Dict(
                    "doctype" => "fallback_autodocs",
                    "documenter_errored" => documenter_errored,
                    "installable" => true,
                    "using_failed" => mod == nothing,
                    "success" => true
                )
            end
            return Dict(
                "doctype" => "fallback_autodocs",
                "documenter_errored" => documenter_errored,
                "installable" => true,
                "using_failed" => mod == nothing,
                "success" => false
            )
        end
    end
end

function build_legacy_documenter(packagespec, docdir)
    open(joinpath(docdir, "Project.toml"), "w") do io
        println(io, """
            [deps]
            Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"

            [compat]
            Documenter = "~0.20"
        """)
    end
    build_documenter(packagespec, docdir)
end

function build_documenter(packagespec, docdir)
    pkgdir = normpath(joinpath(docdir, ".."))
    cd(pkgdir) do
        docsproject = joinpath(docdir, "Project.toml")
        docsmanifest = joinpath(docdir, "Manifest.toml")
        if !isfile(docsproject)
            return build_legacy_documenter(packagespec, docdir)
        end

        # fix permissions to allow us to add the main pacakge to the docs project
        chmod(docsproject, 0o660)
        isfile(docsmanifest) && chmod(docsmanifest, 0o660)

        chmod(joinpath(pkgdir, "Project.toml"), 0o660)
        isfile(joinpath(pkgdir, "Manifest.toml")) && chmod(joinpath(pkgdir, "Manifest.toml"), 0o660)

        rundcocumenter = joinpath(@__DIR__, "rundocumenter.jl")

        makefile = joinpath(docdir, "make.jl")
        if !isfile(makefile)
            @warn("No `make.jl` file found in `$(docdir)`. Trying the first `*.jl` file instead.")
            jlfiles = filter(file -> endswith(file, ".jl"), readdir(docdir))
            if isempty(jlfiles)
                @error("No `*.jl` files found in `$(docdir)`. Aborting.")
                return nothing
            end
            makefile = joinpath(docdir, jlfiles[1])
            @info("Using $(makefile) to generate Documenter docs.")
        end
        _, builddir = fix_makefile(makefile)

        cmd = ```
            $(first(Base.julia_cmd()))
                --project="$(docdir)"
                --compiled-modules=no
                -O0
                $(rundcocumenter)
                $(pkgdir)
                $(makefile)
            ```

        try
            run(cmd)

            return builddir
        catch err
            @error("Failed to evaluate specified make.jl-file.", exception=err)
            return nothing
        end
    end
end

function build_readme_docs(pkgname, pkgroot, docsdir, mod, src_prefix, href_prefix)
    @info("Generating readme-only fallback docs.")

    if pkgroot === nothing || !ispath(pkgroot)
        @error("Julia could not find the package directory. Aborting.")
        return
    end

    pkgloads = mod !== nothing

    readme = find_readme(pkgroot)
    doc_src = joinpath(docsdir, "src")
    mkpath(doc_src)
    index = joinpath(doc_src, "index.md")

    render_html(readme, index, src_prefix, href_prefix; documenter = true)

    if !isfile(index)
        open(index, "w") do io
            println(io, """
            # $pkgname
            """)
        end
    end

    pages = ["Readme" => "index.md"]
    modules = :(Module[Module()])

    if pkgloads
        @info("Deploying `autodocs`.")
        add_autodocs(doc_src, mod)
        push!(pages, "Docstrings" => "autodocs.md")
        modules = :(Module[$mod])
    end

    @eval Module() begin
        using Pkg
        Pkg.add("Documenter")
        using Documenter
        makedocs(
            format = Documenter.HTML(),
            sitename = "$($pkgname).jl",
            modules = $(modules),
            root = $(docsdir),
            pages = $(pages)
        )
    end

    build_dir = joinpath(docsdir, "build")
    if ispath(build_dir)
        return build_dir
    end
    return nothing
end

function find_readme(pkgroot)
    allfiles = readdir(pkgroot)
    # look for readme.md/readme first
    for file in allfiles
        if lowercase(file) in ("readme.md", "readme")
            readme = joinpath(pkgroot, file)
            if isfile(readme)
                return readme
            end
        end
    end
    # if we didn't find one of those, take the first file that contains "readme"
    for file in allfiles
        if occursin("readme", lowercase(file))
            readme = joinpath(pkgroot, file)
            if isfile(readme)
                return readme
            end
        end
    end
end

function add_autodocs(docsdir, mod)
    open(joinpath(docsdir, "autodocs.md"), "w") do io
        println(io, """
        ```@autodocs
        Modules = [$(Symbol(mod))]
        ```
        """)
    end
end

function monkeypatchdocsearch(packagespec, buildpath)
    uuid = packagespec.uuid
    name = packagespec.name
    if !(get(ENV, "DISABLE_CENTRALIZED_SEARCH", false) in ("true", "1", 1))
        searchjs = joinpath(buildpath, "assets", "search.js")
        if isfile(searchjs)
            @info "monkey patching search.js for $(name)"
            rm(searchjs, force=true)
            template = String(read(joinpath(@__DIR__, "search.js.template")))
            template = replace(template, "{{{UUID}}}" => String(uuid))
            open(searchjs, "w") do io
                print(io, template)
            end
        end
    end
end

function copy_package_source(packagespec, buildpath)
    outpath = joinpath(buildpath, "_packagesource")
    try
        mktempdir() do envdir
            if packagespec.name == "julia"
                Downloads.download("https://github.com/JuliaLang/julia/releases/download/v$(packagespec.version)/julia-$(packagespec.version).tar.gz", joinpath(envdir, string(packagespec.uuid)*".tar.gz"))
                run(`tar -xzf $(joinpath(envdir, string(packagespec.uuid)*".tar.gz"))  -C $envdir`)
                docs_path = joinpath(envdir, packagespec.name*"-"*string(packagespec.version), "doc", "_build", "html", "en")
                src_path = joinpath(envdir, packagespec.name*"-"*string(packagespec.version))
                cp(docs_path , joinpath(buildpath, "_docsbuild"), force = true)
                cp(src_path, outpath, force = true)
                return outpath
            end
            pkgname = packagespec.name
            installable = try_install_package(packagespec, envdir)

            if !installable
                @error("Package not installable. Can't get source code.")
                return
            end
            pkgfile = Base.find_package(pkgname)
            pkgroot = normpath(joinpath(pkgfile, "..", ".."))

            @info("Copying source code for $(pkgname).")
            if isdir(pkgroot)
                cp(pkgroot, outpath; force=true)
            end
            @info("Done copying source code for $(pkgname).")
            return outpath
        end
    catch err
        @error("Error trying to copy package source.", exception=err)
    end
end

function render_html(input, output, src_prefix="", href_prefix=""; documenter = false)
    if input === nothing
        @error("Package doesn't seem to have a readme. ")

        return
    end

    io = IOBuffer()
    GithubMarkdown.rendergfm(io, input; documenter = false)

    allow_class_on_code = deepcopy(HTMLSanitizer.WHITELIST)
    allow_class_on_code[:attributes]["code"] = ["class"]

    out = sanitize(String(take!(io)), prettyprint = false, whitelist = allow_class_on_code)

    out = postprocess_html_readme(out; src_prefix = src_prefix, href_prefix = href_prefix)
    out = out.children[2]

    print(io, out)
    out = String(take!(io))

    out = replace(out, r"^" => "")
    out = replace(out, r"$" => "")

    open(output, "w") do io
        documenter && println(io, "````````````@raw html")
        println(io, out)
        documenter && println(io, "````````````")
    end
end

function render_readme_html(pkgroot, buildpath, src_prefix="", href_prefix=""; documenter = false)
    outpath = joinpath(buildpath, "_readme")
    try
        readme = find_readme(pkgroot)

        mkpath(outpath)
        readmehtml = joinpath(outpath, "readme.html")

        @info("Rendering readme to HTML.")
        render_html(readme, readmehtml, src_prefix, href_prefix; documenter = false)
        @info("Done rendering readme to HTML.")

        return readmehtml
    catch err
        @error("Error trying to render readme to HTML.", exception=err)
    end
end

function should_rewrite_url(url)
    if occursin(r"^\.?\.?//?"i, url)
        return true
    elseif startswith(url, '#')
        return false
    else
        return !occursin(r"^\w+://"i, url)
    end
end

function replace_url(el, attr, url_prefix; fix_github = false)
    old_src = getattr(el, attr)
    if should_rewrite_url(old_src)
        setattr!(el, attr, string(url_prefix, old_src))
    end
    if fix_github
        url = getattr(el, attr)
        setattr!(el, attr, github_to_raw(url))
    end
end

function github_to_raw(url)
    m = match(r"(https?)\:\/\/(?:www\.)?github.com/(.*?)/(.*?)/blob/(.*)", url)
    if m !== nothing
        return string(m[1], "://raw.githubusercontent.com/", m[2], "/", m[3], "/", m[4])
    end
    return url
end

function postprocess_html_readme(html; src_prefix="", href_prefix="")
    doc = parsehtml(html, preserve_whitespace = true).root

    if !endswith(src_prefix, '/')
        src_prefix = string(src_prefix, '/')
    end
    if !endswith(href_prefix, '/')
        href_prefix = string(href_prefix, '/')
    end

    header_refs = Set([])
    for el in PreOrderDFS(doc)
        if el isa HTMLElement
            if Gumbo.tag(el) in [:h1, :h2, :h3, :h4, :h5]
                heading = replace(text(el), r"\s" => "-")
                heading = replace(heading, r"[^\w-]" => "")

                while heading in header_refs
                    new_heading = replace(heading, r"_\d+$" => s -> begin
                        string('_', parse(Int, s[2:end]) + 1)
                    end)
                    if new_heading == heading
                        heading = string(heading, "_1")
                    else
                        heading = new_heading
                    end
                end

                orig_content = deepcopy(el.children)
                a_el = HTMLElement{:a}(orig_content, el, Dict(
                    "href" => string('#', heading),
                    "class" => "docs-heading-anchor"
                ))

                for c in orig_content
                    c.parent = a_el
                end

                el.children = [a_el]

                setattr!(el, "id", heading)
                push!(header_refs, heading)
            elseif Gumbo.tag(el) == :code
                if Gumbo.tag(el.parent) == :pre && length(el.children) == 1 && typeof(el.children[1]) == HTMLText
                    highlight_syntax_html(el)
                end
            elseif hasattr(el, "src")
                replace_url(el, "src", src_prefix, fix_github = true)
            elseif hasattr(el, "href")
                replace_url(el, "href", href_prefix)
            end
        end
    end

    return doc
end

function highlight_syntax_html(el)
    lexer = nothing
    content = Gumbo.text(el)
    m = match(r"language\-(.+)\b", hasattr(el, "class") ? getattr(el, "class") : "")
    if m ≠ nothing
        # language specified
        m = m[1]
        if m == "julia" || m == "jl"
            lexer = Lexers.JuliaLexer
        elseif m == "julia-console" || m == "jldoctest" || m == "jl-doctest" || m == "jl-console" || m == "jl-repl"
            lexer = Lexers.JuliaConsoleLexer
        elseif m == "matlab"
            lexer = Lexers.MatlabLexer
        elseif m == "r"
            lexer = Lexers.RLexer
        elseif m == "fortran"
            lexer = Lexers.FortranLexer
        end
    else
        # no language, but looks like REPL content
        if startswith(content, "julia>")
            lexer = Lexers.JuliaConsoleLexer
        end
    end

    if lexer ≠ nothing
        io = IOBuffer()
        Highlights.highlight(io, MIME("text/html"), content, lexer)
        highlighted_content = String(take!(io))
        parsed = Gumbo.parsehtml(highlighted_content, preserve_whitespace = true).root.children[2].children[1]
        parsed.parent = el.parent
        el.children = [parsed]
    end
end


text:
"""
    read_block(filename::String, blockname::String;
                                info::InfoLine=InfoLine(),
                                parttype::Integer=-1)

Reads a block in a snapshot with given name. Names are case sensitive.

# Examples
```julia
julia> pos_info = InfoLine("POS", Float32, 1, [1, 1, 1, 1, 1, 1])
[...]
julia> gas_pos = read_block(filename, "POS", info=pos_info, parttype=0)
[...]
```
"""
function read_block(filename::String, blockname::String;
                    parttype::Integer=-1,
                    block_position::Integer=-1,
                    info::Union{Nothing,InfoLine}=nothing,
                    h::Union{Nothing,SnapshotHeader}=nothing)

    # the file is actually a filebase -> return concatenated arrays
    if !isfile(filename)
        return read_block_subsnaps(filename, blockname; parttype, info, h)
    end

    # we need to know which particle to read
    if parttype == -1
        error("Please specify particle type!")
    end

    if isnothing(h)
        # read header - super fast and needed for flexibility
        h = head_to_obj(filename)
    end

    # check if particle type is present
    if h.nall[parttype+1] == UInt32(0)
        error("Particle Type $parttype not present in simulation!")
    end

    # check if info is present
    if isnothing(info)
        info = check_info(filename, blockname)
    end

    # check if a block position is supplied
    if block_position == -1

        block_position = check_block_position(filename, blockname)

        if block_position == -1 || (info.is_present[parttype+1] == 0)
            # if no mass block is present we can read it from the header
            if blockname == "MASS"
                block = Array{info.data_type,1}(undef, h.npart[parttype+1])
                block .= h.massarr[parttype+1]
                return block
            else
                error("Block $blockname not present for particle type $parttype !")
            end
        end
    end

    f = open(filename)

    seek(f,block_position)

    for i ∈ 1:size(info.is_present,1)
        p = position(f)

        if info.is_present[i] == Int32(1)
            if i == (parttype+1)
                block = read_block_data(f, info.data_type, info.n_dim, h.npart[i])
                close(f)
                return block
            else
                seek(f, p + ( sizeof(info.data_type)*info.n_dim*h.npart[i] ))
            end # if i == (parttype+1)
        end # info.is_present[i] == Int32(1)
    end # i ∈ 1:size(info.is_present,1)

    close(f)

end

"""
    read_block_subsnaps(filename::String, blockname::String;
                                info::InfoLine=InfoLine(),
                                parttype::Integer=-1)

Reads a block over all sub-snapshots. Names are case sensitive.

# Examples
```julia
julia> pos_info = InfoLine("POS", Float32, 1, [1, 1, 1, 1, 1, 1])
[...]
julia> gas_pos = read_block(filename, "POS", info=pos_info, parttype=0)
[...]
```
"""
function read_block_subsnaps(filebase::String, blockname::String;
                             parttype::Integer=-1,
                             info::Union{Nothing,InfoLine}=nothing,
                             h::Union{Nothing,SnapshotHeader}=nothing)

    # the file is actually a filebase -> return concatenated arrays
    if !isfile(filebase * ".0")
        error("Neither $filebase nor $(filebase * ".0") present!")
    end

    # we need to know which particle to read
    if parttype == -1
        error("Please specify particle type!")
    end

    if isnothing(h)
        # read header - super fast and needed for flexibility
        h = head_to_obj(filebase * ".0")
    end

    # check if info is present
    if isnothing(info)
        info = check_info(filebase * ".0", blockname)
    end

    return read_subsnaps(filebase, blockname, parttype, info, h)

end

"""
    read_subsnaps(filebase::String, blockname::String, parttype::Integer,
                          info::InfoLine, h_global::SnapshotHeader)

Reads a block over distributed files and returns it in one large Array.
"""
function read_subsnaps(filebase::String, blockname::String, parttype::Integer,
                          info::InfoLine, h_global::SnapshotHeader)

    # allocate empty array for block
    if info.n_dim > 1
        block = Array{info.data_type,2}(undef, ( info.n_dim, h_global.nall[parttype+1] ))
    else
        block = Array{info.data_type,1}(undef, h_global.nall[parttype+1])
    end

    # store number of particles that have been read
    N_read = 0

    # loop over all sub-files
    @inbounds for i = 0:(h_global.num_files-1)

        # get filename of the current file
        filename = "$filebase.$i"

        # read header of subfile
        h = read_header(filename)

        # get number of particles for this sub-file
        N_to_read = h.npart[parttype+1]

        # read the block
        if info.n_dim > 1
            block[:, N_read+1:N_read+N_to_read] = read_block(filename, blockname;
                                                            parttype, info, h)
        else
            block[N_read+1:N_read+N_to_read]    = read_block(filename, blockname;
                                                            parttype, info, h)
        end

        # update number of read particles
        N_read += N_to_read
    end

    return block
end

"""
    read_block_data(f::IOStream, data_type::DataType, n_dim::Integer, npart::Integer)

Reads the binary data in a block.
"""
function read_block_data(f::IOStream, data_type::DataType, n_dim::Integer, npart::Integer)

    if n_dim > 1
        return read!(f, Array{data_type,2}(undef, (n_dim, npart)))
    else
        read!(f, Array{data_type,1}(undef, npart))
    end
end

"""
    check_info(filename::String, blockname::String)

Helper function to read INFO block or construct `InfoLine` for MASS block, if no INFO block is present.
Returns a single `InfoLine` struct.
"""
function check_info(filename::String, blockname::String)
    info = read_info(filename)
    if info == 1
        if blockname == "MASS"
            return InfoLine("MASS", Float32, 1, [0, 0, 0, 0, 0, 0])
        else
            error("No Info block in snapshot! Supply InfoLine type!")
        end
    else # info != 1
        for i ∈ 1:size(info,1)
            if info[i].block_name == blockname
                return info[i]
            end # if block found
        end # loop over info
        if isa(info, Array)
            if (blockname == "MASS")
                return InfoLine("MASS", Float32, 1, [0, 0, 0, 0, 0, 0])
            else
                error("Block $blockname not present!")
            end
        end
    end # info == 1
end

"""
    check_block_position(filename::String, blockname::String)

Helper function to find read positions.
"""
function check_block_position(filename::String, blockname::String)
    # read positions of all blocks
    block_positions_dict = get_block_positions(filename)

    # check if the requested block is present
    if haskey(block_positions_dict, blockname)
        return block_positions_dict[blockname]
    else
        if blockname != "MASS"
            # if the block is not present we need error handling!
            error("Requested block $blockname not present!")
        else
            return -1
        end

    end
end

"""
    read_block_with_offset!(data, n_read::Integer, filename::String, pos0::Integer, info::InfoLine,
                                offset::Integer, offset_key::Array{<:Integer},
                                part_per_key::Array{<:Integer} )

Read part of a block to pre-allocated array.
"""
function read_block_with_offset!(data, n_read::Integer, filename::String, pos0::Integer, info::InfoLine,
                                offset::Integer, offset_key::Array{<:Integer},
                                part_per_key::Array{<:Integer} )

    # open the file
    f = open(filename)

    # number of bits in data_type
    len = sizeof(info.data_type) * info.n_dim

    # jump to position of particle type in relevant block
    seek(f, pos0+offset*len)

    # store position in file
    p = position(f)

    n_this_key = n_read

    for i = 1:size(offset_key,1)

        # jump to start of key
        seek(f, p + len*offset_key[i])
        n_this_key += part_per_key[i]

        if info.n_dim == 1
            data[n_read+1:n_this_key]    = read_block_data(f, info.data_type, info.n_dim, part_per_key[i])
        else
            data[:, n_read+1:n_this_key] = read_block_data(f, info.data_type, info.n_dim, part_per_key[i])
        end

        n_read += part_per_key[i]

    end # for

    # close the file
    close(f)

    data
end


text:
push!(LOAD_PATH,"../src/")

using Documenter, AntiIdClosedGroup

makedocs(
    sitename="AntiIdClosedGroup.jl",
    authors = "Felix Höddinghaus",
    pages = [
        "index.md",
        "anti-identity.md",
        "notation.md",
        "generators.md",
        "permutations.md",
        "conjugacy-invariants.md",
        "example.md",
    ],
)

# show documentation locally after generating with
#   python3 -m http.server --bind localhost


text:
using Test
using JET

using IntervalSets

import ClimaCore: Domains, Meshes, Spaces, Fields, Operators
import ClimaCore.Domains: Geometry

# We need to pull these broadcasted expressions out as
# toplevel functions due to how broadcast expressions are
# lowered so JETTest can have a single callsite to analyze.

function opt_InterpolateF2C(face_field)
    I = Operators.InterpolateF2C()
    return I.(identity.(face_field))
end

function opt_WeightedInterpolateF2C(weights, face_field)
    WI = Operators.WeightedInterpolateF2C()
    return identity.(WI.(weights, face_field))
end

function opt_LeftBiasedF2C(face_field)
    LB = Operators.LeftBiasedF2C(left = Operators.SetValue(0.0))
    return LB.(identity.(face_field))
end

function opt_RightBiasedF2C(face_field)
    RB = Operators.RightBiasedF2C(right = Operators.SetValue(0.0))
    return RB.(identity.(face_field))
end

function opt_AdvectionF2F(face_vel, face_field)
    A = Operators.AdvectionF2F()
    return A.(face_vel, identity.(face_field))
end

function opt_FluxCorrectionF2F_Extrapolate(face_vel, face_field)
    FC = Operators.FluxCorrectionF2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    return FC.(face_vel, identity.(face_field))
end

function opt_GradientF2C(face_field)
    ∇ᶜ = Operators.GradientF2C()
    return Geometry.WVector.(∇ᶜ.(sin.(face_field)))
end

function opt_DivergenceF2C(face_field)
    divᶜ = Operators.DivergenceF2C()
    return divᶜ.(Geometry.WVector.(sin.(face_field)))
end

function opt_SetBoundary_SetValue(face_field)
    B = Operators.SetBoundaryOperator(
        left = Operators.SetValue(0.0),
        right = Operators.SetValue(0.0),
    )
    return B.(sin.(face_field))
end

function opt_InterpolateC2F_SetValue(center_field)
    I = Operators.InterpolateC2F(
        left = Operators.SetValue(0.0),
        right = Operators.SetValue(0.0),
    )
    return I.(identity.(center_field))
end

function opt_InterpolateC2F_SetGradient(center_field)
    I = Operators.InterpolateC2F(
        left = Operators.SetGradient(Geometry.WVector(0.0)),
        right = Operators.SetGradient(Geometry.WVector(0.0)),
    )
    return I.(identity.(center_field))
end

function opt_InterpolateC2F_Extrapolate(center_field)
    I = Operators.InterpolateC2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    return I.(identity.(center_field))
end

function opt_WeightedInterpolateC2F_SetValue(weights, center_field)
    WI = Operators.WeightedInterpolateC2F(
        left = Operators.SetValue(0.0),
        right = Operators.SetValue(0.0),
    )
    return identity.(WI.(weights, center_field))
end

function opt_WeightedInterpolateC2F_SetGradient(weights, center_field)
    WI = Operators.WeightedInterpolateC2F(
        left = Operators.SetGradient(Geometry.WVector(0.0)),
        right = Operators.SetGradient(Geometry.WVector(0.0)),
    )
    return identity.(WI.(weights, center_field))
end

function opt_WeightedInterpolateC2F_Extrapolate(weights, center_field)
    WI = Operators.WeightedInterpolateC2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    return identity.(WI.(weights, center_field))
end

function opt_LeftBiasedC2F(center_field)
    LB = Operators.LeftBiasedC2F(left = Operators.SetValue(0.0))
    return LB.(identity.(center_field))
end

function opt_RightBiasedC2F(center_field)
    RB = Operators.RightBiasedC2F(right = Operators.SetValue(0.0))
    return RB.(identity.(center_field))
end

function opt_UpwindBiasedProductC2F_SetValue(face_vel, center_field)
    UB = Operators.UpwindBiasedProductC2F(
        left = Operators.SetValue(0.0),
        right = Operators.SetValue(0.0),
    )
    return UB.(face_vel, identity.(center_field))
end

function opt_UpwindBiasedProductC2F_Extrapolate(face_vel, center_field)
    UB = Operators.UpwindBiasedProductC2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    return UB.(face_vel, identity.(center_field))
end

function opt_AdvectionC2C_SetValue(face_vel, center_field)
    A = Operators.AdvectionC2C(
        left = Operators.SetValue(0.0),
        right = Operators.SetValue(0.0),
    )
    return A.(face_vel, identity.(center_field))
end

function opt_AdvectionC2C_Extrapolate(face_vel, center_field)
    A = Operators.AdvectionC2C(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    return A.(face_vel, identity.(center_field))
end

function opt_FluxCorrectionC2C_Extrapolate(face_vel, center_field)
    FC = Operators.FluxCorrectionC2C(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    return FC.(face_vel, identity.(center_field))
end

function opt_GradientC2F_SetValue(center_field)
    ∇ᶠ = Operators.GradientC2F(
        left = Operators.SetValue(1.0),
        right = Operators.SetValue(-1.0),
    )
    return Geometry.WVector.(∇ᶠ.(cos.(center_field)))
end

function opt_GradientC2F_SetGradient(center_field)
    ∇ᶠ = Operators.GradientC2F(
        left = Operators.SetGradient(Geometry.WVector(0.0)),
        right = Operators.SetGradient(Geometry.WVector(0.0)),
    )
    return Geometry.WVector.(∇ᶠ.(cos.(center_field)))
end

function opt_DivergenceC2F_SetValue(center_field)
    divᶠ = Operators.DivergenceC2F(
        left = Operators.SetValue(Geometry.WVector(0.0)),
        right = Operators.SetValue(Geometry.WVector(0.0)),
    )
    return divᶠ.(Geometry.WVector.(sin.(center_field)))
end

function opt_DivergenceC2F_SetDivergence(center_field)
    # DivergenceC2F, SetDivergence
    divᶠ = Operators.DivergenceC2F(
        left = Operators.SetDivergence(0.0),
        right = Operators.SetDivergence(0.0),
    )
    return divᶠ.(Geometry.WVector.(cos.(center_field)))
end

function opt_CurlC2F_SetValue(center_field)
    # DivergenceC2F, SetDivergence
    curlᶠ = Operators.CurlC2F(
        left = Operators.SetValue(Geometry.Covariant1Vector(0.0)),
        right = Operators.SetValue(Geometry.Covariant1Vector(0.0)),
    )
    return curlᶠ.(Geometry.Covariant1Vector.(cos.(center_field)))
end

# Test that Julia ia able to optimize Stencil operations v1.7+
@static if @isdefined(var"@test_opt")
    @testset "Scalar Field FiniteDifferenceSpaces optimizations" begin
        for FT in (Float64,)
            domain = Domains.IntervalDomain(
                Geometry.ZPoint{FT}(0.0),
                Geometry.ZPoint{FT}(pi);
                boundary_tags = (:left, :right),
            )
            mesh = Meshes.IntervalMesh(domain; nelems = 16)

            center_space = Spaces.CenterFiniteDifferenceSpace(mesh)
            face_space = Spaces.FaceFiniteDifferenceSpace(center_space)

            faces = getproperty(Fields.coordinate_field(face_space), :z)
            face_values = ones(FT, face_space)
            face_velocities = Geometry.WVector.(face_values)

            centers = getproperty(Fields.coordinate_field(center_space), :z)
            center_values = ones(FT, center_space)
            center_velocities = Geometry.WVector.(center_values)

            # face space operators
            @test_opt sum(ones(FT, face_space))
            @test_opt sum(sin.(faces))

            @test_opt opt_InterpolateF2C(faces)
            @test_opt opt_WeightedInterpolateF2C(face_values, faces)

            @test_opt opt_LeftBiasedF2C(faces)
            @test_opt opt_RightBiasedF2C(faces)

            @test_opt opt_AdvectionF2F(face_velocities, faces)

            @test_opt opt_FluxCorrectionF2F_Extrapolate(
                center_velocities,
                faces,
            )

            @test_opt opt_GradientF2C(faces)
            @test_opt opt_DivergenceF2C(faces)

            @test_opt opt_SetBoundary_SetValue(faces)

            # center space operators
            @test_opt sum(ones(FT, center_space))
            @test_opt sum(sin.(centers))

            @test_opt opt_InterpolateC2F_SetValue(centers)
            @test_opt opt_InterpolateC2F_SetGradient(centers)
            @test_opt opt_InterpolateC2F_Extrapolate(centers)

            @test_opt opt_WeightedInterpolateC2F_SetValue(
                center_values,
                centers,
            )
            @test_opt opt_WeightedInterpolateC2F_SetGradient(
                center_values,
                centers,
            )
            @test_opt opt_WeightedInterpolateC2F_Extrapolate(
                center_values,
                centers,
            )

            @test_opt opt_LeftBiasedC2F(centers)
            @test_opt opt_RightBiasedC2F(centers)

            @test_opt opt_UpwindBiasedProductC2F_SetValue(
                face_velocities,
                centers,
            )
            @test_opt opt_UpwindBiasedProductC2F_Extrapolate(
                face_velocities,
                centers,
            )

            @test_opt opt_AdvectionC2C_SetValue(face_velocities, centers)
            @test_opt opt_AdvectionC2C_Extrapolate(face_velocities, centers)

            @test_opt opt_FluxCorrectionC2C_Extrapolate(
                face_velocities,
                centers,
            )

            @test_opt opt_GradientC2F_SetValue(centers)
            @test_opt opt_GradientC2F_SetGradient(centers)

            @test_opt opt_DivergenceC2F_SetValue(centers)
            @test_opt opt_DivergenceC2F_SetDivergence(centers)
            @test_opt opt_CurlC2F_SetValue(centers)
        end
    end
end


text:
function movements_mark!(r,side::HorizonSide)
    while (isborder(r,side) ==  false) && (ismarker(r) == false)
        putmarker!(r)
        move!(r,side)
    end
end

function movements_to_side!(r,side::HorizonSide)
    count = 0
    while isborder(r,side) ==  false
        move!(r,side)
        count = count + 1
    end
    return count
end

function movements_to_zero!(r, side::HorizonSide, count)

    for i = 1:count
        move!(r,side)
    end
end

function kvadro!(r)

    count = movements_to_side!(r,Ost)

    for side in [Sud, West, Nord, Ost, Sud]
        movements_mark!(r,side)
    end

    movements_to_zero!(r, West, count)

end

kvadro!(r)
show(r)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
# module for fast slam 1.0
# estimate robot's pose and landmark's position on map

using Distributions, LinearAlgebra, StatsBase

include(joinpath(split(@__FILE__, "src")[1], "src/model/map/map.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/slam/map_particle.jl"))

mutable struct FastSlam1
  particles
  motion_noise_rate_pdf
  dist_dev
  dir_dev
  max_likelihood_particle
  estimated_pose

  function FastSlam1(init_pose; particle_num=100, objects_num=0,
                     motion_noise_stds=Dict("nn"=>0.20, "no"=>0.001, "on"=>0.11, "oo"=>0.20),
                     dist_dev_rate=0.14, dir_dev=0.05)
    self = new()
    self.particles = [MapParticle(init_pose, 1.0/particle_num, objects_num) for i in 1:particle_num]
    self.dist_dev = dist_dev_rate
    self.dir_dev = dir_dev
    self.max_likelihood_particle = self.particles[1]
    self.estimated_pose = self.max_likelihood_particle.pose

    v = motion_noise_stds
    c = diagm(0 => [v["nn"]^2, v["no"]^2, v["on"]^2, v["oo"]^2])
    self.motion_noise_rate_pdf = MvNormal(c)
    return self
  end
end

function motion_update(self::FastSlam1, speed, yaw_rate, time_interval)
  for p in self.particles
    motion_update(p, speed, yaw_rate, time_interval, self.motion_noise_rate_pdf)
  end
end

function set_max_likelihood_pose(self::FastSlam1)
  max_index = argmax([p.weight for p in self.particles])
  self.max_likelihood_particle = self.particles[max_index]
  self.estimated_pose = self.max_likelihood_particle.pose
end

function observation_update(self::FastSlam1, observation)
  for p in self.particles
    observation_update(p, observation, self.dist_dev, self.dir_dev, true)
  end
  set_max_likelihood_pose(self)
  resampling(self)
end

function resampling(self::FastSlam1)
  ws = [(if p.weight <= 0.0 1e-100 else p.weight end) for p in self.particles]

  ps = sample(self.particles, Weights(ws), length(self.particles))
  self.particles = [deepcopy(e) for e in ps]
  for p in self.particles
    p.weight = 1.0/length(self.particles)
  end
end

function draw!(self::FastSlam1)
  k = 0.5 # scale for length of arrows
  # all of particles
  px = [p.pose[1] for p in self.particles]
  py = [p.pose[2] for p in self.particles]
  vx = [cos(p.pose[3])*k for p in self.particles]
  vy = [sin(p.pose[3])*k for p in self.particles]
  quiver!(px, py, quiver=(vx, vy), aspect_ratio=true, color="blue")

  # maximum likelihood particle
  mx = [self.estimated_pose[1]]
  my = [self.estimated_pose[2]]
  mvx = [cos(self.estimated_pose[3])*k]
  mvy = [sin(self.estimated_pose[3])*k]
  quiver!(mx, my, quiver=(mvx, mvy), aspect_ratio=true, color="red")

  # estimated landmarks on map
  # the max likelihood particle
  for elm in self.max_likelihood_particle.map.objects
    draw!(elm)
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using .QuantumOptics

function dirac(io::IO, state::Union{Ket, Bra}, statename="ψ"; header::Bool=false)
    header && summary(io, state)
    data = state.data
    shape = reverse(state.basis.shape)
    if _islatex && isdefined(Main, :IJulia) && Main.IJulia.inited # for IJulia rendering
        if statename == "ψ"
            statename = "\\psi"
        end
        if state isa Ket
            str = "\$" * sprint(io -> print_dirac(io, data, shape, statename)) * "\$"
        else
            str = "\$" * sprint(io -> print_dirac(io, transpose(data), shape, statename)) * "\$"
        end
        display("text/markdown", str)
    else
        if state isa Ket
            print_dirac(io, data, shape, statename)
        else
            print_dirac(io, transpose(data), shape, statename)
        end
    end
end
dirac(state::Union{Ket, Bra}, statename="ψ"; header::Bool=false) = dirac(stdout, state, statename, header=header)

function dirac(io::IO, state::Union{DenseOperator, SparseOperator}, statename="ρ"; header::Bool=false)
    if header
        summary(io, state)
        println(io)
    end
    data = state.data
    lshape = reverse(state.basis_l.shape)
    rshape = reverse(state.basis_r.shape)
    if _islatex && isdefined(Main, :IJulia) && Main.IJulia.inited # for IJulia rendering
        if statename == "ρ"
            statename = "\\rho"
        end
        str = "\$" * sprint(io -> print_dirac(io, data, lshape, rshape, statename)) * "\$"

        display("text/markdown", str)
    else
        print_dirac(io, data, lshape, rshape, statename)
    end
end
dirac(state::Union{DenseOperator, SparseOperator}, statename="ρ"; header::Bool=false) = dirac(stdout, state, statename, header=header)


text:
/* test input */

int main() {
  int x = readInt();
  double y = readDouble();

  printInt(x-5);

  if (y > 42.0 || y < 43.0)
    printString("yay!");
  else
    printString("nay!");
  return 0 ;

}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
export ⋅, normalize, slerp

Base.:(+)(q₁::Quaternion{T}, q₂::Quaternion{U}) where {T,U} = Quaternion{promote_type(T,U)}(q₁.v + q₂.v, q₁.w + q₂.w)

Base.:(-)(q₁::Quaternion{T}, q₂::Quaternion{U}) where {T,U} = Quaternion{promote_type(T,U)}(q₁.v - q₂.v, q₁.w - q₂.w)

function Base.:(*)(q::Quaternion{T}, c::U) where {T, U<:Number}
    @assert !isnan(c)
    Quaternion{promote_type(T,U)}(q.v * c, q.w * c)
end

Base.:(*)(c::U, q::Quaternion{T}) where {T,U} = q * c

function Base.:(/)(q::Quaternion{U}, c::T) where {T<:Number,U}
    @assert c != zero(T)
    q * (one(U) / c)
end

⋅(q₁::Quaternion, q₂::Quaternion) = q₁.v ⋅ q₂.v + q₁.w * q₂.w

normalize(q::Quaternion) = q / √(q ⋅ q)

function slerp(t::T, q₁::Quaternion, q₂::Quaternion) where {T<:Number,U,V}
    cosθ = q₁ ⋅ q₂
    if cosθ > 0.9995
        return normalize((one(T) - t) * q₁ + t * q₂)
    else
        θ = acos(clamp(cosθ, -one(T), one(T)))
        θₚ = θ * t
        q⟂ = normalize(q₂ - q₁ * cosθ)
        return q₁ * cos(θₚ) + q⟂ * sin(θₚ)
    end
end


text:
using Brillouin
using Test

# defines `test_show` to test `show` methods with good error diffs, using DeepDiffs.jl
include("testutils_show.jl")

# ---------------------------------------------------------------------------------------- #
# core API/utilities
include("kpaths.jl")
include("wignerseitz.jl")
include("kpaths_spglib.jl")

# ---------------------------------------------------------------------------------------- #
# optional plot utilities
include("plotlyjs.jl")


text:
using Soss

multiple_poisson_regression = @model (traps, live_in_super, log_sq_foot) begin
    beta ∼ Normal(-0.25, 1)
    beta_super ∼ Normal(-0.5, 1)
    alpha ∼ Normal(log(4), 1)

    complaints ∼ For(eachindex(traps)) do i
        # TODO: Numerically safe version
        λ = exp(alpha + beta * traps[i] + beta_super * live_in_super[i] + log_sq_foot[i])
        Poisson(λ)
    end
end


text:
############################################################
# BinaryCrossEntropyLossLayer
############################################################
@defstruct BinaryCrossEntropyLossLayer Layer (
  name :: AbstractString = "binary-cross-entropy-loss",
  (weight :: AbstractFloat = 1.0, weight >= 0),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == 2),
)
@characterize_layer(BinaryCrossEntropyLossLayer,
  has_loss => true,
  is_sink  => true,
  can_do_bp => true,
)

type BinaryCrossEntropyLossLayerState{T} <: LayerState
  layer :: BinaryCrossEntropyLossLayer
  loss  :: T
end

function setup(backend::Backend, layer::BinaryCrossEntropyLossLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  data_type = eltype(inputs[1])
  state = BinaryCrossEntropyLossLayerState(layer, convert(data_type, 0))
  return state
end

function shutdown(backend::Backend, state::BinaryCrossEntropyLossLayerState)
end

function forward(backend::CPUBackend, state::BinaryCrossEntropyLossLayerState, inputs::Vector{Blob})
  pred = vec(inputs[1].data)
  label = vec(inputs[2].data)
  loss = BLAS.dot(log(pred), label) + BLAS.dot(log1p(-pred), (1-label))

  num = get_num(inputs[1])
  state.loss = state.layer.weight * -loss/num
end

function backward(backend::CPUBackend, state::BinaryCrossEntropyLossLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  diff = diffs[1]
  if isa(diff, CPUBlob)

    # Diffs is df/dloss
    # we want to multiply this by dloss/dpred and dl

    label = inputs[2].data
    pred  = reshape(inputs[1].data, size(label))

    # l = -w sum_i log(p_i) y_i + (log(1-p_i)(1-y_i)
    # dl/dp_i = -w ( y_i/p_i - ((1-y_i)(1-p_i) )
    # dl/dy_i = -w ( log(p_i) - (log(1-p_i) ) = -w log(p_i/(1-p_i))

    n = length(pred)
    a = convert(eltype(pred), -state.layer.weight/get_num(inputs[1]))

    erase!(diff) # is this correct? square-loss does it - should we allow for any incoming diff?
    dl_dpred = (label ./ pred) - ((1-label) ./ (1-pred)) # dloss/d?
    BLAS.axpy!(n, a, dl_dpred, 1, diff.data, 1)
  end
  diff = diffs[2]
  if isa(diff, CPUBlob)
    dl_dlabel = log(pred ./ (1-pred))
    erase!(diff) # is this correct? square-loss does it
    BLAS.axpy!(n, a, dl_dlabel, 1, diff.data, 1)
  end
end


text:
module PrintingTests

using Test
using WeaklyHard

@testset "Print automaton" begin
    g1 = build_automaton(AnyHitConstraint(1, 3))
    @test_nowarn println(Base.DevNull(), 1)
    @test_nowarn println(Base.DevNull(), g1.head)

    g2 = build_automaton(AnyHitConstraint(10, 16))
    @test_nowarn println(Base.DevNull(), g2)
    @test_nowarn println(Base.DevNull(), g2.head)

    g3 = build_automaton(RowHitConstraint(10, 60))
    @test_nowarn println(Base.DevNull(), g3)
    @test_nowarn println(Base.DevNull(), g3.head)
end #testset

@testset "Print BigInt" begin
    g1 = build_automaton(AnyHitConstraint(1, 3))
    seq1 = random_sequence(g1, 200)
    @test_nowarn bitstring(seq1)
    seq2 = random_sequence(g1, 100_000)
    @test_nowarn bitstring(seq2, 10)
end #testset

@testset "Print Constraint" begin
    l1 = AnyHitConstraint(1, 3)
    @test_nowarn println(Base.DevNull(), l1)
    l2 = AnyMissConstraint(1, 3)
    @test_nowarn println(Base.DevNull(), l2)
    l3 = RowHitConstraint(1, 3)
    @test_nowarn println(Base.DevNull(), l3)
    l4 = RowMissConstraint(1, 3)
    @test_nowarn println(Base.DevNull(), l4)
end #testset

end # module


text:
module AO

using DataDeps              # For storage location of data

function fallback_download(remotepath, localdir)
    @assert(isdir(localdir))
    filename = basename(remotepath)  # only works for URLs with filename as last part of name
    localpath = joinpath(localdir, filename)
    Base.download(remotepath, localpath)
    return localpath
end

function register_AO_files()
    register(
        DataDep(
            "AWESOME-OCIM",
            """
            References for OCIM1:
            - John, S. G., Liang, H., Weber, T., Devries, T., Primeau, F., Moore, K., Holzer, M., Mahowald, N., Gardner, W., Mishonov, A., Richardson, M., J., Faugere, Y., and Taburet, G. (2020). AWESOME OCIM: A simple, flexible, and powerful tool for modeling elemental cycling in the oceans. Chemical Geology, 533, 119403. doi: 10.1016/j.chemgeo.2019.119403.
            """,
            "https://github.com/hengdiliang/AWESOME-OCIM/archive/refs/heads/master.zip",
            sha2_256,
            fetch_method = fallback_download,
            post_fetch_method = unpack
        )
    )
    return nothing
end

"""
    download_and_unpack

Downloads and unpacks the AO zip file from the MTEL website.
Returns the local path of where the AO is installed for you.
"""
function download_and_unpack()
    register_AO_files()
    AO_path = datadep"AWESOME-OCIM"
    @info """You are about to download (and unpack) the AWESOME OCIM (AO)
          files as a zip file from GitHub and unpack them into:

              $AO_path

          You can run

              tree $AO_path

          (in shell mode) to check its contents.

          Please check with Seth John (sethjohn@usc.edu) for references
          to cite if you use the other data contained in the AO files,
          e.g., GEOTRACES, WOA, Weber and John, and so on.

          Also note that downloading files like this from GitHub is
          probably not very robust. However, some of the data files in
          the AO repository do not exist anywhere else officially,
          therefore this seems like the best solution at this stage.
          PRs welcome to improve this!
          """
    return AO_path
end

end # module

export AO