text:
module TestRecode
    using Base.Test
    using CategoricalArrays
    using CategoricalArrays: DefaultRefType, CatArray, CatVector
    using NullableArrays

    ## Test recode!, used by recode

    # Test both recoding into x itself and into an uninitialized vector
    # (since for CategoricalVectors possible bugs can happen when working in-place)

    # Recoding from Int to Int
    for x in ([1:10;], CategoricalArray(1:10)),
        y in (similar(x), Array{Int}(size(x)), CategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 1=>100, 2:4=>0, [5; 9:10]=>-1)
        @test y === z
        @test y == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
        if isa(y, CatArray)
            @test levels(y) == [6, 7, 8, 100, 0, -1]
            @test !isordered(y)
        end
    end
    for x in (Nullable{Int}[1:10;], NullableArray(1:10), NullableCategoricalArray(1:10)),
        y in (similar(x), Array{Nullable{Int}}(size(x)),
              NullableArray{Int}(size(x)), NullableCategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 1=>100, 2:4=>0, [5; 9:10]=>-1)
        @test y === z
        @test isequal(y, Nullable{Int}[100, 0, 0, 0, -1, 6, 7, 8, -1, -1])
        if isa(y, CatArray)
            @test levels(y) == [6, 7, 8, 100, 0, -1]
            @test !isordered(y)
        end
    end

    # Recoding from Int to Int with duplicate recoded values
    for x in ([1:10;], CategoricalArray(1:10)),
        y in (similar(x), Array{Int}(size(x)), CategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 1=>100, 2:4=>100, [5; 9:10]=>-1)
        @test y === z
        @test y == [100, 100, 100, 100, -1, 6, 7, 8, -1, -1]
        if isa(y, CatArray)
            @test levels(y) == [6, 7, 8, 100, -1]
            @test !isordered(y)
        end
    end
    for x in (Nullable{Int}[1:10;], NullableArray(1:10), NullableCategoricalArray(1:10)),
        y in (similar(x), Array{Nullable{Int}}(size(x)),
              NullableArray{Int}(size(x)), NullableCategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 1=>100, 2:4=>100, [5; 9:10]=>-1)
        @test y === z
        @test isequal(y, Nullable{Int}[100, 100, 100, 100, -1, 6, 7, 8, -1, -1])
        if isa(y, CatArray)
            @test levels(y) == [6, 7, 8, 100, -1]
            @test !isordered(y)
        end
    end

    # Recoding from Int to Int with unused level
    for x in ([1:10;], CategoricalArray(1:10)),
        y in (similar(x), Array{Int}(size(x)), CategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 1=>100, 2:4=>0, [5; 9:10]=>-1, 100=>1)
        @test y === z
        @test y == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
        if isa(y, CatArray)
            @test levels(y) == [6, 7, 8, 100, 0, -1, 1]
            @test !isordered(y)
        end
    end
    for x in (Nullable{Int}[1:10;], NullableArray(1:10), NullableCategoricalArray(1:10)),
        y in (similar(x), Array{Nullable{Int}}(size(x)),
              NullableArray{Int}(size(x)), NullableCategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 1=>100, 2:4=>0, [5; 9:10]=>-1, 100=>1)
        @test y === z
        @test isequal(y, Nullable{Int}[100, 0, 0, 0, -1, 6, 7, 8, -1, -1])
        if isa(y, CatArray)
            @test levels(y) == [6, 7, 8, 100, 0, -1, 1]
            @test !isordered(y)
        end
    end

    # Recoding from Int to Int with duplicate default
    for x in ([1:10;], CategoricalArray(1:10)),
        y in (similar(x), Array{Int}(size(x)), CategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 100, 1=>100, 2:4=>100, [5; 9:10]=>-1)
        @test y === z
        @test y == [100, 100, 100, 100, -1, 100, 100, 100, -1, -1]
        if isa(y, CatArray)
            @test levels(y) == [100, -1]
            @test !isordered(y)
        end
    end
    for x in (Nullable{Int}[1:10;], NullableArray(1:10), NullableCategoricalArray(1:10)),
        y in (similar(x), Array{Nullable{Int}}(size(x)),
              NullableArray{Int}(size(x)), NullableCategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 100, 1=>100, 2:4=>100, [5; 9:10]=>-1)
        @test y === z
        @test isequal(y, Nullable{Int}[100, 100, 100, 100, -1, 100, 100, 100, -1, -1])
        if isa(y, CatArray)
            @test levels(y) == [100, -1]
            @test !isordered(y)
        end
    end

    # Recoding from Int to Int, with default
    for x in ([1:10;], CategoricalArray(1:10)),
        y in (similar(x), Array{Int}(size(x)), CategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, -10, 1=>100, 2:4=>0, [5; 9:10]=>-1)
        @test y === z
        @test y == [100, 0, 0, 0, -1, -10, -10, -10, -1, -1]
        if isa(y, CatArray)
            @test levels(y) == [100, 0, -1, -10]
            @test !isordered(y)
        end
    end
    for x in (Nullable{Int}[1:10;], NullableArray(1:10), NullableCategoricalArray(1:10)),
        y in (similar(x), Array{Nullable{Int}}(size(x)),
              NullableArray{Int}(size(x)), NullableCategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, -10, 1=>100, 2:4=>0, [5; 9:10]=>-1)
        @test y === z
        @test isequal(y, Nullable{Int}[100, 0, 0, 0, -1, -10, -10, -10, -1, -1])
        if isa(y, CatArray)
            @test levels(y) == [100, 0, -1, -10]
            @test !isordered(y)
        end
    end

    # Recoding from Int to Int, with a first value being Float64
    for x in ([1:10;], CategoricalArray(1:10)),
        y in (similar(x), Array{Int}(size(x)), CategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 1.0=>100, 2:4=>0, [5; 9:10]=>-1)
        @test y == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
        if isa(y, CatArray)
            @test levels(y) == [6, 7, 8, 100, 0, -1]
            @test !isordered(y)
        end
    end
    for x in (Nullable{Int}[1:10;], NullableArray(1:10), NullableCategoricalArray(1:10)),
        y in (similar(x), Array{Nullable{Int}}(size(x)),
              NullableArray{Int}(size(x)), NullableCategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 1.0=>100, 2:4=>0, [5; 9:10]=>-1)
        @test isequal(y, Nullable{Int}[100, 0, 0, 0, -1, 6, 7, 8, -1, -1])
        if isa(y, CatArray)
            @test levels(y) == [6, 7, 8, 100, 0, -1]
            @test !isordered(y)
        end
    end

    # Recoding from Int to Int with overlapping pairs
    for x in ([1:10;], CategoricalArray(1:10)),
        y in (similar(x), Array{Int}(size(x)), CategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 1=>100, 2:4=>0, [5; 9:10]=>-1, 1:10=>0)
        @test y === z
        @test y == [100, 0, 0, 0, -1, 0, 0, 0, -1, -1]
        if isa(y, CatArray)
            @test levels(y) == [100, 0, -1]
            @test !isordered(y)
        end
    end
    for x in (Nullable{Int}[1:10;], NullableArray(1:10), NullableCategoricalArray(1:10)),
        y in (similar(x), Array{Nullable{Int}}(size(x)),
              NullableArray{Int}(size(x)), NullableCategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 1=>100, 2:4=>0, [5; 9:10]=>-1, 1:10=>0)
        @test y === z
        @test isequal(y, Nullable{Int}[100, 0, 0, 0, -1, 0, 0, 0, -1, -1])
        if isa(y, CatArray)
            @test levels(y) == [100, 0, -1]
            @test !isordered(y)
        end
    end

    # Recoding from Int to Int, with changes to levels order
    for x in ([1:10;], CategoricalArray(1:10)),
        y in (similar(x), Array{Int}(size(x)), CategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 1=>100, 2:4=>0, [5; 9:10]=>-1)
        @test y === z
        @test y == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
        if isa(y, CatArray)
            @test levels(y) == [6, 7, 8, 100, 0, -1]
            @test !isordered(y)
        end
    end
    for x in (Nullable{Int}[1:10;], NullableArray(1:10), NullableCategoricalArray(1:10)),
        y in (similar(x), Array{Nullable{Int}}(size(x)),
              NullableArray{Int}(size(x)), NullableCategoricalArray{Int}(size(x)), x)
        z = @inferred recode!(y, x, 1=>100, 2:4=>0, [5; 9:10]=>-1)
        @test y === z
        @test isequal(y, Nullable{Int}[100, 0, 0, 0, -1, 6, 7, 8, -1, -1])
        if isa(y, CatArray)
            @test levels(y) == [6, 7, 8, 100, 0, -1]
            @test !isordered(y)
        end
    end

    # Recoding nullable array to non-nullable categorical array: check that error is thrown
    for x in (NullableArray(["a", Nullable(), "c", "d"]),
              NullableCategoricalArray(["a", Nullable(), "c", "d"]))
        y = Vector{String}(4)
        @test_throws MethodError recode!(y, x, "a", "c"=>"b")

        y = CategoricalVector{String}(4)
        @test_throws NullException recode!(y, x, "a", "c"=>"b")
    end

    # Recoding nullable array with null value and default
    for x in (NullableArray(["a", Nullable(), "c", "d"]),
              NullableCategoricalArray(["a", Nullable(), "c", "d"])),
        y in (similar(x), Array{Nullable{String}}(size(x)),
              NullableArray{String}(size(x)), NullableCategoricalArray{String}(size(x)), x)
        z = @inferred recode!(y, x, "a", "c"=>"b")
        @test y === z
        @test isequal(y, Nullable{String}["a", Nullable(), "b", "a"])
        if isa(y, CatArray)
            @test levels(y) == ["b", "a"]
            @test !isordered(y)
        end
    end

    # Recoding nullable array with null value and no default
    for x in (NullableArray(["a", Nullable(), "c", "d"]),
              NullableCategoricalArray(["a", Nullable(), "c", "d"])),
        y in (similar(x), Array{Nullable{String}}(size(x)),
              NullableArray{String}(size(x)), NullableCategoricalArray{String}(size(x)), x)
        z = @inferred recode!(y, x, "c"=>"b")
        @test y === z
        @test isequal(y, Nullable{String}["a", Nullable(), "b", "d"])
        if isa(y, CatArray)
            @test levels(y) == ["a", "d", "b"]
            @test !isordered(y)
        end
    end

    # Recoding nullable array with null value, no default and with null as a key pair
    for x in (NullableArray(["a", Nullable(), "c", "d"]),
              NullableCategoricalArray(["a", Nullable(), "c", "d"])),
        y in (similar(x), Array{Nullable{String}}(size(x)),
              NullableArray{String}(size(x)), NullableCategoricalArray{String}(size(x)), x)
        z = @inferred recode!(y, x, "a", "c"=>"b", Nullable()=>"d")
        @test y === z
        @test isequal(y, Nullable{String}["a", "d", "b", "a"])
        if isa(y, CatArray)
            @test levels(y) == ["b", "d", "a"]
            @test !isordered(y)
        end
    end

    # Recoding nullable array with null value, no default and with null as a key pair
    for x in (NullableArray(["a", Nullable(), "c", "d"]),
              NullableCategoricalArray(["a", Nullable(), "c", "d"])),
        y in (similar(x), Array{Nullable{String}}(size(x)),
              NullableArray{String}(size(x)), NullableCategoricalArray{String}(size(x)), x)
        z = @inferred recode!(y, x, "c"=>"b", Nullable()=>"d")
        @test y === z
        @test isequal(y, Nullable{String}["a", "d", "b", "d"])
        if isa(y, CatArray)
            @test levels(y) == ["a", "b", "d"]
            @test !isordered(y)
        end
    end

    # Recoding into array with incompatible size
    for x in (NullableArray(["a", Nullable(), "c", "d"]),
              NullableCategoricalArray(["a", Nullable(), "c", "d"])),
        y in (similar(x, 0), Array{Nullable{String}}(0),
              NullableArray{String}(0), NullableCategoricalArray{String}(0))
        @test_throws DimensionMismatch recode!(y, x, "c"=>"b", Nullable()=>"d")
    end

    # Recoding into array with incompatible element type
    for x in ([1:10;], CategoricalArray(1:10)),
        y in (similar(x, String), Array{String}(size(x)), CategoricalArray{String}(size(x)))
        @test_throws ArgumentError recode!(y, x, 1=>"a", 2:4=>"b", [5; 9:10]=>"c")
    end
    for x in (Nullable{Int}[1:10;], NullableArray(1:10), NullableCategoricalArray(1:10)),
        y in (similar(x), Array{Nullable{Int}}(size(x)),
              NullableArray{Int}(size(x)), NullableCategoricalArray{Int}(size(x)))
        res = @test_throws MethodError recode!(y, x, 1=>"a", 2:4=>"b", [5; 9:10]=>"c")
    end

    ## Test in-place recode!()

    # Recoding from Int to Int without default
    for x in ([1:10;], CategoricalArray(1:10))
        z = @inferred recode!(x, 1=>100, 2:4=>0, [5; 9:10]=>-1)
        @test x === z
        @test x == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
        if isa(x, CatArray)
            @test levels(x) == [6, 7, 8, 100, 0, -1]
            @test !isordered(x)
        end
    end
    for x in (Nullable{Int}[1:10;], NullableArray(1:10), NullableCategoricalArray(1:10))
        z = @inferred recode!(x, 1=>100, 2:4=>0, [5; 9:10]=>-1)
        @test x === z
        @test isequal(x, Nullable{Int}[100, 0, 0, 0, -1, 6, 7, 8, -1, -1])
        if isa(x, CatArray)
            @test levels(x) == [6, 7, 8, 100, 0, -1]
            @test !isordered(x)
        end
    end

    # Recoding from Int to Int with default
    for x in ([1:10;], CategoricalArray(1:10))
        z = @inferred recode!(x, 1, 1=>100, 2:4=>0, [5; 9:10]=>-1)
        @test x === z
        @test x == [100, 0, 0, 0, -1, 1, 1, 1, -1, -1]
        if isa(x, CatArray)
            @test levels(x) == [100, 0, -1, 1]
            @test !isordered(x)
        end
    end
    for x in (Nullable{Int}[1:10;], NullableArray(1:10), NullableCategoricalArray(1:10))
        z = @inferred recode!(x, 1, 1=>100, 2:4=>0, [5; 9:10]=>-1)
        @test x === z
        @test isequal(x, Nullable{Int}[100, 0, 0, 0, -1, 1, 1, 1, -1, -1])
        if isa(x, CatArray)
            @test levels(x) == [100, 0, -1, 1]
            @test !isordered(x)
        end
    end

    ## Test recode() promotion

    for x in (1:10, [1:10;], CategoricalArray(1:10),
              NullableArray(1:10), NullableCategoricalArray(1:10))
        # Recoding from Int to Float64 due to a second value being Float64
        y = @inferred recode(x, 1=>100.0, 2:4=>0, [5; 9:10]=>-1)
        if eltype(x) <: Nullable
            @test isequal(y, Nullable{Int}[100, 0, 0, 0, -1, 6, 7, 8, -1, -1])
        else
            @test y == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
        end
        if isa(x, CatArray)
            @test isa(y, CatVector{Float64, DefaultRefType})
            @test levels(y) == [6, 7, 8, 100, 0, -1]
            @test !isordered(y)
        elseif eltype(x) <: Nullable
            @test typeof(y) === NullableVector{Float64}
        else
            @test typeof(y) === Vector{Float64}
        end

        # Recoding from Int to Float64, with Float64 default and all other values Int
        y = @inferred recode(x, -10.0, 1=>100, 2:4=>0, [5; 9:10]=>-1)
        if eltype(x) <: Nullable
            @test isequal(y, Nullable{Int}[100, 0, 0, 0, -1, -10, -10, -10, -1, -1])
        else
            @test y == [100, 0, 0, 0, -1, -10, -10, -10, -1, -1]
        end
        if isa(x, CatArray)
            @test isa(y, CatVector{Float64, DefaultRefType})
            @test levels(y) == [100, 0, -1, -10]
            @test !isordered(y)
        elseif eltype(x) <: Nullable
            @test typeof(y) === NullableVector{Float64}
        else
            @test typeof(y) === Vector{Float64}
        end

        # Recoding from Int to Any
        y = @inferred recode(x, 1=>"a", 2:4=>0, [5; 9:10]=>-1)
        if eltype(x) <: Nullable
            @test isequal(y, Nullable{Any}["a", 0, 0, 0, -1, 6, 7, 8, -1, -1])
        else
            @test y == ["a", 0, 0, 0, -1, 6, 7, 8, -1, -1]
        end
        if isa(x, CatArray)
            @test isa(y, CatVector{Any, DefaultRefType})
            @test levels(y) == [6, 7, 8, "a", 0, -1]
            @test !isordered(y)
        elseif eltype(x) <: Nullable
            @test typeof(y) === NullableVector{Any}
        else
            @test typeof(y) === Vector{Any}
        end

        # Recoding from Int to String, with String default
        y = @inferred recode(x, "d", 1=>"a", 2:4=>"b", [5; 9:10]=>"c")
        if eltype(x) <: Nullable
            @test isequal(y, Nullable{String}["a", "b", "b", "b", "c", "d", "d", "d", "c", "c"])
        else
            @test y == ["a", "b", "b", "b", "c", "d", "d", "d", "c", "c"]
        end
        if isa(x, CatArray)
            @test isa(y, CatVector{String, DefaultRefType})
            @test levels(y) == ["a", "b", "c", "d"]
            @test !isordered(y)
        elseif eltype(x) <: Nullable
            @test typeof(y) === NullableVector{String}
        else
            @test typeof(y) === Vector{String}
        end

        # Recoding from Int to String, with all original levels recoded
        y = @inferred recode(x, 1:4=>"a", [5; 9:10]=>"b", 6:8=>"c")
        if eltype(x) <: Nullable
            @test isequal(y, Nullable{String}["a", "a", "a", "a", "b", "c", "c", "c", "b", "b"])
        else
            @test y == ["a", "a", "a", "a", "b", "c", "c", "c", "b", "b"]
        end
        if isa(x, CatArray)
            @test isa(y, CatVector{String, DefaultRefType})
            @test levels(y) == ["a", "b", "c"]
            @test !isordered(y)
        elseif eltype(x) <: Nullable
            @test typeof(y) === NullableVector{String}
        else
            @test typeof(y) === Vector{String}
        end

        # Recoding from Int to Int/String (i.e. Any), with default String and other values Int
        y = @inferred recode(x, "x", 1=>100, 2:4=>0, [5; 9:10]=>-1)
        if eltype(x) <: Nullable
            @test isequal(y, Nullable{Any}[100, 0, 0, 0, -1, "x", "x", "x", -1, -1])
        else
            @test y == [100, 0, 0, 0, -1, "x", "x", "x", -1, -1]
        end
        if isa(x, CatArray)
            @test isa(y, CatVector{Any, DefaultRefType})
            @test levels(y) == [100, 0, -1, "x"]
            @test !isordered(y)
        elseif eltype(x) <: Nullable
            @test typeof(y) === NullableVector{Any}
        else
            @test typeof(y) === Vector{Any}
        end

        # Recoding from Int to Int/String, without any Int value in pairs
        # and keeping some original Int levels
        # This must fail since we cannot take into account original eltype (whether original
        # levels are kept is only known at run time)
        res = @test_throws ArgumentError recode(x, 1=>"a", 2:4=>"b", [5; 9:10]=>"c")
        @test sprint(showerror, res.value) ==
            "ArgumentError: cannot `convert` value 6 (of type $Int) to type of recoded levels (String). " *
            "This will happen with recode() when not all original levels are recoded " *
            "(i.e. some are preserved) and their type is incompatible with that of recoded levels."

        # Recoding from Int to Nullable{Int} with null default
        y = @inferred recode(x, Nullable(), 1=>100, 2:4=>0, [5; 9:10]=>-1)
        @test isequal(y, [100, 0, 0, 0, -1, Nullable(), Nullable(), Nullable(), -1, -1])
        if isa(x, CatArray)
            @test isa(y, NullableCategoricalVector{Int, DefaultRefType})
            @test levels(y) == [100, 0, -1]
            @test !isordered(y)
        else
            @test typeof(y) === NullableVector{Int}
        end

        # Recoding from Int to Nullable{Int} with nullable non-null default
        y = @inferred recode(x, Nullable(1), 1=>100, 2:4=>0, [5; 9:10]=>-1)
        @test isequal(y, Nullable{Int}[100, 0, 0, 0, -1, 1, 1, 1, -1, -1])
        if isa(x, CatArray)
            @test isa(y, NullableCategoricalVector{Int, DefaultRefType})
            @test levels(y) == [100, 0, -1, 1]
            @test !isordered(y)
        else
            @test typeof(y) === NullableVector{Int}
        end

        # Recoding from Int to Nullable{Int} with null RHS
        y = @inferred recode(x, 1=>Nullable(), 2:4=>0, [5; 9:10]=>-1)
        @test isequal(y, [Nullable(), 0, 0, 0, -1, 6, 7, 8, -1, -1])
        if isa(x, CatArray)
            @test isa(y, NullableCategoricalVector{Int, DefaultRefType})
            @test levels(y) == [6, 7, 8, 0, -1]
            @test !isordered(y)
        else
            @test typeof(y) === NullableVector{Int}
        end
    end

    for x in (["a", "c", "b", "a"], CategoricalArray(["a", "c", "b", "a"]))
        # Recoding from String to String
        y = @inferred recode(x, "c"=>"x", "b"=>"y", "a"=>"z")
        @test y == ["z", "x", "y", "z"]
        if isa(x, CatArray)
            @test typeof(y) === CategoricalVector{String, DefaultRefType}
            @test levels(y) == ["x", "y", "z"]
            @test !isordered(y)
        else
            @test typeof(y) === Vector{String}
        end
    end

    for x in (['a' 'c'; 'b' 'a'], CategoricalArray(['a' 'c'; 'b' 'a']))
        # Recoding a Matrix
        y = @inferred recode(x, 'c'=>'x', 'b'=>'y', 'a'=>'z')
        @test y == ['z' 'x'; 'y' 'z']
        if isa(x, CatArray)
            @test typeof(y) === CategoricalMatrix{Char, DefaultRefType}
            @test levels(y) == ['x', 'y', 'z']
            @test !isordered(y)
        else
            @test typeof(y) === Matrix{Char}
        end
    end

    for x in (10:-1:1, CategoricalArray(10:-1:1))
        # Recoding from Int to Int/String (i.e. Any), with index and levels in different orders
        y = @inferred recode(x, 0, 1=>"a", 2:4=>"c", [5; 9:10]=>"b")
        @test y == ["b", "b", 0, 0, 0, "b", "c", "c", "c", "a"]
        if isa(x, CatArray)
            @test typeof(y) === CategoricalVector{Any, DefaultRefType}
            @test levels(y) == ["a", "c", "b", 0]
            @test !isordered(y)
        else
            @test typeof(y) === Vector{Any}
        end

        # Recoding from Int to String via default, with index and levels in different orders
        y = @inferred recode(x, "x", 1=>"a", 2:4=>"c", [5; 9:10]=>"b")
        @test y == ["b", "b", "x", "x", "x", "b", "c", "c", "c", "a"]
        if isa(x, CatArray)
            @test typeof(y) === CategoricalVector{String, DefaultRefType}
            @test levels(y) == ["a", "c", "b", "x"]
            @test !isordered(y)
        else
            @test typeof(y) === Vector{String}
        end
    end

    # Recoding CategoricalArray with custom reftype
    x = CategoricalVector{Int, UInt8}(1:10)
    y = @inferred recode(x, 1=>100, 2:4=>0, [5; 9:10]=>-1)
    @test y == [100, 0, 0, 0, -1, 6, 7, 8, -1, -1]
    @test typeof(y) === CategoricalVector{Int, UInt8}
    @test levels(y) == [6, 7, 8, 100, 0, -1]
    @test !isordered(y)

    # Recoding ordered CategoricalArray and merging two categories
    x = CategoricalArray(["a", "c", "b", "a"])
    ordered!(x, true)
    y = @inferred recode(x, "c"=>"a")
    @test y == ["a", "a", "b", "a"]
    @test typeof(y) === CategoricalVector{String, DefaultRefType}
    @test levels(y) == ["a", "b"]
    @test isordered(y)

    # Recoding ordered CategoricalArray and merging one category into another (contd.)
    x = CategoricalArray(["a", "c", "b", "a"])
    ordered!(x, true)
    y = @inferred recode(x, "b"=>"c")
    @test y == ["a", "c", "c", "a"]
    @test typeof(y) === CategoricalVector{String, DefaultRefType}
    @test levels(y) == ["a", "c"]
    @test isordered(y)

    # Recoding ordered CategoricalArray with new level which cannot be ordered
    x = CategoricalArray(["a", "c", "b", "a"])
    ordered!(x, true)
    y = @inferred recode(x, "b"=>"d")
    @test y == ["a", "c", "d", "a"]
    @test typeof(y) === CategoricalVector{String, DefaultRefType}
    @test levels(y) == ["a", "c", "d"]
    @test !isordered(y)

    # Recoding ordered CategoricalArray with conflicting orders
    x = CategoricalArray(["a", "c", "b", "a"])
    ordered!(x, true)
    y = @inferred recode(x, "b"=>"b", "a"=>"a")
    @test y == ["a", "c", "b", "a"]
    @test typeof(y) === CategoricalVector{String, DefaultRefType}
    @test levels(y) == ["c", "b", "a"]
    @test !isordered(y)

    # Recoding ordered CategoricalArray with default already in levels
    x = CategoricalArray(["a", "c", "b", "a"])
    ordered!(x, true)
    y = @inferred recode(x, "a", "c"=>"b")
    @test y == ["a", "b", "a", "a"]
    @test typeof(y) === CategoricalVector{String, DefaultRefType}
    @test levels(y) == ["b", "a"]
    @test !isordered(y)

    # Recoding ordered CategoricalArray with default not in levels
    x = CategoricalArray(["d", "c", "b", "d"])
    ordered!(x, true)
    y = @inferred recode(x, "a", "c"=>"b")
    @test y == ["a", "b", "a", "a"]
    @test typeof(y) === CategoricalVector{String, DefaultRefType}
    @test levels(y) == ["b", "a"]
    @test !isordered(y)

    # Recoding NullableCategoricalArray with null values and no default
    x = NullableCategoricalArray(["a", "b", "c", "d"])
    x[2] = Nullable()
    y = @inferred recode(x, "c"=>"b")
    @test isequal(y, ["a", Nullable(), "b", "d"])
    @test typeof(y) === NullableCategoricalVector{String, DefaultRefType}
    @test levels(y) == ["a", "b", "d"]
    @test !isordered(y)

    # Recoding NullableCategoricalArray with null values and non-null default
    x = NullableCategoricalArray(["a", "b", "c", "d"])
    x[2] = Nullable()
    y = @inferred recode(x, "a", "c"=>"b")
    @test isequal(y, ["a", Nullable(), "b", "a"])
    @test typeof(y) === NullableCategoricalVector{String, DefaultRefType}
    @test levels(y) == ["b", "a"]
    @test !isordered(y)
end


text:
struct DummyThermostat <: AbstractThermostat; end

# TODO need to be able to thermostat DoFs individually!!!!
"""
Generic Langevin thermostat for a single particle
"""
mutable struct Langevin <: AbstractLangevin
  T::Float64
  gamma::Float64
  scale::Float64
  seed::UInt64
  rng::AbstractRNG
  rands::Vector{Float64}
  function Langevin(T::Float64,gamma::Float64,seed::UInt64)
    new(T,gamma,sqrt(2.0*T/gamma),seed,MersenneTwister(seed),[])
  end
end
Langevin(T::Float64,gamma::Float64) = Langevin(T,gamma,rand(UInt64))

function applythermostat!(dt::Float64,sys::AbstractThermostattedSystem,langevin::Langevin)
  sys.x .+= (sqrt(dt) * langevin.scale) .* langevin.rands
end

"""
Langevin thermostat for a single particle with different frictions per DoF
"""
mutable struct LangevinND <: AbstractLangevin
  T::Float64
  gamma::Vector{Float64}
  scale::Vector{Float64}
  seed::UInt64
  rng::AbstractRNG
  rands::Vector{Float64}
  function LangevinND(T::Float64,gamma::Vector{Float64},seed::UInt64)
    new(T,gamma,sqrt.(2*T ./ gamma),seed,MersenneTwister(seed),[])
  end
end
LangevinND(T::Float64,gamma::Vector{Float64}) = LangevinND(T,gamma,rand(UInt64))

function applythermostat!(dt::Float64,sys::AbstractThermostattedSystem,langevin::LangevinND)
  sys.x .+= (sqrt(dt) .* langevin.scale .* langevin.rands)
end

"""
Langevin thermostat for a single active particle
"""
mutable struct ActiveLangevin <: AbstractLangevin
  Tt::Float64
  gammat::Float64
  Tr::Float64
  gammar::Float64
  scale::Vector{Float64}
  seed::UInt64
  rng::AbstractRNG
  rands::Vector{Float64}
  function ActiveLangevin(dim::Int64,Tt::Float64,gammat::Float64,Tr::Float64,gammar::Float64,seed::UInt64)
    scale = zeros(Float64,dim)
    scale[1:dim] .= sqrt(2.0*Tt/gammat)
    scale[end] = sqrt(2.0*Tr/gammar)
    new(Tt,gammat,Tr,gammar,scale,seed,MersenneTwister(seed),[])
  end
end
ActiveLangevin(dim::Int64,Tt::Float64,gammat::Float64,Tr::Float64,gammar::Float64) = ActiveLangevin(dim,Tt,gammat,Tr,gammar,rand(UInt64))

function applythermostat!(dt::Float64,sys::AbstractActiveSystem,langevin::ActiveLangevin)
  @. @views sys.x[1:sys.dim-1] .+= (dt*sys.v0) .* [cos(sys.x[end]),sin(sys.x[end])]
  sys.x .+= (sqrt(dt) .* langevin.scale .* langevin.rands)
end

"""
Functions for access to the random number generator
"""
function getseed(langevin::AbstractLangevin)
  return langevin.seed
end

function setseed!(seed::Union{Int64,UInt64},langevin::AbstractLangevin)
  langevin.seed = convert(UInt64,seed)
  langevin.rng = MersenneTwister(langevin.seed)
end

"""
Various functions for generating noise
"""
function gennoises(nsteps::Int64,dim::Int64,rng::AbstractRNG)
  randn(rng,Float64,(nsteps,dim))
end

function gennoises!(nsteps::Int64,dim::Int64,langevin::AbstractLangevin)
  langevin.rands = randn(langevin.rng,Float64,(nsteps,dim))
end

function gennoise!(dim::Int64,langevin::AbstractLangevin)
  langevin.rands = randn(langevin.rng,Float64,dim)
  #langevin.rands .= 1.0
end

#function gennoise!(dim::Int64,langevin::AbstractLangevin)
#  langevin.rands .= 1.0
#  #r1::Float64 = rand(Float64)
#  #r2::Float64 = rand(Float64)
#  #theta::Float64 = 2.0*4.0*atan(1.0)*r1
#  #r::Float64 = sqrt(-2.0*log(r2))
#  #langevin.rands[1] = 1.0#r*cos(theta)
#  #langevin.rands[2] = 1.0#r*sin(theta)
#end

"""
"""
function diffusion!(du::Array{Float64},langevin::Langevin)
  du .= 1.0
end

ext: jl


text:
"""
	Checks if elements of vec1 are contained in vec2 in the same order (possibly with elements in between)
"""
function containedin(vec1::Vector, vec2::Vector)
	max_elements = length(vec1)
	vec1_index = 1 # keeps track where we are in the firsst vector
	for item in vec2
		if vec1_index > max_elements
			return true
		end

		if item == vec1[vec1_index]
			vec1_index += 1  # increase the index every time we encounter the matching element
		end
	end

	return vec1_index > max_elements
end


text:
@inline function policy(m::KrusellSmith,
                        net_return::Float64, wage::Float64,
                        l_in::Vector{Float64}, g::Vector{Float64}, KF::Array{Float64,2};
                        damp::Float64 = 0.8, tol::Float64 = 1e-5)
    @assert 0 <= damp <= 1

    nw::Int64                = get_setting(m, :nw)
    ns::Int64                = get_setting(m, :ns)
    wgrid::DSGE.Grid         = get_grid(m, :wgrid)
    sgrid::DSGE.Grid         = get_grid(m, :sgrid)

    zhi::Float64             = get_setting(m, :zhi)
    zlo::Float64             = get_setting(m, :zlo)
    smoother::Float64        = get_setting(m, :smoother)
    γ::Float64               = m[:γ]
    β::Float64               = m[:β]

    dist::Float64            = 1.0         # distance between l_in and l_out on decision rules
    l_out::Array{Float64,1}  = copy(l_in)  # initialize l_out

    ap = Array{Float64,1}(undef, nw) # savings
    q  = Array{Float64,2}(undef, nw, nw * ns) # auxilary variable to compute LHS of euler

    count = 1

    while dist > tol && count < 5000

        # compute c(w) given guess for l_in = β*R*E[u'(c_{t+1})]
        c::Array{Float64,1} = min.(l_in.^(-1 / γ), wgrid.points)

        # compute ap(w) given guess for l_in
        ap = wgrid.points - c

        # The term contained inside the (⋅) next to g in equation (9)
        q = euler_lhs(q, ap, wgrid.points, sgrid.points, zhi, zlo, smoother, wage, net_return, nw, ns)

        # Equation (9)
        l_out = β * net_return * q * kron(sgrid.weights .* (g / wage),
                                              wgrid.weights .* (c.^-γ))
        dist = maximum(abs.(l_out - l_in))
        l_in = damp * l_out + (1.0 - damp) * l_in;

        count += 1
    end
    if count == 5000
        @warn "Euler iteration did not converge"
    end

    # Equation (10)
    # You can re-use most of the functional form of the Euler equation to
    # back out the KF for μ(w)
    KF = kolmogorov_forward(q, KF, g, sgrid.weights, wage, nw, ns)

    return ap, l_out, KF
end

function steadystate!(m::KrusellSmith, tol::Float64 = 1e-5, maxit::Int64 = 100, Kdamp::Float64 = 0.1)
    # Load settings and grids
    nw::Int64  = get_setting(m, :nw)

    wgrid::DSGE.Grid = get_grid(m, :wgrid)
    sgrid::DSGE.Grid = get_grid(m, :sgrid)

    zhi::Float64      = get_setting(m, :zhi)
    zlo::Float64      = get_setting(m, :zlo)
    smoother::Float64 = get_setting(m, :smoother)

    α::Float64 = m[:α]
    δ::Float64 = m[:δ]
    β::Float64 = m[:β]

    g = trunc_lognpdf(sgrid, m[:μ_s].value, m[:σ_s].value)      # Density g (of skill, s) evaluated on sgrid

    m[:Lstar] = L::Float64 = quadrature_sum(g, sgrid)           # Aggregate Labor endowment

    K::Float64 = L*(α /((1.0 / β) - 1.0 + δ))^(1.0 / (1.0 - α)) # Aggregate Capital Stock

    l_in = fill(0.3, nw)        # This corresponds to l^i_t = βE_t[R_{t+1} u'(c^i_{t+1})]
                                #                           = βE_t[R_{t+1} min(l^i_{t+1}, (w^i_{t+1}^(-1/γ))]

    ap   = Array{Float64,1}(undef, nw)
    KF   = Array{Float64,2}(undef, nw, nw) # The Kolmogorov Forward Equation for μ
    μ    = Array{Float64,1}(undef, nw)     # The cross-sectional density of cash on hand, w
    newK = 0.0

    diseqm = 1000.
    count = 1

    # Equation (11)
    while abs(diseqm) > tol && count < maxit # clearing markets

        wage = MPL(1.0, L, K, α)
        net_return = MPK(1.0, L, K, α, δ)

        ap, l_in, KF = policy(m, net_return, wage, l_in, g, KF)

        # sspolicy returns the consumption decision rule, a prime
        # decision rule, l_in = updated β R u'(c_{t+1}),
        # KF is the Kolmogorov forward operator and is the map which moves you
        # from cash on hand distribution today to cash on hand dist tomorrow

        # find eigenvalue closest to 1
        D,V = eigen(wgrid.weights[1] * KF)

        if norm(D[1] - 1.0) > 2e-1 # that's the tolerance we are allowing
            @warn "your eigenvalue is too far from 1, something is wrong"
        end
        μ = real(@view(V[:,1]))  # Pick the eigenvector associated with the largest
                                 # eigenvalue and moving it back to values

        μ = μ / (wgrid.weights' * μ)        # Scale of eigenvectors not determinate:
                                            # rescale to integrate to exactly 1
        newK  = wgrid.weights' * (μ .* ap)  # compute excess supply of savings, which is a fn of w
        diseqm = newK - K

        K += Kdamp * diseqm
        count += 1
    end

    m[:lstar] = l_in
    m[:cstar] = min.(m[:lstar].value.^(-1.0/m[:γ]), wgrid.points)
    m[:μstar] = μ
    m[:Kstar] = newK
    m[:KFstar] = KF

    nothing
end

@inline function mollifier(z::Float64,zhi::Float64,
                   zlo::Float64, smoother::Float64; In::Float64=0.443993816237631)
    # mollifier function
    if z < zhi && z > zlo
        temp = (-1.0 + 2.0 * (z - zlo) / (zhi - zlo)) / smoother
        out = (2.0/(zhi - zlo)) * exp(-1.0 / (1.0 - temp*temp)) / (In * smoother);
    else
        out = 0.0
    end
    return out
end

# Derivative of the mollifier
@inline function dmollifier(x::Float64, zhi::Float64, zlo::Float64, smoother::Float64;
                            In::Float64=0.443993816237631)
    temp = (-1.0 + 2.0*(x-zlo)/(zhi-zlo))/smoother
    dy  = -(2*temp./((1 - temp.^2).^2)).*(2/(smoother*(zhi-zlo))).*mollifier(x, zhi, zlo, smoother)
end

@inline function trunc_lognpdf(grid::Grid, mu::Float64, sig::Float64)
    x = grid.points
    xlo = x[1]
    xhi = x[end]

    # truncated log normal pdf
    logn = LogNormal(mu, sig)
    return (x .- xlo .>= -eps()) .* (x .- xhi .<= eps()) .* pdf.(logn,x) / (cdf(logn,xhi) .- cdf(logn,xlo))
end

# ω_t
@inline function MPL(Z::Float64, L::Float64, K::Float64, α::Float64)
    (1.0 - α) * Z * (K^α) * (L^(-α))  # Marginal product of labor
end

# r_t + (1 - δ)
@inline function MPK(Z::Float64, L::Float64, K::Float64, α::Float64, δ::Float64)
    α * Z * (K^(α - 1.0)) * (L^(1.0 - α)) + (1.0 - δ)   # Marginal product of capital + 1-δ
end

####### Is this function ever used? #########
#=
# Re-computing
@inline function kolmogorov_forward(g::Vector{Float64}, sweights::Vector{Float64},
                            ap::Vector{Float64}, wgrid::Vector{Float64}, sgrid::Vector{Float64},
                            zhi::Float64, zlo::Float64, smoother::Float64,
                            wage::Float64, R::Float64, nw::Int64, ns::Int64)

    KF = Array{Float64, 2}(nw, nw)   # transition matrix mapping todays dist of w to w'
    for iw=1:nw
        for iwp=1:nw
            sumns = 0.0
            cursum = (wgrid[iwp] - R*ap[iw]) / wage
            for isp=1:ns
                sumns += mollifier(cursum - sgrid[isp], zhi, zlo, smoother) * (g[isp]) * sweights[isp]
            end
            @inbounds KF[iwp,iw] = sumns / wage
        end
    end
    return KF
end
=#

# Indexing into q since its already been computed
@inline function kolmogorov_forward(q::Matrix{Float64}, KF::Array{Float64,2}, g::Vector{Float64},
                                    sweights::Vector{Float64}, wage::Float64, nw::Int64, ns::Int64)
    g_sweights_prod::Array{Float64,1} = g / wage .* sweights
    nw_array = nw * (0:ns-1)
    for iw = 1:nw
        for iwp = 1:nw
            @inbounds KF[iwp, iw] = sum(q[iw, nw_array .+ iwp] .* g_sweights_prod)
        end
    end
    return KF
end

@inline function euler_lhs(q::Matrix{Float64}, ap::Vector{Float64}, wgrid::Vector{Float64},
                   sgrid::Vector{Float64}, zhi::Float64, zlo::Float64, smoother::Float64,
                   wage::Float64, net_return::Float64, nw::Int64, ns::Int64)
    for iwp = 1:nw
        for iw = 1:nw
            wgrid_term = (wgrid[iwp] - net_return * ap[iw]) / wage
            for iss = 1:ns
                @inbounds q[iw, nw * (iss - 1) + iwp] = mollifier(wgrid_term - sgrid[iss], zhi, zlo, smoother)
            end
        end
    end
    return q
end


text:
using .Ranger.Custom

mutable struct GameScene <: Ranger.AbstractScene
    base::Nodes.NodeData
    transform::Nodes.TransformProperties{Float64}

    replacement::Ranger.AbstractScene

    # Collection of layers. The first is always a background layer
    children::Array{Ranger.AbstractNode,1}

    function GameScene(world::Ranger.World, name::String)
        obj = new()

        obj.base = Nodes.NodeData(name, Nodes.NodeNil(), world)
        obj.replacement = Nodes.SceneNil()
        obj.transform = Nodes.TransformProperties{Float64}()
        obj.children = Array{Ranger.AbstractNode,1}[]

        obj
    end
end

function build(scene::GameScene, world::Ranger.World)
    layer = GameLayer(world, "GameLayer", scene)
    push!(scene.children, layer)
    build(layer, world);

    cross = Custom.CrossNode(world, "CrossNode", scene)
    Nodes.set_nonuniform_scale!(cross.transform, Float64(world.view_width), Float64(world.view_height))
    cross.color = Rendering.LightGray()
    push!(scene.children, cross)

    text = Custom.VectorTextNode(world, "VectorTextNode", scene)
    Nodes.set_scale!(text, 25.0)
    # Nodes.set_rotation_in_degrees!(text, 45.0)
    text.color = RangerGame.orange
    Custom.set_text!(text, "RANGER IS A GO!")
    # Custom.set_text!(text, "=,./!:_-")
    push!(scene.children, text)

    # Note: If you apply a scale that is the size of the view-space
    # then the rest of your code will need to define all mesh data
    # in unit-space.
    # set_nonuniform_scale!(scene, world.view_width, world.view_height)
    # # Bake in the transform rather repeatedly perform in draw()
    # calc_transform!(scene.transform)

    # # This node/scene is never dirtied at the scene level
    Nodes.set_dirty!(scene, false);
end

function Nodes.transition(scene::GameScene)
    Scenes.NO_ACTION
end

function Nodes.get_replacement(scene::GameScene)
    scene.replacement
end

# --------------------------------------------------------
# Grouping
# --------------------------------------------------------
function Nodes.get_children(node::GameScene)
    node.children
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(pager),String})
end


text:
using CapacityExpansion
using Clp
using Plots
optimizer=Clp.Optimizer
state="GER_1"
ts_input_data = load_timeseries_data_provided(state; T=24, years=[2016])

plot(ts_input_data.data["solar-germany"], legend=false, linestyle=:dot, xlabel="Time [h]", ylabel="Solar availability factor [%]")
savefig(joinpath(dirname(@__FILE__),"preparing_clust_data_load.svg"))

ts_clust_data = run_clust(ts_input_data;method="kmeans",representation="centroid",n_init=50,n_clust=5).best_results
plot(ts_clust_data.data["solar-germany"], legend=false, linestyle=:solid, width=3, xlabel="Time [h]", ylabel="Solar availability factor [%]")
savefig(joinpath(dirname(@__FILE__),"preparing_clust_data_agg.svg"))

cep_data = load_cep_data_provided(state)

co2_result = run_opt(ts_clust_data,cep_data,optimizer;descriptor="co2",co2_limit=500) #hide

# use the get variable set in order to get the labels: indicate the variable as "CAP" and the set-number as 1 to receive those set values
variable=co2_result.variables["CAP"]
labels=axes(variable,"tech")

data=variable[:,:,"germany"]
# use the data provided for a simple bar-plot without a legend
bar(data,title="Cap", xticks=(1:length(labels),labels),legend=false, ylabel="Capacity [MW]", xlabel="technologies", color="orange")
savefig(joinpath(dirname(@__FILE__),"opt_cep_cap_plot.svg"))


text:
######################################################################
# includes stuff I don't want to need distributions for...
######################################################################

# only handles stuff in 2d
function my_pdf(x::Vector{Float64}, m::Vector{Float64}, S::Matrix{Float64})
	return my_pdf((x[1],x[2]), m, S)
end

function my_pdf(x::NTuple{2,Float64},m::Vector{Float64},S::Matrix{Float64})
	dx = x[1] - m[1]
	dy = x[2] - m[2]

	Si = inv(S)
	e_stuff = dx*Si[1,1]*dx + dx*Si[1,2]*dy + dy*Si[2,1]*dx + dy*Si[2,2]*dy
	return det(sqrtm(Si)) * exp(-0.5 * e_stuff) / (2.0 * pi)
end

function my_cdf(sigma::Float64, x::Float64)
	temp = x / (sigma*sqrt(2.0))
	return 0.5 * (1. + erf(temp))
end

# Fits an angle into -180 to 180
# Assume the angle is within -360 to 360
function fit_180(angle::Float64)
	if angle > 180.0
		#angle = 360.0 - angle   # omfg this is how i did it before. wrong.
		angle -= 360.0
	elseif angle < -180.0
		angle += 360.0
	end
	return angle
end

######################################################################
# Find the true angle between UAV and jammer.
#
# Parameters:
#  xp - location of vehicle
#  theta - location of jammer
#
# Returns true angle, measured from north, in degrees.
function true_bearing(xp::NTuple{2,Float64}, theta)
	xr = theta[1] - xp[1]
	yr = theta[2] - xp[2]
	return mod(rad2deg(atan2(xr,yr)), 360)
end
true_bearing(p::Pose, theta) = true_bearing( (p[1], p[2]), theta)
function true_bearing(xp, theta::Vector{Float64})
	return true_bearing(xp, (theta[1], theta[2]))
end
function true_bearing(xp::Pose, theta::Vector{Float64})
	return true_bearing((xp[1], xp[2]), (theta[1], theta[2]))
end

export get_distance
function get_distance(p1, p2)
    dx = p1[1] - p2[1]
    dy = p1[2] - p2[2]

    return sqrt(dx*dx + dy*dy)
end
function get_distance2(p1, p2)
    dx = p1[1] - p2[1]
    dy = p1[2] - p2[2]

    return dx*dx + dy*dy
end

######################################################################
# I don't remember what these are for at all
function angle_mean(v)
	n = length(v)
	sin_sum = 0.0
	cos_sum = 0.0
	for i = 1:n
		sin_sum += sin(v[i])
		cos_sum += cos(v[i])
	end
	return mod(rad2deg(atan2(sin_sum, cos_sum)), 360.)
end

function angle_mean(v,w)
	n = length(v)
	sin_sum = 0.0
	cos_sum = 0.0
	for i = 1:n
		sin_sum += w[i] * sin(v[i])
		cos_sum += w[i] * cos(v[i])
	end
	return mod(rad2deg(atan2(sin_sum, cos_sum)), 360.)
end

######################################################################
# I don't really remember using this
# Takes a vector of 36 angles and turns it into 360 via linear interpolation
function makenorm()
	data = readcsv("norm.csv")
	values = vec(data[:,2])
	new_values = zeros(360)
	for i = 1:35
		val_start = 10*(i-1)+1
		new_values[val_start] = values[i]
		for j = 1:9
			new_values[val_start + j] = (10.0-j)*values[i] + j*values[i+1]
			new_values[val_start + j] /= 10.0
		end
	end
	# now handle the last one
	i = 36
	val_start = 10*(i-1)+1
	new_values[val_start] = values[i]
	for j = 1:9
		new_values[val_start + j] = (10.0-j)*values[i] + j*values[1]
		new_values[val_start + j] /= 10.0
	end
	return new_values
end


text:
"""
### Cutout sub-regions of the data base of DataSetType
- select shape of a region
- select size of a region (w/o intersecting cells)
- give the spatial center (with units) of the data relative to the full box
- relate the coordinates to a direction (x,y,z)
- inverse the selected region

```julia
subregion(dataobject::DataSetType, shape::Symbol=:cuboid;
            xrange::Array{<:Any,1}=[missing, missing],  # cuboid
            yrange::Array{<:Any,1}=[missing, missing],  # cuboid
            zrange::Array{<:Any,1}=[missing, missing],  # cuboid

            radius::Real=0.,              # cylinder, sphere
            height::Real=0.,              # cylinder
            direction::Symbol=:z,           # cylinder

            center::Array{<:Any,1}=[0.,0.,0.],   # all
            range_unit::Symbol=:standard,          # all
            cell::Bool=true,                        # hydro and gravity
            inverse::Bool=false,                    # all
            verbose::Bool=verbose_mode)             # all
```

#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "DataSetType"
- **`shape`:** select between region shapes: :cuboid, :cylinder/:disc, :sphere
##### Predefined/Optional Keywords:
**For cuboid region, related to a given center:**
- **`xrange`:** the range between [xmin, xmax] in units given by argument `range_unit` and relative to the given `center`; zero length for xmin=xmax=0. is converted to maximum possible length
- **`yrange`:** the range between [ymin, ymax] in units given by argument `range_unit` and relative to the given `center`; zero length for ymin=ymax=0. is converted to maximum possible length
- **`zrange`:** the range between [zmin, zmax] in units given by argument `range_unit` and relative to the given `center`; zero length for zmin=zmax=0. is converted to maximum possible length

**For cylindrical region, related to a given center:**
- **`radius`:** the radius between [0., radius] in units given by argument `range_unit` and relative to the given `center`
- **`height`:** the hight above and below a plane [-height, height] in units given by argument `range_unit` and relative to the given `center`
- **`direction`:** todo

**For spherical region, related to a given center:**
- **`radius`:** the radius between [0., radius] in units given by argument `range_unit` and relative to the given `center`

**Keywords related to all region shapes**
- **`range_unit`:** the units of the given ranges: :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`center`:** in units given by argument `range_unit`; by default [0., 0., 0.]; the box-center can be selected by e.g. [:bc], [:boxcenter], [value, :bc, :bc], etc..
- **`inverse`:** inverse the region selection = get the data outside of the region
- **`cell`:** take intersecting cells of the region boarder into account (true) or only the cells-centers within the selected region (false)
- **`verbose`:** print timestamp, selected vars and ranges on screen; default: set by the variable `verbose_mode`

"""

function subregion(dataobject::DataSetType, shape::Symbol=:cuboid;
    xrange::Array{<:Any,1}=[missing, missing],  # cuboid
    yrange::Array{<:Any,1}=[missing, missing],  # cuboid
    zrange::Array{<:Any,1}=[missing, missing],  # cuboid

    radius::Real=0.,              # cylinder, sphere
    height::Real=0.,              # cylinder
    direction::Symbol=:z,           # cylinder

    center::Array{<:Any,1}=[0.,0.,0.],   # all
    range_unit::Symbol=:standard,          # all
    cell::Bool=true,                        # hydro and gravity
    inverse::Bool=false,                    # all
    verbose::Bool=verbose_mode)             # all

    # subregion = wrapper over all subregion functions
    if shape == :cuboid
        if typeof(dataobject) == HydroDataType
            return subregioncuboid(dataobject,
                        xrange=xrange, yrange=yrange, zrange=zrange,
                        center=center,
                        range_unit=range_unit,
                        cell=cell,
                        inverse=inverse,
                        verbose=verbose)
        else
            return subregioncuboid(dataobject,
                        xrange=xrange, yrange=yrange, zrange=zrange,
                        center=center,
                        range_unit=range_unit,
                        inverse=inverse,
                        verbose=verbose)
        end

    elseif shape == :cylinder || shape == :disc
        if typeof(dataobject) == HydroDataType
            return subregioncylinder(dataobject,
                            radius=radius,
                            height=height,
                            center=center,
                            range_unit=range_unit,
                            direction=direction,
                            cell=cell,
                            inverse=inverse,
                            verbose=verbose)
        else
            return subregioncylinder(dataobject,
                            radius=radius,
                            height=height,
                            center=center,
                            range_unit=range_unit,
                            direction=direction,
                            inverse=inverse,
                            verbose=verbose)
        end

    elseif shape == :sphere
        if typeof(dataobject) == HydroDataType
            return subregionsphere(dataobject,
                            radius=radius,
                            center=center,
                            range_unit=range_unit,
                            cell=cell,
                            inverse=inverse,
                            verbose=verbose)
        else
            return subregionsphere(dataobject,
                            radius=radius,
                            center=center,
                            range_unit=range_unit,
                            inverse=inverse,
                            verbose=verbose)
        end
    end

end


text:
# Autogenerated wrapper script for HepMC3_jll for powerpc64le-linux-gnu-cxx11
export libHepMC3, libHepMC3search

JLLWrappers.@generate_wrapper_header("HepMC3")
JLLWrappers.@declare_library_product(libHepMC3, "libHepMC3.so.3")
JLLWrappers.@declare_library_product(libHepMC3search, "libHepMC3search.so.4")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libHepMC3,
        "lib/libHepMC3.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libHepMC3search,
        "lib/libHepMC3search.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -5.992552298143351 <= q <= 4.750344613921084)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.47930083070023954 + -0.2811056889538768*1/(1 + exp(-(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3])))) + -0.8918376604733185*1/(1 + exp(-(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3])))) + -0.8083246377220759*1/(1 + exp(-(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))))))) + 1/(1 + exp(-(-0.740550943643274 + -0.007926005977438866*1/(1 + exp(-(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3])))) + 0.328386397895855*1/(1 + exp(-(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3])))) + 0.6266373732109729*1/(1 + exp(-(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))))))) + 1/(1 + exp(-(-0.3322291763851517 + -0.3879865017411772*1/(1 + exp(-(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3])))) + 0.13951112377587105*1/(1 + exp(-(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3])))) + -0.16669374723112762*1/(1 + exp(-(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
export HMat;

@enum ADMISSTYPES ADMISS_STRONG ADMISS_STANDARD ADMISS_MILD ADMISS_WEAK;
@enum BLOCKTYPES LOWRANK DENSE HMAT;

type HMat{T<:Number}
    # variables
    height::    Int
    width::     Int
    level::     Int
    trg::       Array{Int,1}
    src::       Array{Int,1}
    blockType:: BLOCKTYPES
    childHMat:: Array{HMat{T},2}
    UMat::      Array{T,2}
    VMat::      Array{T,2}
    DMat::      Array{T,2}
    # global settings
    EPS::       Float64
    MAXRANK::   Int
    MINN::      Int

    HMat() = new()
end

function admiss(x,y,type_admiss)
    if ((type_admiss == ADMISS_STRONG) || (type_admiss == ADMISS_STANDARD))
        if maximum(x-y) > 1
            return true;
        else
            return false;
        end
    elseif type_admiss == ADMISS_MILD
        if sum(abs(x-y)) > 1
            return true;
        else
            return false;
        end
    elseif type_admiss == ADMISS_WEAK
        if x != y
            return true;
        else
            return false;
        end
    end
end

function Z2Cmapper(n,minn,tensor)
    if maximum(n) <= minn
        return reshape(tensor,prod(n));
    end
    nl = ifloor(n/2);
    nr = n-nl;
    nv = zeros(n);
    nrange = Array(Any,2,1);
    offset = 0;
    order = vec(zeros(Int64,prod(n),1));
    for i = 0:2^2-1
        for d = 0:2-1
            nv[d+1] = nl[d+1]*(1-mod(i>>d,2))+nr[d+1]*mod(i>>d,2);
            nrange[d+1] = nl[d+1]*mod(i>>d,2)+(1:ifloor(nv[d+1]));
        end
        order[offset+(1:prod(nv))] = Z2Cmapper(nv,minn,tensor[nrange[1],nrange[2]]);
        offset = offset + prod(nv);
    end
    return order;
end

function svdtrunc(A,eps,mR)
    (U, S, V) = svd(full(A),thin=true);
    if minimum(size(A))>1 && S[1] > 1e-15
        idx = find(find(S/S[1].>eps).<=mR);
        return U[:,idx], S[idx], V[:,idx];
    else
        return U, S, V;
    end
end

function HMatd2h{T}(D::AbstractMatrix{T}, nTrg, nSrc, type_admiss, idxTrg, idxSrc, level, EPS, MaxRank, minn)
    node = HMat{T}();
    node.height = prod(nTrg);
    node.width = prod(nSrc);
    node.trg = idxTrg;
    node.src = idxSrc;
    node.level = level;

    node.EPS = EPS;
    node.MAXRANK = MaxRank;
    node.MINN = minn;

    if admiss(idxTrg,idxSrc,type_admiss)
        (Utmp, Stmp, Vtmp) = svdtrunc(D,EPS,MaxRank);
        node.blockType = LOWRANK;
        if length(Stmp) > 0
            node.UMat = Utmp.*sqrt(Stmp)';
            node.VMat = Vtmp.*sqrt(Stmp)';
        else
            node.UMat = Utmp;
            node.VMat = Vtmp;
        end
    elseif maximum(nTrg) <= minn || maximum(nSrc) <= minn
        node.blockType = DENSE;
        node.DMat = full(D);
    else
        node.blockType = HMAT;
        node.childHMat = Array(HMat{T},4,4);
        trg = 2*idxTrg;
        src = 2*idxSrc;
        toffset = 0;
        for tx = 0:1, ty = 0:1
            trg = 2*idxTrg + [tx,ty];
            (txlen,tylen) = ifloor(nTrg/2).*(1-[tx,ty]) + iceil(nTrg/2).*[tx,ty];
            tRange = toffset + (1:txlen*tylen);
            toffset = toffset + txlen*tylen;

            soffset = 0;
            for sx = 0:1, sy = 0:1
                src = 2*idxSrc + [sx,sy];
                (sxlen,sylen) = ifloor(nSrc/2).*(1-[sx,sy]) + iceil(nSrc/2).*[sx,sy];
                sRange = soffset + (1:sxlen*sylen);
                soffset = soffset + sxlen*sylen;
                node.childHMat[tx*2+ty+1,sx*2+sy+1] = HMatd2h(D[tRange,sRange],[txlen,tylen],[sxlen,sylen],type_admiss,trg,src,level+1,EPS,MaxRank,minn);
            end
        end
    end
    return node;
end

include("../src/hcopy.jl");
include("../src/hadjoint.jl");
include("../src/hempty.jl");
include("../src/hidentity.jl");
include("../src/hmatvec.jl");
include("../src/hcompress.jl");
include("../src/huncompress.jl");
include("../src/hscale.jl");
include("../src/hadd.jl");
include("../src/hnorm.jl");
include("../src/hmul.jl");
include("../src/hinv.jl");
include("../src/hNewton.jl");

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
module MusicTransformer

using CUDA
using Flux
using Transformers

include("layers.jl")
include("models.jl")
include("musicencoders.jl")
include("generate.jl")
include("pretrain.jl")
include("datasets.jl")

device = Flux.cpu

if CUDA.has_cuda()
    device = Flux.gpu
    println("GPU is enabled.")
    CUDA.allowscalar(false)
end

function call_registers()
    register_configs(pretrained_configs)
    register_datasets(dataset_configs)
end

@init call_registers()

end


text:
using Test, LinearAlgebra
using MethodOfLines

@testset "finite-difference weights from fornberg(1988) & fornberg(2020)" begin
    order = 2; z = 0.0; x = [-1, 0, 1.0];
    @test MethodOfLines.calculate_weights(order, z, x) == [1,-2,1]  # central difference of second-derivative with unit-step

    order = 1; z = 0.0; x = [-1., 1.0];
    @test MethodOfLines.calculate_weights(order, z, x) == [-0.5,0.5] # central difference of first-derivative with unit step

    order = 1; z = 0.0; x = [0, 1];
    @test MethodOfLines.calculate_weights(order, z, x) == [-1, 1] # forward difference

    order = 1; z = 1.0; x = [0, 1];
    @test MethodOfLines.calculate_weights(order, z, x) == [-1, 1] # backward difference

    # forward-diff of third derivative with order of accuracy == 3
    order = 3; z = 0.0; x = [0,1,2,3,4,5]
    @test MethodOfLines.calculate_weights(order, z, x) == [-17/4, 71/4, -59/2, 49/2, -41/4, 7/4]

    order = 3; z = 0.0; x = collect(-3:3)
    d, e = MethodOfLines.calculate_weights(order, z, x;dfdx = true)
    @test d ≈ [-167/18000, -963/2000, -171/16,0,171/16,963/2000,167/18000]
    @test e ≈ [-1/600,-27/200,-27/8,-49/3,-27/8,-27/200,-1/600]

    order = 3; z = 0.0; x = collect(-4:4)
    d, e = MethodOfLines.calculate_weights(order, z, x;dfdx = true)
    @test d ≈ [-2493/5488000, -12944/385875, -87/125 ,-1392/125,0,1392/125,87/125,12944/385875,2493/5488000]
    @test e ≈ [-3/39200,-32/3675,-6/25,-96/25,-205/12, -96/25, -6/25,-32/3675,-3/39200]
end


text:
#==============================================================================#
# sign.jl
#
# AWS Request Signing.
#
# Copyright Sam O'Connor 2014 - All rights reserved
#==============================================================================#

import Nettle: digest, hexdigest

using URIParser

function sign!(r::AWSRequest, t = now(Dates.UTC))

    if r[:service] == "sdb"
        sign_aws2!(r, t)
    elseif r[:service] == "ses"
        sign_aws3!(r, t)
    else
        sign_aws4!(r, t)
    end
end

# Create AWS Signature Version 2 Authentication query parameters.
# http://docs.aws.amazon.com/general/latest/gr/signature-version-2.html

function sign_aws2!(r::AWSRequest, t)

    uri = URI(r[:url])

    query = Dict{AbstractString,AbstractString}()
    for elem in split(r[:content], '&', keep=false)
        (n, v) = split(elem, "=")
        query[n] = unescape(v)
    end

    r[:headers]["Content-Type"] =
        "application/x-www-form-urlencoded; charset=utf-8"

    query["AWSAccessKeyId"] = r[:creds].access_key_id
    query["Expires"] = Dates.format(t + Dates.Second(120),
                                    "yyyy-mm-ddTHH:MM:SSZ")
    query["SignatureVersion"] = "2"
    query["SignatureMethod"] = "HmacSHA256"
    if r[:creds].token != ""
        query["SecurityToken"] = r[:creds].token
    end

    query = [(k, query[k]) for k in sort(collect(keys(query)))]

#FIXME  see https://github.com/JuliaWeb/URIParser.jl/pull/31
#    to_sign = "POST\n$(uri.host)\n$(uri.path)\n$(format_query_str(query))"
    to_sign = "POST\n$(uri.host)\n$(uri.path)\n$(escape_with(format_query_str(query), "()"))"

    secret = r[:creds].secret_key
    push!(query, ("Signature", digest("sha256", secret, to_sign)
                               |> base64encode |> strip))

    r[:content] = format_query_str(query)
end

# Create AWS3 Authentication Headers.
# http://docs.aws.amazon.com/ses/latest/DeveloperGuide/query-interface-authentication.html

function sign_aws3!(r::AWSRequest, t)

    date = string(Dates.format(t, "e, d u yyyy H:M:S"), " GMT")

    r[:headers]["Date"] = date
    r[:headers]["X-Amzn-Authorization"] = string(
        "AWS3-HTTPS AWSAccessKeyId=", r[:creds].access_key_id, ",",
        "Algorithm=HmacSHA256,",
        "Signature=", digest("sha256", r[:creds].secret_key, date)
                      |> base64encode |> strip)
end

# Create AWS Signature Version 4 Authentication Headers.
# http://docs.aws.amazon.com/general/latest/gr/signature-version-4.html

function sign_aws4!(r::AWSRequest, t)

    # ISO8601 date/time strings for time of request...
    date = Dates.format(t,"yyyymmdd")
    datetime = Dates.format(t,"yyyymmddTHHMMSSZ")

    # Authentication scope...
    scope = [date, r[:region], r[:service], "aws4_request"]

    # Signing key generated from today's scope string...
    signing_key = string("AWS4", r[:creds].secret_key)
    for element in scope
        signing_key = digest("sha256", signing_key, element)
    end

    # Authentication scope string...
    scope = join(scope, "/")

    # SHA256 hash of content...
    content_hash = hexdigest("sha256", r[:content])

    # HTTP headers...
    delete!(r[:headers], "Authorization")
    merge!(r[:headers], Dict(
        "x-amz-content-sha256" => content_hash,
        "x-amz-date"           => datetime,
        "Content-MD5"          => base64encode(digest("md5", r[:content]))
    ))
    if r[:creds].token != ""
        r[:headers]["x-amz-security-token"] = r[:creds].token
    end

    # Sort and lowercase() Headers to produce canonical form...
    canonical_headers = ["$(lowercase(k)):$(strip(v))" for (k,v) in r[:headers]]
    signed_headers = join(sort([lowercase(k) for k in keys(r[:headers])]), ";")

    # Sort Query String...
    uri = URI(r[:url])
    query = query_params(uri)
    query = [(k, query[k]) for k in sort(collect(keys(query)))]

    # Create hash of canonical request...
    canonical_form = string(r[:verb], "\n",
                            escape_with(uri.path, path_esc_chars), "\n",
                            format_query_str(query), "\n",
                            join(sort(canonical_headers), "\n"), "\n\n",
                            signed_headers, "\n",
                            content_hash)
    canonical_hash = hexdigest("sha256", canonical_form)

    # Create and sign "String to Sign"...
    string_to_sign = "AWS4-HMAC-SHA256\n$datetime\n$scope\n$canonical_hash"
    signature = hexdigest("sha256", signing_key, string_to_sign)

    # Append Authorization header...
    r[:headers]["Authorization"] = string(
        "AWS4-HMAC-SHA256 ",
        "Credential=$(r[:creds].access_key_id)/$scope, ",
        "SignedHeaders=$signed_headers, ",
        "Signature=$signature"
    )
end

const path_esc_chars = filter(c->c!='/', URIParser.unescaped)

#==============================================================================#
# End of file.
#==============================================================================#


text:
function run_simulation(set_sizes; fun=conjunctive_set, kwargs...)
    results = DataFrame[]
    for n in set_sizes
        experiment = Experiment(set_size=n, populate_visicon=fun,
            n_trials=10^4)
        run_condition!(experiment; kwargs...)
        df = DataFrame(experiment.data)
        df_present = filter(x->x.target_present == :present, df)
        hit_rate = mean(df_present.response .== :present)
        temp = by(df, [:target_present,:response], :rt=>mean)
        temp[!,:distractors] .= n
        temp[!,:hit_rate] .= hit_rate
        push!(results, temp)
    end
    return vcat(results...)
end


text:
# Autogenerated wrapper script for P4est_jll for armv7l-linux-gnueabihf
export libp4est, libsc

using Zlib_jll
JLLWrappers.@generate_wrapper_header("P4est")
JLLWrappers.@declare_library_product(libp4est, "libp4est.so.0")
JLLWrappers.@declare_library_product(libsc, "libsc.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        libp4est,
        "lib/libp4est.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsc,
        "lib/libsc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    FourierV <: AbstractVApproximator{Int}

    struct FourierV <: AbstractVApproximator{Int}
        weights::Vector{Float64}
    end

Using Fourier cosine basis to approximate the state value.
`weights` is the featur vector.

See more details at Section (9.5.2) on Page 211 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
struct FourierV <: AbstractVApproximator{Int}
    weights::Vector{Float64}
end

"""
    FourierV(order::Int)

By specifying the `order`, feature vector will be initialized with 0.
"""
FourierV(order::Int) = FourierV(zeros(Float64, order+1))

function (f::FourierV)(s)
    sum(w * cos((i-1) * π * s) for (i, w) in enumerate(f.weights))
end

function update!(f::FourierV, s, e)
    for i in 1:length(f.weights)
        f.weights[i] += e * cos((i-1) * π * s)
    end
end


text:
import StatsBase: sample

# Model selection

# Taken from https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2017.0009
"""
	$(SIGNATURES)

Computes the Akaike Information Criterion (AIC) given the free parameters `k` for the data `X` and its
estimate `Y` of the model. `likelihood` can be any function of `X` and `Y`.
"""
function AIC(k::Int64, X::AbstractArray, Y::AbstractArray; likelihood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return 2*k - 2*log(likelihood(X, Y))
end

# Taken from https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2017.0009
"""
	$(SIGNATURES)

Computes the Akaike Information Criterion compensated for finite samples (AICC) given the free parameters `k` for the data `X` and its
estimate `Y` of the model. `likelihood` can be any function of `X` and `Y`.
"""
function AICC(k::Int64, X::AbstractMatrix, Y::AbstractMatrix; likelihood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return AIC(k, X, Y, likelihood = likelihood)+ 2*(k+1)*(k+2)/(size(X)[2]-k-2)
end

function AICC(k::Int64, X::AbstractVector, Y::AbstractVector; likelihood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return AIC(k, X, Y, likelihood = likelihood)+ 2*(k+1)*(k+2)/(length(X)-k-2)
end

# Double check on that
# Taken from https://www.immagic.com/eLibrary/ARCHIVES/GENERAL/WIKIPEDI/W120607B.pdf
"""
	$(SIGNATURES)

Computes Bayes Information Criterion (BIC) given the free parameters `k` for the data `X` and its
estimate `Y` of the model. `likelihood` can be any function of `X` and `Y`.
"""
function BIC(k::Int64, X::AbstractMatrix, Y::AbstractMatrix; likelihood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return - 2*log(likelihood(X, Y)) + k*log(size(X)[2])
end

function BIC(k::Int64, X::AbstractVector, Y::AbstractVector; likelihood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return - 2*log(likelihood(X, Y)) + k*log(length(X))
end

# Optimal Shrinkage for data in presence of white noise
# See D. L. Donoho and M. Gavish, "The Optimal Hard Threshold for Singular
# Values is 4/sqrt(3)", http://arxiv.org/abs/1305.5870
# Code taken from https://github.com/erichson/optht

function optimal_svht(m::Int64, n::Int64; known_noise::Bool = false)
    @assert m/n > 0
    @assert m/n <= 1

    β = m/n
    ω = (8*β) / (β+1+sqrt(β^2+14β+1))
    c = sqrt(2*(β+1)+ω)

    if known_noise
        return c
    else
        median = median_marcenko_pastur(β)
        return c / sqrt(median)
    end
end

function marcenko_pastur_density(t, lower, upper, beta)
    sqrt((upper-t).*(t-lower))./(2π*beta*t)
end

function incremental_marcenko_pastur(x, beta, gamma)
    @assert beta <= 1
    upper = (1+sqrt(beta))^2
    lower = (1-sqrt(beta))^2

    @inline marcenko_pastur(x) = begin
        if (upper-x)*(x-lower) > 0
            return marcenko_pastur_density(x, lower, upper, beta)
        else
            return zero(eltype(x))
        end
    end

    if gamma ≈ zero(eltype(gamma))
        i, ϵ = quadgk(x->(x^gamma)*marcenko_pastur(x), x, upper)
        return i
    else
        i, ϵ = quadgk(x->marcenko_pastur(x), x, upper)
        return i
    end
end

function median_marcenko_pastur(beta)
    @assert 0 < beta <= 1
    upper = (1+sqrt(beta))^2
    lower = (1-sqrt(beta))^2
    change = true
    x = ones(eltype(upper), 5)
    y = similar(x)
    while change && (upper - lower > 1e-5)
        x = range(lower, upper, length = 5)
        for (i,xi) in enumerate(x)
            y[i] = one(eltype(x)) - incremental_marcenko_pastur(xi, beta, 0)
        end
        any(y .< 0.5) ? lower = maximum(x[y .< 0.5]) : change = false
        any(y .> 0.5) ? upper = minimum(x[y .> 0.5]) : change = false
    end
    return (lower+upper)/2
end

"""
    $(SIGNATURES)

Compute a feature reduced version of the data array `X` via thresholding the
singular values by computing the [optimal threshold for singular values](http://arxiv.org/abs/1305.5870).
"""
function optimal_shrinkage(X::AbstractArray{T, 2}) where T <: Number
    m,n = minimum(size(X)), maximum(size(X))
    U, S, V = svd(X)
    τ = optimal_svht(m,n)
    inds = S .>= τ*median(S)
    return U[:, inds]*Diagonal(S[inds])*V[:, inds]'
end

"""
    $(SIGNATURES)

Compute a feature reduced version of the data array `X` inplace via thresholding the
singular values by computing the [optimal threshold for singular values](http://arxiv.org/abs/1305.5870).
"""
function optimal_shrinkage!(X::AbstractArray{T, 2}) where T <: Number
    m,n = minimum(size(X)), maximum(size(X))
    U, S, V = svd(X)
    τ = optimal_svht(m,n)
    inds = S .>= τ*median(S)
    X .= U[:, inds]*Diagonal(S[inds])*V[:, inds]'
    return
end

"""
	($SIGNATURES)

Randomly selects `n` bursts of data with size `samplesize` from the data `X`.

Randomly selects `n` bursts of data with size `samplesize` from the data `X` and `Y`.

Randomly selects `n` bursts of data within a time window `period` from the data `X`. The time information
has to be provided in `t`.
"""
@inline function burst_sampling(x::AbstractArray, samplesize::Int64, bursts::Int64)
    @assert size(x)[end] >= samplesize*bursts "Bursting impossible. Please provide more data or reduce bursts or samplesize."
    inds = sample(1:size(x)[end]-samplesize, bursts, replace = false)
    inds = sort(unique(vcat([collect(i:i+samplesize) for i in inds]...)))
    return resample(x, inds)
end

@inline function burst_sampling(x::AbstractArray, y::AbstractArray, samplesize::Int64, bursts::Int64)
    @assert size(x)[end] >= samplesize*bursts "Bursting impossible. Please provide more data or reduce bursts or samplesize"
    @assert size(x)[end] == size(y)[end]
    inds = sample(1:size(x)[end]-samplesize, bursts, replace = false)
    inds = sort(unique(vcat([collect(i:i+samplesize) for i in inds]...)))
    return resample(x, inds), resample(y, inds)
end

@inline function burst_sampling(x::AbstractArray, t::AbstractVector, period::T, bursts::Int64) where T <: AbstractFloat
    @assert period > zero(typeof(period)) "Sampling period has to be positive."
    @assert size(x)[end] == size(t)[end] "Provide consistent data."
    @assert bursts >= 1 "Number of bursts has to be positive."
    @assert t[end]-t[1]>= period*bursts "Bursting impossible. Please provide more data or reduce bursts or samplesize"
    t_ids = zero(eltype(t)) .<= t .- period  .<= t[end] .- 2*period
    samplesize = Int64(floor(period/(t[end]-t[1])*length(t)))
    inds = sample(collect(1:length(t))[t_ids], bursts, replace = false)
    inds = sort(unique(vcat([collect(i:i+samplesize) for i in inds]...)))
    return resample(x, inds), resample(t, inds)
end

"""
	$(SIGNATURES)

Returns the subsampled `X` with only every `n`-th entry.

Returns the subsampled `X` with a a minimum period of `dt` between two data points. `t` provides the
time information.
"""
@inline function subsample(x::AbstractVector, frequency::Int64)
    @assert frequency > 0 "Sampling frequency has to be positive."
    return x[1:frequency:end]
end

@inline function subsample(x::AbstractArray, frequency::Int64)
    @assert frequency > 0 "Sampling frequency has to be positive."
    return x[:, 1:frequency:end]
end

@inline function subsample(x::AbstractArray, t::AbstractVector, period::T) where T <: AbstractFloat
    @assert period > zero(typeof(period)) "Sampling period has to be positive."
    @assert size(x)[end] == size(t)[end] "Provide consistent data."
    @assert t[end]-t[1]>= period "Subsampling impossible. Sampling period exceeds time window."
    idx = Int64[1]
    t_now = t[1]
    @inbounds for (i, t_current) in enumerate(t)
        if t_current - t_now >= period
            push!(idx, i)
            t_now = t_current
        end
    end
    return resample(x, idx), resample(t, idx)
end

@inline function resample(x::AbstractArray{T,1}, indx::AbstractArray{Int64}) where T <: Number
    @assert maximum(indx) <= length(x) "Sampling index has to be consistent with array dimensions."
    @assert minimum(indx) >= 1 "Sampling index has to be consistent with array dimensions."
    return x[indx]
end

@inline function resample(x::AbstractArray{T,2}, indx::AbstractArray{Int64}) where T <: Number
    @assert maximum(indx) <= size(x, 2) "Sampling index has to be consistent with array dimensions."
    @assert minimum(indx) >= 1 "Sampling index has to be consistent with array dimensions."
    return x[:, indx]
end


text:
using DownstreamFail
using Test

@testset "Downstream Fail" begin
    @test true == false
end


text:
abstract type AbstractNode{T<:AbstractFloat,O<:Intersectable{T}} <: Intersectable{T} end
abstract type AbstractTree{T<:AbstractFloat,O} end

mutable struct BVH{T,O,Tree<:AbstractTree{T,O},Node<:AbstractNode{T,O}}
    tree::Tree
    root::Node
end

aabb(n::Node) where {Node<:AbstractNode} = n.bbox
Base.intersect(ray::Ray{T}, n::Node) where {T,Node<:AbstractNode{T}} = intersect(ray, n.bbox)

function traverse_reduce(tree::Tree, node::Node,
                         binary::N, binop::Function = +,
                         leafop::Function = n -> one(N)) where {Tree<:AbstractTree, Node<:AbstractNode, N}
    if is_inner(node)
        binary + binop(traverse_reduce(tree, tree[node.left], binary, binop, leafop),
                       traverse_reduce(tree, tree[node.right], binary, binop, leafop))
    else
        leafop(node)
    end
end

depth(bvh::B) where {B<:BVH} = traverse_reduce(bvh.tree, bvh.root, 1, max)
count_leaves(bvh::B) where {B<:BVH} = traverse_reduce(bvh.tree, bvh.root, 0)
count_objs(bvh::B) where {B<:BVH} = traverse_reduce(bvh.tree, bvh.root, 0, +, n -> length(objects(n)))

Base.show(io::IO, bvh::B) where {B<:BVH} =
    write(io, "$B of depth $(depth(bvh)) with $(count_leaves(bvh)) leaf nodes with $(count_objs(bvh)) objects")

export BVH


text:
using DependencyParser.DependencyParsing

function build_updated_gold_state(config::Configuration; system::ParsingSystem = ArcEagerSystem())
  tree = build_gold_tree()

  GoldState(tree, config, system)
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
using BinomialGPU
using CUDA

using BenchmarkTools
using Test

@testset "BinomialGPU.jl" begin
    @testset "constant parameters" begin
        n = 128
        p = 0.5

        @testset "A of dim $(length(Adims))" for Adims in [[2,], [2, 4], [2, 4, 8], [2, 4, 8, 16]]
            A = CUDA.zeros(Int, Tuple(Adims))
            @test rand_binomial!(A, count = n, prob = p) isa CuArray{Int}
            @test minimum(rand_binomial!(A, count = n, prob = p)) >= 0
            @test maximum(rand_binomial!(A, count = n, prob = p)) <= n
        end
    end

    @testset "parameter arrays" begin
        @testset "A of dim $(length(Adims))" for Adims in [[2,], [2, 4], [2, 4, 8], [2, 4, 8, 16]]
            A = CUDA.zeros(Int, Tuple(Adims))
            @testset "count of dim 0, prob of dim $j" for j in eachindex(Adims)
                pdim = Adims[1:j]
                n = 128
                ps = CUDA.rand(pdim...)
                @test rand_binomial!(A, count = n, prob = ps) isa CuArray{Int}
                @test minimum(rand_binomial!(A, count = n, prob = ps)) >= 0
                @test minimum(n .- rand_binomial!(A, count = n, prob = ps)) >= 0
            end

            @testset "count of dim $i, prob of dim 0" for i in eachindex(Adims)
                ndim = Adims[1:i]
                ns = CUDA.fill(128, Tuple(ndim))
                p = 0.5
                @test rand_binomial!(A, count = ns, prob = p) isa CuArray{Int}
                @test minimum(rand_binomial!(A, count = ns, prob = p)) >= 0
                @test minimum(ns .- rand_binomial!(A, count = ns, prob = p)) >= 0
            end

            @testset "count of dim $i, prob of dim $j" for i in eachindex(Adims), j in eachindex(Adims)
                ndim = Adims[1:i]
                pdim = Adims[1:j]
                ns = CUDA.fill(128, Tuple(ndim))
                ps = CUDA.rand(pdim...)
                n = 128
                p = 0.5
                @test rand_binomial!(A, count = ns, prob = ps) isa CuArray{Int}
                @test minimum(rand_binomial!(A, count = ns, prob = ps)) >= 0
                @test minimum(ns .- rand_binomial!(A, count = ns, prob = ps)) >= 0

                # wrong size in the last dimension
                for k in 1:i, l in 1:j
                    ndim[k] += 1
                    ns = CUDA.fill(128, Tuple(ndim))
                    ps = CUDA.rand(pdim...)
                    @test_throws DimensionMismatch rand_binomial!(A, count = ns, prob = ps)
                    ndim[k] -= 1
                    pdim[l] += 1
                    ns = CUDA.fill(128, Tuple(ndim))
                    ps = CUDA.rand(pdim...)
                    @test_throws DimensionMismatch rand_binomial!(A, count = ns, prob = ps)
                    pdim[l] -= 1
                end

                # wrong number of dimensions
                if i == length(Adims)
                    push!(ndim, 32)
                    ns = CUDA.fill(128, Tuple(ndim))
                    ps = CUDA.rand(pdim...)
                    @test_throws DimensionMismatch rand_binomial!(A, count = ns, prob = ps)
                    pop!(ndim)
                end
                if j == length(Adims)
                    push!(pdim, 32)
                    ns = CUDA.fill(128, Tuple(ndim))
                    ps = CUDA.rand(pdim...)
                    @test_throws DimensionMismatch rand_binomial!(A, count = ns, prob = ps)
                    pop!(pdim)
                end
            end
        end
    end

    @testset "bad parameter values" begin
        # bad parameter values default
        A = CUDA.zeros(Int, 256)
        @test rand_binomial!(A, count = -1, prob = 0.5) == CUDA.zeros(256) # negative counts are equivalent to zero
        @test rand_binomial!(A, count = 2, prob = -0.1) == CUDA.zeros(256) # negative probabilities are equivalent to zero
        @test rand_binomial!(A, count = 2, prob = 1.5) == CUDA.fill(2, 256) # probabilities greater than 1 are equivalent to 1
        @test_throws MethodError rand_binomial!(A, count = 5., prob = 0.5) # non-integer counts throw an error
    end

    @testset "benchmarks" begin
        # benchmarks
        A = CUDA.zeros(Int, 1024, 1024)
        n = 128
        p = 0.5
        ns = CUDA.fill(128, (1024, 1024))
        ps = CUDA.rand(1024, 1024)
        println("")
        println("Benchmarking constant parameter array: should run in less than 2ms on an RTX20xx card")
        display(@benchmark CUDA.@sync rand_binomial!($A, count = $n, prob = $p))
        println("")
        println("Benchmarking full parameter array: should run in less than 2ms on an RTX20xx card")
        display(@benchmark CUDA.@sync rand_binomial!($A, count = $ns, prob = $ps))
        println("")
    end
end


text:
function get_tendencies!(ξ_tend, D_tend, Tₐ_tend, pₛ_tend, tr_tend, j2)

    # =========================================================================
    # Computation of grid-point tendencies (converted to spectral at the end of
    # grtend)
    # =========================================================================

    get_grid_point_tendencies!(ξ_tend, D_tend, Tₐ_tend, pₛ_tend, tr_tend, 1, j2)

    # =========================================================================
    # Computation of spectral tendencies
    # =========================================================================

    if α < 0.5
        get_spectral_tendencies!(D_tend, Tₐ_tend, pₛ_tend, j2)
    else
        get_spectral_tendencies!(D_tend, Tₐ_tend, pₛ_tend, 1)

        # Implicit correction
        implicit_terms!(D_tend, Tₐ_tend, pₛ_tend)
    end
end

# Compute non-linear tendencies in grid point space from dynamics and
# physical parametrizations, and convert them to spectral tendencies
# dF/dt = T_dyn(F(J2)) + T_phy(F(J1))
#   Input:  j1 = time level index for physical tendencies
#           j2 = time level index for dynamical tendencies
#   Output: ξ_tend = spectral tendency of vorticity
#           D_tend = spectral tendency of divergence
#           Tₐ_tend   = spectral tendency of temperature
#           pₛ_tend  = spectral tendency of log(p_s)
#           tr_tend  = spectral tendency of tracers
function get_grid_point_tendencies!(ξ_tend, D_tend, Tₐ_tend, pₛ_tend, tr_tend, j1, j2)
    # =========================================================================
    # Convert prognostics to grid point space
    # =========================================================================

    u_grid = zeros(RealType, nlon, nlat, nlev)
    v_grid = zeros(RealType, nlon, nlat, nlev)
    Tₐ_grid = zeros(RealType, nlon, nlat, nlev)
    ξ_grid = zeros(RealType, nlon, nlat, nlev)
    D_grid = zeros(RealType, nlon, nlat, nlev)
    tr_grid = zeros(RealType, nlon, nlat, nlev, n_trace)
    dumc = zeros(Complex{RealType}, mx, nx, 2)

    for k in 1:nlev
        ξ_grid[:,:,k] = spec_to_grid(ξ[:,:,k,j2])
        D_grid[:,:,k] = spec_to_grid(D[:,:,k,j2])
        Tₐ_grid[:,:,k]   = spec_to_grid(Tₐ[:,:,k,j2])

        for itr in 1:n_trace
            tr_grid[:,:,k,itr] = spec_to_grid(tr[:,:,k,j2,itr])
        end

        uvspec!(ξ[:,:,k,j2], D[:,:,k,j2], @view(dumc[:,:,1]), @view(dumc[:,:,2]))
        v_grid[:,:,k] = spec_to_grid(dumc[:,:,2], scale=true)
        u_grid[:,:,k] = spec_to_grid(dumc[:,:,1], scale=true)

        for j in 1:nlat
            for i in 1:nlon
                ξ_grid[i,j,k] += f[j]
            end
        end
    end

    umean = zeros(RealType, nlon, nlat)
    vmean = zeros(RealType, nlon, nlat)
    dmean = zeros(RealType, nlon, nlat)

    for k in 1:nlev
        umean += u_grid[:,:,k]*dhs[k]
        vmean += v_grid[:,:,k]*dhs[k]
        dmean += D_grid[:,:,k]*dhs[k]
    end

    # Compute tendency of log(surface pressure)
    # pₛ(1,1,j2)=zero
    grad!(pₛ[:,:,j2], @view(dumc[:,:,1]), @view(dumc[:,:,2]))
    px = spec_to_grid(dumc[:,:,1], scale=true)
    py = spec_to_grid(dumc[:,:,2], scale=true)

    pₛ_tend = grid_to_spec(-umean.*px - vmean.*py)
    pₛ_tend[1,1] = Complex{RealType}(zero)

    # Compute "vertical" velocity
    σ_tend = zeros(RealType, nlon, nlat, nlev+1)
    σ_m = zeros(RealType, nlon, nlat, nlev+1)

    # (The following combination of terms is utilized later in the
    #     temperature equation)
    puv = zeros(RealType, nlon, nlat, nlev)
    for k in 1:nlev
        puv[:,:,k] = (u_grid[:,:,k] - umean).*px + (v_grid[:,:,k] - vmean).*py
    end

    for k in 1:nlev
        σ_tend[:,:,k+1] = σ_tend[:,:,k] - dhs[k]*(puv[:,:,k] + D_grid[:,:,k] - dmean)
        σ_m[:,:,k+1] = σ_m[:,:,k] - dhs[k]*puv[:,:,k]
    end

    # Subtract part of temperature field that is used as reference for implicit terms
    t_grid_anom = zeros(RealType, nlon, nlat, nlev)

    for k in 1:nlev
        t_grid_anom[:,:,k] = Tₐ_grid[:,:,k] .- tref[k]
    end

    # Zonal wind tendency
    temp = zeros(RealType, nlon, nlat, nlev+1)
    u_tend = zeros(RealType, nlon, nlat, nlev)

    for k in 2:nlev
        temp[:,:,k] = σ_tend[:,:,k].*(u_grid[:,:,k] - u_grid[:,:,k-1])
    end

    for k in 1:nlev
        u_tend[:,:,k] = v_grid[:,:,k].*ξ_grid[:,:,k] - t_grid_anom[:,:,k]*R.*px
            - (temp[:,:,k+1] + temp[:,:,k])*dhsr[k]
    end

    # Meridional wind tendency
    v_tend = zeros(RealType, nlon, nlat, nlev)

    for k in 2:nlev
        temp[:,:,k] = σ_tend[:,:,k].*(v_grid[:,:,k] - v_grid[:,:,k-1])
    end

    for k in 1:nlev
        v_tend[:,:,k] = -u_grid[:,:,k].*ξ_grid[:,:,k] - t_grid_anom[:,:,k]*R.*py
            - (temp[:,:,k+1] + temp[:,:,k])*dhsr[k]
    end

    # Temperature tendency
    Tₐ_tend_grid = zeros(RealType, nlon, nlat, nlev)

    for k in 2:nlev
        temp[:,:,k] = σ_tend[:,:,k].*(t_grid_anom[:,:,k] - t_grid_anom[:,:,k-1])
            + σ_m[:,:,k].*(tref[k] - tref[k-1])
    end

    for k in 1:nlev
        Tₐ_tend_grid[:,:,k] = t_grid_anom[:,:,k].*D_grid[:,:,k]
            - (temp[:,:,k+1] + temp[:,:,k])*dhsr[k]
            + fsgr[k]*t_grid_anom[:,:,k].*(σ_tend[:,:,k+1] + σ_tend[:,:,k])
            + tref3[k]*(σ_m[:,:,k+1] + σ_m[:,:,k])
            + akap*(Tₐ_grid[:,:,k].*puv[:,:,k] - t_grid_anom[:,:,k].*dmean)
    end

    # Tracer tendency
    tr_tend_grid = zeros(RealType, nlon, nlat, nlev, n_trace)
    for itr in 1:n_trace
        for k in 2:nlev
            temp[:,:,k] = σ_tend[:,:,k].*(tr_grid[:,:,k,itr] - tr_grid[:,:,k-1,itr])
        end

        temp[:,:,2:3] .= zero

        for k in 1:nlev
            tr_tend_grid[:,:,k,itr] = tr_grid[:,:,k,itr].*D_grid[:,:,k]
                - (temp[:,:,k+1] + temp[:,:,k])*dhsr[k]
        end
    end

    # =========================================================================
    # Convert tendencies to spectral space
    # =========================================================================

    for k in 1:nlev
        #  Convert u and v tendencies to vor and div spectral tendencies
        #  vdspec takes a grid u and a grid v and converts them to
        #  spectral vor and div
        vdspec!(u_tend[:,:,k], v_tend[:,:,k],
            @view(ξ_tend[:,:,k]), @view(D_tend[:,:,k]), true)

        # Divergence tendency
        # add -lapl(0.5*(u**2+v**2)) to div tendency
        D_tend[:,:,k] = D_tend[:,:,k]
            - ∇²(grid_to_spec(half*(u_grid[:,:,k].^two + v_grid[:,:,k].^two)))

        # Temperature tendency
        # and add div(vT) to spectral t tendency
        vdspec!(-u_grid[:,:,k].*t_grid_anom[:,:,k], -v_grid[:,:,k].*t_grid_anom[:,:,k],
            dumc[:,:,1], @view(Tₐ_tend[:,:,k]), true)
        Tₐ_tend[:,:,k] += grid_to_spec(Tₐ_tend_grid[:,:,k])

        # tracer tendency
        for itr in 1:n_trace
            vdspec!(-u_grid[:,:,k].*tr_grid[:,:,k,itr], -v_grid[:,:,k].*tr_grid[:,:,k,itr],
                dumc[:,:,1], @view(tr_tend[:,:,k,itr]), true)
            tr_tend[:,:,k,itr] += grid_to_spec(tr_tend_grid[:,:,k,itr])
        end
    end
end

# Compute spectral tendencies of divergence, temperature  and log(surface pressure)
# Input/output : D_tend = divergence tendency (spectral)
#                Tₐ_tend   = temperature tendency (spectral)
#                pₛ_tend  = tendency of log_surf.pressure (spectral)
#                j2    = time level index (1 or 2)
function get_spectral_tendencies!(D_tend, Tₐ_tend, pₛ_tend, j2)
    # Vertical mean div and pressure tendency
    dmeanc = zeros(Complex{RealType}, mx, nx)
    for k in 1:nlev
        dmeanc = dmeanc + D[:,:,k,j2]*dhs[k]
    end

    pₛ_tend = pₛ_tend - dmeanc
    pₛ_tend[1,1] = Complex{RealType}(zero)

    # Sigma-dot "velocity" and temperature tendency
    σ_tendc = zeros(Complex{RealType}, mx, nx, nlev+1)

    for k in 1:nlev - 1
        σ_tendc[:,:,k+1] = σ_tendc[:,:,k] - dhs[k]*(D[:,:,k,j2] - dmeanc)
    end

    dumk = zeros(Complex{RealType}, mx, nx, nlev+1)

    for k in 2:nlev
        dumk[:,:,k] = σ_tendc[:,:,k]*(tref[k] - tref[k-1])
    end

    for k in 1:nlev
        Tₐ_tend[:,:,k] -= (dumk[:,:,k+1] + dumk[:,:,k])*dhsr[k]
            + tref3[k]*(σ_tendc[:,:,k+1] + σ_tendc[:,:,k]) - tref2[k]*dmeanc
    end

    # Geopotential and divergence tendency
    get_geopotential(Tₐ[:,:,:,j2])
    for k in 1:nlev
        D_tend[:,:,k] -= ∇²(ϕ[:,:,k] + R*tref[k]*pₛ[:,:,j2])
    end
end


text:
# routines for front end.

# fits a positive function via SDP2.
function fitnDdensity(  y::Vector{T},
                        X::Vector{Vector{T}},
                        μ::T, # regularization weight.
                        θ::KT,
                        zero_tol::T,
                        max_iters::Int ) where {T,KT}
    #
    # fit.
    K = constructkernelmatrix(X, θ)
    α_solution = fitpdfviaSDP(y, K, μ, max_iters)
    α = clamp.(α_solution, zero_tol, Inf)
    η = RKHSProblemType(α, X, θ, μ)

    return η, α
end

# function fitnDdensityRiemannian(  y::Vector{T},
#                                 X::Vector{Vector{T}},
#                                 μ::T, # regularization weight.
#                                 θ::KT,
#                                 zero_tol::T,
#                                 max_iter::Int;
#                                 α_initial = ones(T, length(y)),     ## initial guess.
#                                 verbose_flag = false,
#                                 max_iter_tCG = 100,
#                                 ρ_lower_acceptance = 0.2, # recommended to be less than 0.25
#                                 ρ_upper_acceptance = 5.0,
#                                 minimum_TR_radius::T = 1e-3,
#                                 maximum_TR_radius::T = 10.0,
#                                 norm_df_tol = 1e-5,
#                                 objective_tol = 1e-5,
#                                 avg_Δf_tol = 0.0, #1e-12 #1e-5
#                                 avg_Δf_window = 10,
#                                 max_idle_update_count = 50,
#                                 g::Function = pp->1.0/(dot(pp,pp)+1.0),
#                                 𝑟 = 1e-2 ) where {T,KT}
#     #
#     # fit.
#     K = constructkernelmatrix(X, θ)
#     α_solution, f_α_array_unused, norm_df_array_unused,
#         num_iters_unused = RiemannianOptim.solveRKHSℝpproblem(y, K, μ;
#                             max_iter = max_iter,
#                             α_initial = α_initial,
#                             verbose_flag = verbose_flag,
#                             max_iter_tCG = max_iter_tCG,
#                             ρ_lower_acceptance = ρ_lower_acceptance,
#                             ρ_upper_acceptance = ρ_upper_acceptance,
#                             minimum_TR_radius = minimum_TR_radius,
#                             maximum_TR_radius = maximum_TR_radius,
#                             norm_df_tol = norm_df_tol,
#                             objective_tol = objective_tol,
#                             avg_Δf_tol = avg_Δf_tol,
#                             avg_Δf_window = avg_Δf_window,
#                             max_idle_update_count = max_idle_update_count,
#                             𝑟 = 𝑟,
#                             g = g)

#     α = clamp.(α_solution, zero_tol, Inf)
#     η = RKHSProblemType(α, X, θ, μ)

#     return η, α
# end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
# Some rough definitions of some common traits/interfaces of the
# standard library.  Note that most are single parameter traits.

export Eq, Cmp, IsIterable, IsCollection, IsIterColl, IsIndexable, IsAssociative, Arith

@traitdef Eq{X,Y} begin
    # note anything is part of Eq as ==(::Any,::Any) is defined
    ==(X,Y) -> Bool
end

@traitdef Cmp{X,Y} <: Eq{X,Y} begin
    isless(X,Y) -> Bool
    # automatically provides:
    # >, <, >=, <=
end

@traitdef IsIterable{X}  begin
    # type-functions based on return_type:
    State = Base.return_types(start, (X,))[1]
    Item =  Base.return_types(next, (X,State))[1][1] # use eltype instead

    # interface functions
    start(X) -> State
    next(X, State) -> Item, State
    done(X, State) -> Bool
    # automatically provides:
    # zip, enumerated, in, map, reduce, ...
end

# general collections
@traitdef IsCollection{X} <: IsIterable{X} begin
    El = eltype(X)

    isempty(X) -> Bool
    length(X) -> Integer
    eltype(X) -> Type{El}
    eltype(Type{X}) -> Type{El}
end

@traitdef IsIterColl{X} <: IsCollection{X} begin # iterable collection
    empty!(X) # -> X # ToDo: fix after updating return types
end

@traitdef IsIndexable{X} <:IsCollection{X} begin
    El = eltype(X) # note, for Associative this would be a tuple, which is wrong: see base_fixes.jl
    I = indextype(X)

    getindex(X, I) -> El
    setindex!(X, El, I)

    length(X) -> Integer
    # automatically provided:
    # size(X) -> Tuple
    # size(X,Integer) -> Integer
end

@traitdef IsAssociative{X} <: IsIndexable{X} begin
    V = eltype(X)
    K = indextype(X)

    # note, ObjectId dict is not part of this interface
    haskey(X, Any)
    get(X, Any, Any)
    get(Function, X, Any)
    get!(X, Any, Any)
    get!(Function, X, Any)
    getkey(X, Any, Any)
    delete!(X, Any) -> X
    pop!(X, Any)
    pop!(X, Any, Any)
    # merge(X, Any...) -> X
    # merge!(X, Any...)
    # provieds
    # keys(X) -> Base.KeyIterator
    # values(X) -> Base.ValueIterator
end

@traitdef Arith{X,Y} begin
    Z = promote_type(X,Y)
    D = Tuple{X,Y}<:Tuple{Integer, Integer} ? Float64 : Z

     # note, promote_type is defined for (Type{Any},Type{Any}), so no
     # need to include it here:
    +(X,Y) -> Z
    -(X,Y) -> Z
    *(X,Y) -> Z
    /(X,Y) -> D
end


text:
module OmegaCore

using Reexport
using Spec

include("util/util.jl")         # General utilities
@reexport using .Util

include("traits.jl")
@reexport using .Traits

include("tagging/tagging.jl")      # Tags
using .Tagging

include("rng.jl")               # Random number generation
using ..RNG

include("ids/ids.jl")           # IDs
@reexport using .IDS

include("basis/basis.jl")
@reexport using .Basis

include("var/var.jl")           # Random / Parameteric Variables
@reexport using .Var

include("space/space.jl")       # Probability / Paramter Spaces
@reexport using .Space

# include("ciid.jl")            # Conditional Independence
# @reexport using .CIID

include("interventions/interventions.jl")         # Causal interventions
@reexport using .Interventions

include("rand.jl")            # Higher order inference
@reexport using .Rand

# include("cassette.jl")

include("conditioning.jl")         # Conditioning variables
@reexport using .Conditioning

include("Higher/higher.jl")            # Higher order inference
@reexport using .Higher

include("queries.jl")           # Query Templtes
@reexport using .Queries

include("sample.jl")            # Sample
@reexport using .Sample

include("trackerror.jl")
@reexport using .TrackError

include("proposal/proposal.jl")            # Log density
@reexport using .Proposal

include("solve.jl")               # Satisfy
@reexport using .Solver

include("seen.jl")
@reexport using .SeenVars

include("logenergy.jl")
@reexport using .LogEnergies

# Basic Inference methods
# include("fail.jl")              # Fails when conditions are not satisfied
# @reexport using .Fail

include("rejection.jl")         # Rejection sampling Inference
@reexport using .OmegaRejectionSample

# include("pointwise.jl")
# @reexport using .Pointwise

include("mem.jl")
@reexport using .Memoize

include("syntax/syntax.jl")
@reexport using .Syntax

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect