text:
module TestGrouping

using Test, DataFrames, Random, Statistics
const ≅ = isequal

# Call groupby, checking that groups field is consistent with other fields
# (since == and isequal do not use it)
function groupby_checked(df::AbstractDataFrame, keys, args...; kwargs...)
    gd = groupby(df, keys, args...; kwargs...)
    for i in 1:length(gd)
        @assert findall(==(i), gd.groups) == gd.idx[gd.starts[i]:gd.ends[i]]
    end
    gd
end

@testset "parent" begin
    df = DataFrame(a = [1, 1, 2, 2], b = [5, 6, 7, 8])
    gd = groupby(df, :a)
    @test parent(gd) === df
end

@testset "colwise" begin
    Random.seed!(1)
    df = DataFrame(a = repeat(Union{Int, Missing}[1, 3, 2, 4], outer=[2]),
                   b = repeat(Union{Int, Missing}[2, 1], outer=[4]),
                   c = Vector{Union{Float64, Missing}}(randn(8)))

    missingfree = DataFrame([collect(1:10)], [:x1])

    @testset "::Function, ::AbstractDataFrame" begin
        cw = colwise(sum, df)
        answer = [20, 12, -0.4283098098931877]
        @test isa(cw, Vector{Real})
        @test size(cw) == (ncol(df),)
        @test cw == answer

        cw = colwise(sum, missingfree)
        answer = [55]
        @test isa(cw, Array{Int, 1})
        @test size(cw) == (ncol(missingfree),)
        @test cw == answer
    end

    @testset "::Function, ::GroupedDataFrame" begin
        gd = groupby_checked(DataFrame(A = [:A, :A, :B, :B], B = 1:4), :A)
        @test colwise(length, gd) == [[2,2], [2,2]]
    end

    @testset "::Vector, ::AbstractDataFrame" begin
        cw = colwise([sum], df)
        answer = [20 12 -0.4283098098931877]
        @test isa(cw, Array{Real, 2})
        @test size(cw) == (length([sum]),ncol(df))
        @test cw == answer

        cw = colwise([sum, minimum], missingfree)
        answer = reshape([55, 1], (2,1))
        @test isa(cw, Array{Int, 2})
        @test size(cw) == (length([sum, minimum]), ncol(missingfree))
        @test cw == answer

        cw = colwise([Vector{Union{Int, Missing}}], missingfree)
        answer = reshape([Vector{Union{Int, Missing}}(1:10)], (1,1))
        @test isa(cw, Array{Vector{Union{Int, Missing}},2})
        @test size(cw) == (1, ncol(missingfree))
        @test cw == answer

        @test_throws MethodError colwise(["Bob", :Susie], DataFrame(A = 1:10, B = 11:20))
    end

    @testset "::Vector, ::GroupedDataFrame" begin
        gd = groupby_checked(DataFrame(A = [:A, :A, :B, :B], B = 1:4), :A)
        @test colwise([length], gd) == [[2 2], [2 2]]
    end

    @testset "::Tuple, ::AbstractDataFrame" begin
        cw = colwise((sum, length), df)
        answer = Any[20 12 -0.4283098098931877; 8 8 8]
        @test isa(cw, Array{Real, 2})
        @test size(cw) == (length((sum, length)), ncol(df))
        @test cw == answer

        cw = colwise((sum, length), missingfree)
        answer = reshape([55, 10], (2,1))
        @test isa(cw, Array{Int, 2})
        @test size(cw) == (length((sum, length)), ncol(missingfree))
        @test cw == answer

        cw = colwise((CategoricalArray, Vector{Union{Int, Missing}}), missingfree)
        answer = reshape([CategoricalArray(1:10), Vector{Union{Int, Missing}}(1:10)],
                         (2, ncol(missingfree)))
        @test typeof(cw) == Array{AbstractVector,2}
        @test size(cw) == (2, ncol(missingfree))
        @test cw == answer

        @test_throws MethodError colwise(("Bob", :Susie), DataFrame(A = 1:10, B = 11:20))
    end

    @testset "::Tuple, ::GroupedDataFrame" begin
        gd = groupby_checked(DataFrame(A = [:A, :A, :B, :B], B = 1:4), :A)
        @test colwise((length), gd) == [[2,2],[2,2]]
    end
end

@testset "by, groupby and map(::Function, ::GroupedDataFrame)" begin
    Random.seed!(1)
    df = DataFrame(a = repeat(Union{Int, Missing}[1, 3, 2, 4], outer=[2]),
                   b = repeat(Union{Int, Missing}[2, 1], outer=[4]),
                   c = repeat([0, 1], outer=[4]),
                   x = Vector{Union{Float64, Missing}}(randn(8)))

    f1(df) = DataFrame(xmax = maximum(df.x))
    f2(df) = (xmax = maximum(df.x),)
    f3(df) = maximum(df.x)
    f4(df) = [maximum(df.x), minimum(df.x)]
    f5(df) = reshape([maximum(df.x), minimum(df.x)], 2, 1)
    f6(df) = [maximum(df.x) minimum(df.x)]
    f7(df) = (x2 = df.x.^2,)
    f8(df) = DataFrame(x2 = df.x.^2)

    for cols in ([:a, :b], [:b, :a], [:a, :c], [:c, :a],
                 [1, 2], [2, 1], [1, 3], [3, 1],
                 [true, true, false, false], [true, false, true, false])
        colssym = names(df[cols])
        hcatdf = hcat(df[cols], df, makeunique=true)
        nms = names(hcatdf)
        res = unique(df[cols])
        res.xmax = [maximum(df[(df[colssym[1]] .== a) .& (df[colssym[2]] .== b), :x])
                    for (a, b) in zip(res[colssym[1]], res[colssym[2]])]
        res2 = unique(df[cols])[repeat(1:4, inner=2), :]
        res2.x1 = collect(Iterators.flatten(
            [[maximum(df[(df[colssym[1]] .== a) .& (df[colssym[2]] .== b), :x]),
              minimum(df[(df[colssym[1]] .== a) .& (df[colssym[2]] .== b), :x])]
             for (a, b) in zip(res[colssym[1]], res[colssym[2]])]))
        res3 = unique(df[cols])
        res3.x1 = [maximum(df[(df[colssym[1]] .== a) .& (df[colssym[2]] .== b), :x])
                   for (a, b) in zip(res[colssym[1]], res[colssym[2]])]
        res3.x2 = [minimum(df[(df[colssym[1]] .== a) .& (df[colssym[2]] .== b), :x])
                   for (a, b) in zip(res[colssym[1]], res[colssym[2]])]
        res4 = df[cols]
        res4.x2 = df.x.^2
        shcatdf = sort(hcatdf, colssym)
        sres = sort(res, colssym)
        sres2 = sort(res2, colssym)
        sres3 = sort(res3, colssym)
        sres4 = sort(res4, colssym)

        # by() without groups sorting
        @test sort(by(df, cols, identity), colssym) == shcatdf
        @test sort(by(df, cols, df -> df[1, :]), colssym) ==
            shcatdf[.!nonunique(shcatdf, colssym), :]
        @test by(df, cols, f1) == res
        @test by(df, cols, f2) == res
        @test rename(by(df, cols, f3), :x1 => :xmax) == res
        @test by(df, cols, f4) == res2
        @test by(df, cols, f5) == res2
        @test by(df, cols, f6) == res3
        @test sort(by(df, cols, f7), colssym) == sres4
        @test sort(by(df, cols, f8), colssym) == sres4

        # by() with groups sorting
        @test by(df, cols, identity, sort=true) == shcatdf
        @test by(df, cols, df -> df[1, :], sort=true) ==
            shcatdf[.!nonunique(shcatdf, colssym), :]
        @test by(df, cols, f1, sort=true) == sres
        @test by(df, cols, f2, sort=true) == sres
        @test rename(by(df, cols, f3, sort=true), :x1 => :xmax) == sres
        @test by(df, cols, f4, sort=true) == sres2
        @test by(df, cols, f5, sort=true) == sres2
        @test by(df, cols, f6, sort=true) == sres3
        @test by(df, cols, f7, sort=true) == sres4
        @test by(df, cols, f8, sort=true) == sres4

        @test by(df, [:a], f1) == by(df, :a, f1)
        @test by(df, [:a], f1, sort=true) == by(df, :a, f1, sort=true)

        # groupby() without groups sorting
        gd = groupby_checked(df, cols)
        @test names(parent(gd))[gd.cols] == colssym
        @test sort(combine(identity, gd), colssym) ==
            sort(combine(gd), colssym) ==
            shcatdf
        @test combine(f1, gd) == res
        @test combine(f2, gd) == res
        @test rename(combine(f3, gd), :x1 => :xmax) == res
        @test combine(f4, gd) == res2
        @test combine(f5, gd) == res2
        @test combine(f6, gd) == res3
        @test sort(combine(f7, gd), colssym) == sort(res4, colssym)
        @test sort(combine(f8, gd), colssym) == sort(res4, colssym)

        # groupby() with groups sorting
        gd = groupby_checked(df, cols, sort=true)
        @test names(parent(gd))[gd.cols] == colssym
        for i in 1:length(gd)
            @test all(gd[i][colssym[1]] .== sres[i, colssym[1]])
            @test all(gd[i][colssym[2]] .== sres[i, colssym[2]])
        end
        @test combine(identity, gd) == combine(gd) == shcatdf
        @test combine(f1, gd) == sres
        @test combine(f2, gd) == sres
        @test rename(combine(f3, gd), :x1 => :xmax) == sres
        @test combine(f4, gd) == sres2
        @test combine(f5, gd) == sres2
        @test combine(f6, gd) == sres3
        @test combine(f7, gd) == sres4
        @test combine(f8, gd) == sres4

        # map() without and with groups sorting
        for sort in (false, true)
            gd = groupby_checked(df, cols, sort=sort)
            v = map(d -> d[[:x]], gd)
            @test length(gd) == length(v)
            nms = [colssym; :x]
            @test v[1] == gd[1][nms]
            @test v[1] == gd[1][nms] &&
                v[2] == gd[2][nms] &&
                v[3] == gd[3][nms] &&
                v[4] == gd[4][nms]
            @test names(parent(v))[v.cols] == colssym
            v = map(f1, gd)
            @test vcat(v[1], v[2], v[3], v[4]) == by(f1, df, cols, sort=sort)
            v = map(f2, gd)
            @test vcat(v[1], v[2], v[3], v[4]) == by(f2, df, cols, sort=sort)
            v = map(f3, gd)
            @test vcat(v[1], v[2], v[3], v[4]) == by(f3, df, cols, sort=sort)
            v = map(f4, gd)
            @test vcat(v[1], v[2], v[3], v[4]) == by(f4, df, cols, sort=sort)
            v = map(f5, gd)
            @test vcat(v[1], v[2], v[3], v[4]) == by(f5, df, cols, sort=sort)
            v = map(f5, gd)
            @test vcat(v[1], v[2], v[3], v[4]) == by(f5, df, cols, sort=sort)
            v = map(f6, gd)
            @test vcat(v[1], v[2], v[3], v[4]) == by(f6, df, cols, sort=sort)
            v = map(f7, gd)
            @test vcat(v[1], v[2], v[3], v[4]) == by(f7, df, cols, sort=sort)
            v = map(f8, gd)
            @test vcat(v[1], v[2], v[3], v[4]) == by(f8, df, cols, sort=sort)
        end
    end

    # test number of potential combinations higher than typemax(Int32)
    N = 2000
    df2 = DataFrame(v1 = levels!(categorical(rand(1:N, 100)), collect(1:N)),
                    v2 = levels!(categorical(rand(1:N, 100)), collect(1:N)),
                    v3 = levels!(categorical(rand(1:N, 100)), collect(1:N)))
    df2b = mapcols(Vector{Int}, df2)
    @test groupby_checked(df2, [:v1, :v2, :v3]) ==
        groupby_checked(df2b, [:v1, :v2, :v3])

    # grouping empty table
    @test groupby_checked(DataFrame(A=Int[]), :A).starts == Int[]
    # grouping single row
    @test groupby_checked(DataFrame(A=Int[1]), :A).starts == Int[1]

    # issue #960
    x = CategoricalArray(collect(1:20))
    df = DataFrame(v1=x, v2=x)
    groupby_checked(df, [:v1, :v2])

    df2 = by(e->1, DataFrame(x=Int64[]), :x)
    @test size(df2) == (0, 1)
    @test sum(df2[:x]) == 0

    # Check that reordering levels does not confuse groupby
    df = DataFrame(Key1 = CategoricalArray(["A", "A", "B", "B", "B", "A"]),
                   Key2 = CategoricalArray(["A", "B", "A", "B", "B", "A"]),
                   Value = 1:6)
    gd = groupby_checked(df, :Key1)
    @test length(gd) == 2
    @test gd[1] == DataFrame(Key1="A", Key2=["A", "B", "A"], Value=[1, 2, 6])
    @test gd[2] == DataFrame(Key1="B", Key2=["A", "B", "B"], Value=[3, 4, 5])
    gd = groupby_checked(df, [:Key1, :Key2])
    @test length(gd) == 4
    @test gd[1] == DataFrame(Key1="A", Key2="A", Value=[1, 6])
    @test gd[2] == DataFrame(Key1="A", Key2="B", Value=2)
    @test gd[3] == DataFrame(Key1="B", Key2="A", Value=3)
    @test gd[4] == DataFrame(Key1="B", Key2="B", Value=[4, 5])
    # Reorder levels, add unused level
    levels!(df[:Key1], ["Z", "B", "A"])
    levels!(df[:Key2], ["Z", "B", "A"])
    gd = groupby_checked(df, :Key1)
    @test gd == groupby_checked(df, :Key1, skipmissing=true)
    @test length(gd) == 2
    @test gd[1] == DataFrame(Key1="B", Key2=["A", "B", "B"], Value=[3, 4, 5])
    @test gd[2] == DataFrame(Key1="A", Key2=["A", "B", "A"], Value=[1, 2, 6])
    gd = groupby_checked(df, [:Key1, :Key2])
    @test gd == groupby_checked(df, [:Key1, :Key2], skipmissing=true)
    @test length(gd) == 4
    @test gd[1] == DataFrame(Key1="B", Key2="B", Value=[4, 5])
    @test gd[2] == DataFrame(Key1="B", Key2="A", Value=3)
    @test gd[3] == DataFrame(Key1="A", Key2="B", Value=2)
    @test gd[4] == DataFrame(Key1="A", Key2="A", Value=[1, 6])
    # Make first level unused too
    replace!(df.Key1, "A"=>"B")
    gd = groupby_checked(df, :Key1)
    @test length(gd) == 1
    @test gd[1] == DataFrame(Key1="B", Key2=["A", "B", "A", "B", "B", "A"], Value=1:6)
    gd = groupby_checked(df, [:Key1, :Key2])
    @test gd == groupby_checked(df, [:Key1, :Key2])
    @test length(gd) == 2
    @test gd[1] == DataFrame(Key1="B", Key2="B", Value=[2, 4, 5])
    @test gd[2] == DataFrame(Key1="B", Key2="A", Value=[1, 3, 6])

    # Check that CategoricalArray column is preserved when returning a value...
    res = combine(d -> DataFrame(x=d[1, :Key2]), groupby_checked(df, :Key1))
    @test res.x isa CategoricalVector{String}
    res = combine(d -> (x=d[1, :Key2],), groupby_checked(df, :Key1))
    @test res.x isa CategoricalVector{String}
    # ...and when returning an array
    res = combine(d -> DataFrame(x=d[:Key1]), groupby_checked(df, :Key1))
    @test res.x isa CategoricalVector{String}

    # Check that CategoricalArray and String give a String...
    res = combine(d -> d.Key1 == ["A", "A"] ? DataFrame(x=d[1, :Key1]) : DataFrame(x="C"),
                 groupby_checked(df, :Key1))
    @test res.x isa Vector{String}
    res = combine(d -> d.Key1 == ["A", "A"] ? (x=d[1, :Key1],) : (x="C",),
                  groupby_checked(df, :Key1))
    @test res.x isa Vector{String}
    # ...even when CategoricalString comes second
    res = combine(d -> d.Key1 == ["B", "B"] ? DataFrame(x=d[1, :Key1]) : DataFrame(x="C"),
                  groupby_checked(df, :Key1))
    @test res.x isa Vector{String}
    res = combine(d -> d.Key1 == ["B", "B"] ? (x=d[1, :Key1],) : (x="C",),
                  groupby_checked(df, :Key1))
    @test res.x isa Vector{String}

    df = DataFrame(x = [1, 2, 3], y = [2, 3, 1])

    # Test function returning DataFrameRow
    res = by(d -> DataFrameRow(d, 1, :), df, :x)
    @test res == DataFrame(x=df.x, x_1=df.x, y=df.y)

    # Test function returning Tuple
    res = by(d -> (sum(d.y),), df, :x)
    @test res == DataFrame(x=df.x, x1=tuple.([2, 3, 1]))

    # Test with some groups returning empty data frames
    @test by(d -> d.x == [1] ? DataFrame(z=[]) : DataFrame(z=1), df, :x) ==
        DataFrame(x=[2, 3], z=[1, 1])
    v = map(d -> d.x == [1] ? DataFrame(z=[]) : DataFrame(z=1), groupby_checked(df, :x))
    @test length(v) == 2
    @test vcat(v[1], v[2]) == DataFrame(x=[2, 3], z=[1, 1])

    # Test that returning values of different types works with NamedTuple
    res = by(d -> d.x == [1] ? 1 : 2.0, df, :x)
    @test res.x1 isa Vector{Float64}
    @test res.x1 == [1, 2, 2]
    # Two columns need to be widened at different times
    res = by(d -> (a=d.x == [1] ? 1 : 2.0, b=d.x == [3] ? missing : "a"), df, :x)
    @test res.a isa Vector{Float64}
    @test res.a == [1, 2, 2]
    @test res.b isa Vector{Union{String,Missing}}
    @test res.b ≅ ["a", "a", missing]
    # Corner case: two columns need to be widened at the same time
    res = by(d -> (a=d.x == [1] ? 1 : 2.0, b=d.x == [1] ? missing : "a"), df, :x)
    @test res.a isa Vector{Float64}
    @test res.a == [1, 2, 2]
    @test res.b isa Vector{Union{String,Missing}}
    @test res.b ≅ [missing, "a", "a"]

    # Test that returning values of different types works with DataFrame
    res = by(d -> DataFrame(x1 = d.x == [1] ? 1 : 2.0), df, :x)
    @test res.x1 isa Vector{Float64}
    @test res.x1 == [1, 2, 2]
    # Two columns need to be widened at different times
    res = by(d -> DataFrame(a=d.x == [1] ? 1 : 2.0, b=d.x == [3] ? missing : "a"), df, :x)
    @test res.a isa Vector{Float64}
    @test res.a == [1, 2, 2]
    @test res.b isa Vector{Union{String,Missing}}
    @test res.b ≅ ["a", "a", missing]
    # Corner case: two columns need to be widened at the same time
    res = by(d -> DataFrame(a=d.x == [1] ? 1 : 2.0, b=d.x == [1] ? missing : "a"), df, :x)
    @test res.a isa Vector{Float64}
    @test res.a == [1, 2, 2]
    @test res.b isa Vector{Union{String,Missing}}
    @test res.b ≅ [missing, "a", "a"]

    # Test return values with columns in different orders
    @test by(d -> d.x == [1] ? (x1=1, x2=3) : (x2=2, x1=4), df, :x) ==
        DataFrame(x=1:3, x1=[1, 4, 4], x2=[3, 2, 2])
    @test by(d -> d.x == [1] ? DataFrame(x1=1, x2=3) : DataFrame(x2=2, x1=4), df, :x) ==
        DataFrame(x=1:3, x1=[1, 4, 4], x2=[3, 2, 2])

    # Test with NamedTuple with columns of incompatible lengths
    @test_throws DimensionMismatch by(d -> (x1=[1], x2=[3, 4]), df, :x)
    @test_throws DimensionMismatch by(d -> d.x == [1] ? (x1=[1], x2=[3]) :
                                                        (x1=[1], x2=[3, 4]), df, :x)

    # Test with incompatible return values
    @test_throws ArgumentError by(d -> d.x == [1] ? (x1=1,) : DataFrame(x1=1), df, :x)
    @test_throws ArgumentError by(d -> d.x == [1] ? NamedTuple() : (x1=1), df, :x)
    @test_throws ArgumentError by(d -> d.x == [1] ? 1 : DataFrame(x1=1), df, :x)
    @test_throws ArgumentError by(d -> d.x == [1] ? DataFrame() : DataFrame(x1=1), df, :x)
    # Special case allowed due to how implementation works
    @test by(d -> d.x == [1] ? 1 : (x1=1), df, :x) == by(d -> 1, df, :x)

    # Test that columns names and types are respected for empty input
    df = DataFrame(x=Int[], y=String[])
    res = by(d -> 1, df, :x)
    @test size(res) == (0, 1)
    @test res.x isa Vector{Int}

    # Test with empty data frame
    df = DataFrame(x=[], y=[])
    gd = groupby_checked(df, :x)
    @test combine(df -> sum(df.x), gd) == DataFrame(x=[])
    res = map(df -> sum(df.x), gd)
    @test length(res) == 0
    @test res.parent == DataFrame(x=[])
end

xv = ["A", missing, "B", "B", "A", "B", "A", "A"]
yv = ["B", "A", "A", missing, "A", missing, "A", "A"]
@testset "grouping with missings" for x in (xv,
                                            categorical(xv),
                                            levels!(categorical(xv), ["A", "B", "X"]),
                                            levels!(categorical(xv), ["X", "A", "B"])),
                                      y in (yv,
                                            categorical(yv),
                                            levels!(categorical(yv), ["A", "B", "X"]),
                                            levels!(categorical(yv), ["X", "A", "B"]))
    df = DataFrame(Key1 = x, Key2 = y, Value = 1:8)

    @testset "sort=false, skipmissing=false" begin
        gd = groupby_checked(df, :Key1)
        @test length(gd) == 3
        if df.Key1 isa CategoricalVector # Order of missing changes
            @test gd[1] == DataFrame(Key1="A", Key2=["B", "A", "A", "A"], Value=[1, 5, 7, 8])
            @test gd[2] ≅ DataFrame(Key1="B", Key2=["A", missing, missing], Value=[3, 4, 6])
            @test gd[3] ≅ DataFrame(Key1=missing, Key2="A", Value=2)
        else
            @test gd[1] == DataFrame(Key1="A", Key2=["B", "A", "A", "A"], Value=[1, 5, 7, 8])
            @test gd[2] ≅ DataFrame(Key1=missing, Key2="A", Value=2)
            @test gd[3] ≅ DataFrame(Key1="B", Key2=["A", missing, missing], Value=[3, 4, 6])
        end

        gd = groupby_checked(df, [:Key1, :Key2])
        @test length(gd) == 5
        if df.Key1 isa CategoricalVector && df.Key2 isa CategoricalVector
            @test gd[1] ≅ DataFrame(Key1="A", Key2="A", Value=[5, 7, 8])
            @test gd[2] == DataFrame(Key1="A", Key2="B", Value=1)
            @test gd[3] == DataFrame(Key1="B", Key2="A", Value=3)
            @test gd[4] ≅ DataFrame(Key1="B", Key2=missing, Value=[4, 6])
            @test gd[5] ≅ DataFrame(Key1=missing, Key2="A", Value=2)
        else
            @test gd[1] == DataFrame(Key1="A", Key2="B", Value=1)
            @test gd[2] ≅ DataFrame(Key1=missing, Key2="A", Value=2)
            @test gd[3] == DataFrame(Key1="B", Key2="A", Value=3)
            @test gd[4] ≅ DataFrame(Key1="B", Key2=missing, Value=[4, 6])
            @test gd[5] ≅ DataFrame(Key1="A", Key2="A", Value=[5, 7, 8])
        end
    end

    @testset "sort=false, skipmissing=true" begin
        gd = groupby_checked(df, :Key1, skipmissing=true)
        @test length(gd) == 2
        @test gd[1] == DataFrame(Key1="A", Key2=["B", "A", "A", "A"], Value=[1, 5, 7, 8])
        @test gd[2] ≅ DataFrame(Key1="B", Key2=["A", missing, missing], Value=[3, 4, 6])

        gd = groupby_checked(df, [:Key1, :Key2], skipmissing=true)
        @test length(gd) == 3
        if df.Key1 isa CategoricalVector && df.Key2 isa CategoricalVector
            @test gd[1] == DataFrame(Key1="A", Key2="A", Value=[5, 7, 8])
            @test gd[2] == DataFrame(Key1="A", Key2="B", Value=1)
            @test gd[3] == DataFrame(Key1="B", Key2="A", Value=3)
        else
            @test gd[1] == DataFrame(Key1="A", Key2="B", Value=1)
            @test gd[2] == DataFrame(Key1="B", Key2="A", Value=3)
            @test gd[3] == DataFrame(Key1="A", Key2="A", Value=[5, 7, 8])
        end
    end

    @testset "sort=true, skipmissing=false" begin
        gd = groupby_checked(df, :Key1, sort=true)
        @test length(gd) == 3
        @test gd[1] == DataFrame(Key1="A", Key2=["B", "A", "A", "A"], Value=[1, 5, 7, 8])
        @test gd[2] ≅ DataFrame(Key1="B", Key2=["A", missing, missing], Value=[3, 4, 6])
        @test gd[3] ≅ DataFrame(Key1=missing, Key2="A", Value=2)

        gd = groupby_checked(df, [:Key1, :Key2], sort=true)
        @test length(gd) == 5
        @test gd[1] ≅ DataFrame(Key1="A", Key2="A", Value=[5, 7, 8])
        @test gd[2] == DataFrame(Key1="A", Key2="B", Value=1)
        @test gd[3] == DataFrame(Key1="B", Key2="A", Value=3)
        @test gd[4] ≅ DataFrame(Key1="B", Key2=missing, Value=[4, 6])
        @test gd[5] ≅ DataFrame(Key1=missing, Key2="A", Value=2)
    end

    @testset "sort=true, skipmissing=true" begin
        gd = groupby_checked(df, :Key1, sort=true, skipmissing=true)
        @test length(gd) == 2
        @test gd[1] == DataFrame(Key1="A", Key2=["B", "A", "A", "A"], Value=[1, 5, 7, 8])
        @test gd[2] ≅ DataFrame(Key1="B", Key2=["A", missing, missing], Value=[3, 4, 6])

        gd = groupby_checked(df, [:Key1, :Key2], sort=true, skipmissing=true)
        @test length(gd) == 3
        @test gd[1] == DataFrame(Key1="A", Key2="A", Value=[5, 7, 8])
        @test gd[2] == DataFrame(Key1="A", Key2="B", Value=1)
        @test gd[3] == DataFrame(Key1="B", Key2="A", Value=3)
    end
end

# We need many rows so that optimized CategoricalArray method is used
xv = rand(["A", "B", missing], 100)
yv = rand(["A", "B", missing], 100)
zv = rand(["A", "B", missing], 100)
@testset "grouping with three keys" for x in (xv,
                                              categorical(xv),
                                              levels!(categorical(xv), ["A", "B", "X"]),
                                              levels!(categorical(xv), ["X", "A", "B"])),
                                        y in (yv,
                                              categorical(yv),
                                              levels!(categorical(yv), ["A", "B", "X"]),
                                              levels!(categorical(yv), ["X", "A", "B"])),
                                        z in (zv,
                                              categorical(zv),
                                              levels!(categorical(zv), ["A", "B", "X"]),
                                              levels!(categorical(zv), ["X", "A", "B"]))
    df = DataFrame(Key1 = x, Key2 = y, Key3 = z, Value = string.(1:100))
    dfb = mapcols(Vector{Union{String, Missing}}, df)

    @test groupby_checked(df, [:Key1, :Key2, :Key3], sort=true) ≅
        groupby_checked(dfb, [:Key1, :Key2, :Key3], sort=true)
    @test groupby_checked(df, [:Key1, :Key2, :Key3], sort=true, skipmissing=true) ≅
          groupby_checked(dfb, [:Key1, :Key2, :Key3], sort=true, skipmissing=true)

    if df.Key1 isa CategoricalVector &&
       df.Key2 isa CategoricalVector &&
       df.Key3 isa CategoricalVector
        @test groupby_checked(df, [:Key1, :Key2, :Key3], sort=true) ≅
            groupby_checked(df, [:Key1, :Key2, :Key3], sort=false)
        @test groupby_checked(df, [:Key1, :Key2, :Key3], sort=true, skipmissing=true) ≅
            groupby_checked(df, [:Key1, :Key2, :Key3], sort=false, skipmissing=true)
    end
end

@testset "by, combine and map with pair interface" begin
    vexp = x -> exp.(x)
    Random.seed!(1)
    df = DataFrame(a = repeat([1, 3, 2, 4], outer=[2]),
                   b = repeat([2, 1], outer=[4]),
                   c = rand(Int, 8))

    # Only test that different by syntaxes work,
    # and rely on tests below for deeper checks
    @test by(:c => sum, df, :a) ==
        by(df, :a, :c => sum) ==
        by(df, :a, (:c => sum,)) ==
        by(df, :a, [:c => sum]) ==
        by(df, :a, c_sum = :c => sum) ==
        by(d -> (c_sum=sum(d.c),), df, :a)
        by(df, :a, d -> (c_sum=sum(d.c),))

    @test by(:c => vexp, df, :a) ==
        by(df, :a, :c => vexp) ==
        by(df, :a, (:c => vexp,)) ==
        by(df, :a, [:c => vexp]) ==
        by(df, :a, c_function = :c => vexp) ==
        by(d -> (c_function=vexp(d.c),), df, :a)
        by(df, :a, d -> (c_function=vexp(d.c),))

    @test by(df, :a, :b => sum, :c => sum) ==
        by(df, :a, (:b => sum, :c => sum,)) ==
        by(df, :a, [:b => sum, :c => sum]) ==
        by(df, :a, b_sum = :b => sum, c_sum = :c => sum) ==
        by(d -> (b_sum=sum(d.b), c_sum=sum(d.c)), df, :a)
        by(df, :a, d -> (b_sum=sum(d.b), c_sum=sum(d.c)))

    @test by(df, :a, :b => vexp, :c => identity) ==
        by(df, :a, (:b => vexp, :c => identity,)) ==
        by(df, :a, [:b => vexp, :c => identity]) ==
        by(df, :a, b_function = :b => vexp, c_identity = :c => identity) ==
        by(d -> (b_function=vexp(d.b), c_identity=identity(d.c)), df, :a)
        by(df, :a, d -> (b_function=vexp(d.b), c_identity=identity(d.c)))

    gd = groupby(df, :a)

    # Only test that different combine syntaxes work,
    # and rely on tests below for deeper checks
    @test combine(:c => sum, gd) ==
        combine(gd, :c => sum) ==
        combine(gd, (:c => sum,)) ==
        combine(gd, [:c => sum]) ==
        combine(gd, c_sum = :c => sum) ==
        combine(:c => x -> (c_sum=sum(x),), gd) ==
        combine(gd, :c => x -> (c_sum=sum(x),)) ==
        combine(d -> (c_sum=sum(d.c),), gd) ==
        combine(gd, d -> (c_sum=sum(d.c),))

    @test combine(:c => vexp, gd) ==
        combine(gd, :c => vexp) ==
        combine(gd, (:c => vexp,)) ==
        combine(gd, [:c => vexp]) ==
        combine(gd, c_function = :c => vexp) ==
        combine(:c => x -> (c_function=exp.(x),), gd) ==
        combine(gd, :c => x -> (c_function=exp.(x),)) ==
        combine(d -> (c_function=exp.(d.c),), gd)
        combine(gd, d -> (c_function=exp.(d.c),))

    @test combine(gd, :b => sum, :c => sum) ==
        combine(gd, (:b => sum, :c => sum,)) ==
        combine(gd, [:b => sum, :c => sum]) ==
        combine(gd, b_sum = :b => sum, c_sum = :c => sum) ==
        combine((:b, :c) => x -> (b_sum=sum(x.b), c_sum=sum(x.c)), gd) ==
        combine(gd, (:b, :c) => x -> (b_sum=sum(x.b), c_sum=sum(x.c))) ==
        combine(d -> (b_sum=sum(d.b), c_sum=sum(d.c)), gd) ==
        combine(gd, d -> (b_sum=sum(d.b), c_sum=sum(d.c)))

    @test combine(gd, :b => vexp, :c => identity) ==
        combine(gd, (:b => vexp, :c => identity,)) ==
        combine(gd, [:b => vexp, :c => identity]) ==
        combine(gd, b_function = :b => vexp, c_identity = :c => identity) ==
        combine((:b, :c) => x -> (b_function=vexp(x.b), c_identity=x.c), gd) ==
        combine(gd, (:b, :c) => x -> (b_function=vexp(x.b), c_identity=x.c)) ==
        combine(d -> (b_function=vexp(d.b), c_identity=d.c), gd) ==
        combine(gd, d -> (b_function=vexp(d.b), c_identity=d.c))

    for f in (map, combine)
        for col in (:c, 3)
            @test f(col => sum, gd) == f(d -> (c_sum=sum(d.c),), gd)
            @test f(col => x -> sum(x), gd) == f(d -> (c_function=sum(d.c),), gd)
            @test f(col => x -> (z=sum(x),), gd) == f(d -> (z=sum(d.c),), gd)
            @test f(col => x -> DataFrame(z=sum(x),), gd) == f(d -> (z=sum(d.c),), gd)
            @test f(col => identity, gd) == f(d -> (c_identity=d.c,), gd)
            @test f(col => x -> (z=x,), gd) == f(d -> (z=d.c,), gd)

            @test f((xyz = col => sum,), gd) ==
                f(d -> (xyz=sum(d.c),), gd)
            @test f((xyz = col => x -> sum(x),), gd) ==
                f(d -> (xyz=sum(d.c),), gd)
            @test f((xyz = col => x -> (sum(x),),), gd) ==
                f(d -> (xyz=(sum(d.c),),), gd)
            @test_throws ArgumentError f((xyz = col => x -> (z=sum(x),),), gd)
            @test_throws ArgumentError f((xyz = col => x -> DataFrame(z=sum(x),),), gd)
            @test_throws ArgumentError f((xyz = col => x -> (z=x,),), gd)
            @test_throws ArgumentError f(col => x -> (z=1, xzz=[1]), gd)

            for wrap in (vcat, tuple)
                @test f(wrap(col => sum), gd) ==
                    f(d -> (c_sum=sum(d.c),), gd)
                @test f(wrap(col => x -> sum(x)), gd) ==
                    f(d -> (c_function=sum(d.c),), gd)
                @test f(wrap(col => x -> (sum(x),)), gd) ==
                    f(d -> (c_function=(sum(d.c),),), gd)
                @test_throws ArgumentError f(wrap(col => x -> (z=sum(x),)), gd)
                @test_throws ArgumentError f(wrap(col => x -> DataFrame(z=sum(x),)), gd)
                @test_throws ArgumentError f(wrap(col => x -> (z=x,)), gd)
                @test_throws ArgumentError f(wrap(col => x -> (z=1, xzz=[1])), gd)
            end
        end
        for cols in ((:b, :c), [:b, :c], (2, 3), 2:3, [2, 3], [false, true, true])
            @test f(cols => x -> (y=exp.(x.b), z=x.c), gd) ==
                f(d -> (y=exp.(d.b), z=d.c), gd)
            @test f(cols => x -> [exp.(x.b) x.c], gd) ==
                f(d -> [exp.(d.b) d.c], gd)

            @test f((xyz = cols => x -> sum(x.b) + sum(x.c),), gd) ==
                f(d -> (xyz=sum(d.b) + sum(d.c),), gd)
            if eltype(cols) === Bool
                cols2 = [[false, true, false], [false, false, true]]
                @test_throws MethodError f((xyz = cols[1] => sum, xzz = cols2[2] => sum), gd)
                @test_throws MethodError f((xyz = cols[1] => sum, xzz = cols2[1] => sum), gd)
                @test_throws MethodError f((xyz = cols[1] => sum, xzz = cols2[2] => x -> first(x)), gd)
            else
                cols2 = cols
                @test f((xyz = cols2[1] => sum, xzz = cols2[2] => sum), gd) ==
                    f(d -> (xyz=sum(d.b), xzz=sum(d.c)), gd)
                @test f((xyz = cols2[1] => sum, xzz = cols2[1] => sum), gd) ==
                    f(d -> (xyz=sum(d.b), xzz=sum(d.b)), gd)
                @test f((xyz = cols2[1] => sum, xzz = cols2[2] => x -> first(x)), gd) ==
                    f(d -> (xyz=sum(d.b), xzz=first(d.c)), gd)
                @test_throws ArgumentError f((xyz = cols2[1] => vexp, xzz = cols2[2] => sum), gd)
            end

            @test_throws ArgumentError f(cols => x -> (y=exp.(x.b), z=sum(x.c)), gd)
            @test_throws ArgumentError f((xyz = cols2 => x -> DataFrame(y=exp.(x.b), z=sum(x.c)),), gd)
            @test_throws ArgumentError f((xyz = cols2 => x -> [exp.(x.b) x.c],), gd)

            for wrap in (vcat, tuple)
                @test f(wrap(cols => x -> sum(x.b) + sum(x.c)), gd) ==
                    f(d -> sum(d.b) + sum(d.c), gd)

                if eltype(cols) === Bool
                    cols2 = [[false, true, false], [false, false, true]]
                    @test f(wrap(cols2[1] => x -> sum(x.b), cols2[2] => x -> sum(x.c)), gd) ==
                        f(d -> (x1=sum(d.b), x2=sum(d.c)), gd)
                    @test f(wrap(cols2[1] => x -> sum(x.b), cols2[2] => x -> first(x.c)), gd) ==
                        f(d -> (x1=sum(d.b), x2=first(d.c)), gd)
                else
                    cols2 = cols
                    @test f(wrap(cols[1] => sum, cols[2] => sum), gd) ==
                        f(d -> (b_sum=sum(d.b), c_sum=sum(d.c)), gd)
                    @test f(wrap(cols[1] => sum, cols[2] => x -> first(x)), gd) ==
                        f(d -> (b_sum=sum(d.b), c_function=first(d.c)), gd)
                    @test_throws ArgumentError f(wrap(cols2[1] => vexp, cols2[2] => sum), gd)
                end

                @test_throws ArgumentError f(wrap(cols => x -> DataFrame(y=exp.(x.b), z=sum(x.c))), gd)
                @test_throws ArgumentError f(wrap(cols => x -> [exp.(x.b) x.c]), gd)
            end
        end
    end
end

struct TestType end
Base.isless(::TestType, ::Int) = true
Base.isless(::Int, ::TestType) = false
Base.isless(::TestType, ::TestType) = false

@testset "combine with aggregation functions" begin
    Random.seed!(1)
    df = DataFrame(a = rand(1:5, 20), x1 = rand(Int, 20), x2 = rand(Complex{Int}, 20))

    for f in (sum, prod, maximum, minimum, mean, var, std, first, last, length)
        gd = groupby_checked(df, :a)

        res = combine(gd, y = :x1 => f)
        expected = combine(gd, y = :x1 => x -> f(x))
        @test res ≅ expected
        @test typeof(res.y) == typeof(expected.y)

        for T in (Union{Missing, Int}, Union{Int, Int8},
                  Union{Missing, Int, Int8})
            df.x3 = Vector{T}(df.x1)
            gd = groupby_checked(df, :a)
            res = combine(gd, y = :x3 => f)
            expected = combine(gd, y = :x3 => x -> f(x))
            @test res ≅ expected
            @test typeof(res.y) == typeof(expected.y)
        end

        f === length && continue

        df.x3 = allowmissing(df.x1)
        df.x3[1] = missing
        gd = groupby_checked(df, :a)
        res = combine(gd, y = :x3 => f)
        expected = combine(gd, y = :x3 => x -> f(x))
        @test res ≅ expected
        @test typeof(res.y) == typeof(expected.y)
        res = combine(gd, y = :x3 => f∘skipmissing)
        expected = combine(gd, y = :x3 => x -> f(collect(skipmissing(x))))
        @test res ≅ expected
        @test typeof(res.y) == typeof(expected.y)
    end
    # Test complex numbers
    for f in (sum, prod, mean, var, std, first, last, length)
        gd = groupby_checked(df, :a)

        res = combine(gd, y = :x2 => f)
        expected = combine(gd, y = :x2 => x -> f(x))
        @test res ≅ expected
        @test typeof(res.y) == typeof(expected.y)
    end
    # Test CategoricalArray
    for f in (maximum, minimum, first, last, length),
        (T, m) in ((Int, false),
                   (Union{Missing, Int}, false), (Union{Missing, Int}, true))
        df.x3 = CategoricalVector{T}(df.x1)
        m && (df.x3[1] = missing)
        gd = groupby_checked(df, :a)
        res = combine(gd, y = :x3 => f)
        expected = combine(gd, y = :x3 => x -> f(x))
        @test res ≅ expected
        @test typeof(res.y) == typeof(expected.y)

        f === length && continue

        res = combine(gd, y = :x3 => f∘skipmissing)
        expected = combine(gd, y = :x3 => x -> f(collect(skipmissing(x))))
        @test res == expected
        @test typeof(res.y) == typeof(expected.y)
        if m
            gd[1].x3 = missing
            @test_throws ArgumentError combine(gd, y = :x3 => f∘skipmissing)
        end
    end
    @test combine(gd, y = :x1 => maximum, z = :x2 => sum) ==
        combine(gd, y = :x1 => x -> maximum(x), z = :x2 => x -> sum(x))
    # Test floating point corner cases
    df = DataFrame(a = [1, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6],
                   x1 = [0.0, 1.0, 2.0, NaN, NaN, NaN, Inf, Inf, Inf, 1.0, NaN, 0.0, -0.0])

    for f in (sum, prod, maximum, minimum, mean, var, std, first, last, length),
        sort in (false, true),
        skip in (false, true)
        gd = groupby_checked(df, :a, sort=sort, skipmissing=skip)

        res = combine(gd, y = :x1 => f)
        expected = combine(gd, y = :x1 => x -> f(x))
        @test res ≅ expected
        @test typeof(res.y) == typeof(expected.y)

        f === length && continue

        df.x3 = allowmissing(df.x1)
        df.x3[1] = missing
        gd = groupby_checked(df, :a, sort=sort, skipmissing=skip)
        res = combine(gd, y = :x3 => f)
        expected = combine(gd, y = :x3 => x -> f(x))
        @test res ≅ expected
        @test typeof(res.y) == typeof(expected.y)
        res = combine(gd, y = :x3 => f∘skipmissing)
        expected = combine(gd, y = :x3 => x -> f(collect(skipmissing(x))))
        @test res ≅ expected
        @test typeof(res.y) == typeof(expected.y)
    end

    df = DataFrame(x = [1, 1, 2, 2], y = Any[1, 2.0, 3.0, 4.0])
    res = by(df, :x, z = :y => maximum)
    @test res.z isa Vector{Float64}
    @test res.z == by(df, :x, z = :y => x -> maximum(x)).z

    # Test maximum when no promotion rule exists
    df = DataFrame(x = [1, 1, 2, 2], y = [1, TestType(), TestType(), TestType()])
    gd = groupby_checked(df, :x)
    for f in (maximum, minimum)
        res = combine(gd, z = :y => maximum)
        @test res.z isa Vector{Any}
        @test res.z == by(df, :x, z = :y => x -> maximum(x)).z
    end
end

@testset "iteration protocol" begin
    gd = groupby_checked(DataFrame(A = [:A, :A, :B, :B], B = 1:4), :A)
    for v in gd
        @test size(v) == (2,2)
    end
end

@testset "getindex" begin
    df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                   b = 1:8)
    gd = groupby_checked(df, :a)
    @test gd[1] isa SubDataFrame
    @test gd[1] == view(df, [1, 5], :)
    @test_throws BoundsError gd[5]
    if VERSION < v"1.0.0-"
        @test gd[true] == gd[1]
    else
        @test_throws ArgumentError gd[true]
    end
    @test_throws MethodError gd["a"]
    gd2 = gd[[true, false, false, false]]
    @test length(gd2) == 1
    @test gd2[1] == gd[1]
    @test_throws BoundsError gd[[true, false]]
    gd3 = gd[:]
    @test gd3 isa GroupedDataFrame
    @test length(gd3) == 4
    @test gd3 == gd
    for i in 1:4
        @test gd3[i] == gd[i]
    end
    gd4 = gd[[1,2]]
    @test gd4 isa GroupedDataFrame
    @test length(gd4) == 2
    for i in 1:2
        @test gd4[i] == gd[i]
    end
    @test_throws BoundsError gd[1:5]
end

@testset "== and isequal" begin
    df1 = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                    b = 1:8)
    df2 = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
                    b = [1:7;missing])
    gd1 = groupby_checked(df1, :a)
    gd2 = groupby_checked(df2, :a)
    @test gd1 == gd1
    @test isequal(gd1, gd1)
    @test ismissing(gd1 == gd2)
    @test !isequal(gd1, gd2)
    @test ismissing(gd2 == gd2)
    @test isequal(gd2, gd2)
    df1.c = df1.a
    df2.c = df2.a
    @test gd1 != groupby_checked(df2, :c)
    df2[7, :b] = 10
    @test gd1 != gd2
    df3 = DataFrame(a = repeat([1, 2, 3, missing], outer=[2]),
                    b = 1:8)
    df4 = DataFrame(a = repeat([1, 2, 3, missing], outer=[2]),
                    b = [1:7;missing])
    gd3 = groupby_checked(df3, :a)
    gd4 = groupby_checked(df4, :a)
    @test ismissing(gd3 == gd4)
    @test !isequal(gd3, gd4)
    gd3 = groupby_checked(df3, :a, skipmissing = true)
    gd4 = groupby_checked(df4, :a, skipmissing = true)
    @test gd3 == gd4
    @test isequal(gd3, gd4)
end

@testset "show" begin
    function capture_stdout(f::Function)
        oldstdout = stdout
        rd, wr = redirect_stdout()
        f()
        str = String(readavailable(rd))
        redirect_stdout(oldstdout)
        size = displaysize(rd)
        close(rd)
        close(wr)
        str, size
    end

    df = DataFrame(A = Int64[1:4;], B = ["x\"", "∀ε>0: x+ε>x", "z\$", "A\nC"],
                   C = Float32[1.0, 2.0, 3.0, 4.0])
    gd = groupby_checked(df, :A)
    io = IOContext(IOBuffer(), :limit=>true)
    show(io, gd)
    str = String(take!(io.io))
    @test str == """
    GroupedDataFrame with 4 groups based on key: A
    First Group (1 row): A = 1
    │ Row │ A     │ B      │ C       │
    │     │ Int64 │ String │ Float32 │
    ├─────┼───────┼────────┼─────────┤
    │ 1   │ 1     │ x"     │ 1.0     │
    ⋮
    Last Group (1 row): A = 4
    │ Row │ A     │ B      │ C       │
    │     │ Int64 │ String │ Float32 │
    ├─────┼───────┼────────┼─────────┤
    │ 1   │ 4     │ A\\nC   │ 4.0     │"""
    show(io, gd, allgroups=true)
    str = String(take!(io.io))
    @test str == """
    GroupedDataFrame with 4 groups based on key: A
    Group 1 (1 row): A = 1
    │ Row │ A     │ B      │ C       │
    │     │ Int64 │ String │ Float32 │
    ├─────┼───────┼────────┼─────────┤
    │ 1   │ 1     │ x\"     │ 1.0     │
    Group 2 (1 row): A = 2
    │ Row │ A     │ B           │ C       │
    │     │ Int64 │ String      │ Float32 │
    ├─────┼───────┼─────────────┼─────────┤
    │ 1   │ 2     │ ∀ε>0: x+ε>x │ 2.0     │
    Group 3 (1 row): A = 3
    │ Row │ A     │ B      │ C       │
    │     │ Int64 │ String │ Float32 │
    ├─────┼───────┼────────┼─────────┤
    │ 1   │ 3     │ z\$     │ 3.0     │
    Group 4 (1 row): A = 4
    │ Row │ A     │ B      │ C       │
    │     │ Int64 │ String │ Float32 │
    ├─────┼───────┼────────┼─────────┤
    │ 1   │ 4     │ A\\nC   │ 4.0     │"""

    # Test two-argument show
    str1, dsize = capture_stdout() do
        show(gd)
    end
    io = IOContext(IOBuffer(), :limit=>true, :displaysize=>dsize)
    show(io, gd)
    str2 = String(take!(io.io))
    @test str1 == str2

    @test sprint(show, "text/html", gd) ==
        "<p><b>GroupedDataFrame with 4 groups based on key: A</b></p>" *
        "<p><i>First Group (1 row): A = 1</i></p><table class=\"data-frame\">" *
        "<thead><tr><th></th><th>A</th><th>B</th><th>C</th></tr><tr><th></th>" *
        "<th>Int64</th><th>String</th><th>Float32</th></tr></thead>" *
        "<tbody><tr><th>1</th><td>1</td><td>x\"</td><td>1.0</td></tr></tbody>" *
        "</table><p>⋮</p><p><i>Last Group (1 row): A = 4</i></p>" *
        "<table class=\"data-frame\"><thead><tr><th></th><th>A</th><th>B</th><th>C</th></tr>" *
        "<tr><th></th><th>Int64</th><th>String</th><th>Float32</th></tr></thead>" *
        "<tbody><tr><th>1</th><td>4</td><td>A\\nC</td><td>4.0</td></tr></tbody></table>"

    @test sprint(show, "text/latex", gd) == """
        GroupedDataFrame with 4 groups based on key: A

        First Group (1 row): A = 1

        \\begin{tabular}{r|ccc}
        \t& A & B & C\\\\
        \t\\hline
        \t& Int64 & String & Float32\\\\
        \t\\hline
        \t1 & 1 & x" & 1.0 \\\\
        \\end{tabular}

        \$\\dots\$

        Last Group (1 row): A = 4

        \\begin{tabular}{r|ccc}
        \t& A & B & C\\\\
        \t\\hline
        \t& Int64 & String & Float32\\\\
        \t\\hline
        \t1 & 4 & A\\textbackslash{}nC & 4.0 \\\\
        \\end{tabular}
        """

    gd = groupby(DataFrame(a=[Symbol("&")], b=["&"]), [1,2])
    @test sprint(show, gd) === """
        GroupedDataFrame with 1 group based on keys: a, b
        Group 1 (1 row): a = :&, b = "&"
        │ Row │ a      │ b      │
        │     │ Symbol │ String │
        ├─────┼────────┼────────┤
        │ 1   │ &      │ &      │"""

    @test sprint(show, "text/html", gd) ==
        "<p><b>GroupedDataFrame with 1 group based on keys: a, b</b></p><p><i>" *
        "First Group (1 row): a = :&, b = \"&\"</i></p>" *
        "<table class=\"data-frame\"><thead><tr><th></th><th>a</th><th>b</th></tr>" *
        "<tr><th></th><th>Symbol</th><th>String</th></tr></thead><tbody><tr><th>1</th>" *
        "<td>&</td><td>&</td></tr></tbody></table>"

    @test sprint(show, "text/latex", gd) == """
        GroupedDataFrame with 1 group based on keys: a, b

        First Group (1 row): a = :\\&, b = "\\&"

        \\begin{tabular}{r|cc}
        \t& a & b\\\\
        \t\\hline
        \t& Symbol & String\\\\
        \t\\hline
        \t1 & \\& & \\& \\\\
        \\end{tabular}
        """
end

end # module


text:
import StatsBase: sample, WeightVec

# test logistic loss

## generate data
srand(1);
m,n,k = 1000,1000,3;
kfit = k+1
# variance of measurement
sigmasq = .1

# coordinates of covariates
X_real = randn(m,k)
# directions of observations
Y_real = randn(k,n)

XY = X_real*Y_real;
A = zeros(Int, (m, n))
logistic(x) = 1/(1+exp(-x))
for i=1:m
	for j=1:n
		A[i,j] = (logistic(XY[i,j]) >= rand()) ? true : false
	end
end

# and the model
losses = LogisticLoss()
rx, ry = QuadReg(.1), QuadReg(.1);
glrm = GLRM(A,losses,rx,ry,kfit)
gg = GGLRM(A, losses, rx, ry, kfit)
fit!(glrm);
fit!(gg);
@test_approx_eq_eps(whole_objective(gg, gg.X'gg.Y), whole_objective(gg, glrm.X'glrm.Y), 0.1*length(A))

ext: jl


text:
# Rohan McLure, Australian National University (2019)
# Distributed implementation of the Transpose Parallel Research Kernel
@everywhere using ArrayChannels

function main()
    # Parameters
    local iterations, order
    if length(ARGS) != 2 || nprocs() == 1
        println("Usage: julia -p <# workers> transpose.jl <# iterations> <# order>")
        exit(1)
    else
        argv = map(x -> parse(Int, x), ARGS)
        iterations, order = argv
    end

    if order % nworkers() != 0
        println("Number of workers must divide order of matrix")
        exit(1)
    end

    # Rank setup
    rank = nworkers()

    # Dimension setup
    std_dim = Int(order / rank)

    # Setup channels
    in_channel = ArrayChannel(Float64, workers(), std_dim, std_dim)
    out_channel = ArrayChannel(Float64, workers(), std_dim, std_dim)
    dims = (order, std_dim)

    futures = Vector{Future}(undef, rank)
    @sync for k in 1 : rank
        futures[k] = @spawnat workers()[k] work(dims, in_channel, out_channel, iterations)
    end

    epsilon = 1.e-8
    results = map(fetch, futures)
    abserr = sum(map(x -> x[1], results))
    time = maximum(map(x -> x[2], results))
    nbytes = 2 * order^2 * 8

    if abserr < epsilon
        println("Solution validates")
    else
        println("error ",abserr, " exceeds threshold ",epsilon)
        println("ERROR: solution did not validate")
    end

    avgtime = time/iterations
    println("Rate (MB/s): ",1.e-6*nbytes/avgtime, " Avg time (s): ", avgtime)
end

# Local transpose operation
@inline @everywhere function do_transpose!(A, B, order)
    for jt in 1 : 32 : order
        for it in 1 : 32 : order
            for j in jt : min(jt+32-1, order)
                for i in it : min(it+32-1, order)
                    @inbounds B[i,j] += A[j,i]
                    @inbounds A[j,i] += 1.0
                end
            end
        end
    end
end

# Same function as do_transpose!, but just reassigns output
@inline @everywhere function write_block!(A, B, order)
    for jt in 1 : 32 : order
        for it in 1 : 32 : order
            for j in jt : min(jt+32-1, order)
                for i in it : min(it+32-1, order)
                    @inbounds B[i,j] = A[j,i]
                    @inbounds A[j,i] += 1.0
                end
            end
        end
    end
end

@everywhere function work(
    dims::Tuple{Int,Int}, # First component will be order
    inbox::ArrayChannel,
    outbox::ArrayChannel,
    iterations::Int
)
    # Global reference to inbox' data
    in_rrid = inbox.rrid

    # Arrays underpinning channels
    in_data = inbox.buffer.src
    out_data = outbox.buffer.src

    # Set up local process data
    id = myid()
    n = nworkers()
    order, width_local = dims

    # Local 'column blocks' of the dist. matrix
    A = zeros(dims...)
    B = zeros(dims...)
    for i in 1:order
        for j in 1:width_local
            jl = j + (id-2)*width_local
            A[i,j] = order * (jl - 1) + (i - 1)
        end
    end

    # Views into the diagonal block
    local_frame_A = view(A, (id-2)*width_local+1:(id-1)*width_local, :)
    local_frame_B = view(B, (id-2)*width_local+1:(id-1)*width_local, :)

    local t0, t1
    for k in 0 : iterations
        if k == 1
            t0 = time_ns()
        end

        do_transpose!(local_frame_A, local_frame_B, width_local)

        for phase in 1 : n-1
            send_idx = (id-2 + phase) % n + 1
            recv_idx = (id-2 + n - phase) % n + 1

            frame = view(A, (send_idx-1)*width_local+1:send_idx*width_local, :)
            write_block!(frame, out_data, width_local)

            @sync begin
                @async begin
                    frame = view(B, (recv_idx-1)*width_local+1:recv_idx*width_local, :)
                    take!(inbox, recv_idx+1)
                    frame .+= in_data
                end
                @async begin
                    put!(outbox, send_idx+1, in_rrid)
                end
            end
        end
    end
    t1 = time_ns()

    # Validate the local solution - sending the error back to master
    addit = (0.5*iterations) * (iterations+1)
    abserr = 0.0
    for i in 1:order
        for j in 1:width_local
            jl = j + (id-2)*width_local
            temp = (order * (i-1) + (jl-1)) * (iterations+1)
            @inbounds abserr = abserr + abs(B[i,j] - (temp+addit))
        end
    end
    return (abserr, (t1-t0)*1.e-9)
end

main()


text:
using PkgBenchmark
results = readresults(joinpath(@__DIR__, "results.json"))
export_markdown(joinpath(@__DIR__, "results.md"), results)


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function xpress_model(MAXTIME, MIPRELSTOP)
    m = direct_model(Xpress.Optimizer(MAXTIME=-MAXTIME, MIPRELSTOP=MIPRELSTOP, OUTPUTLOG = 0)) #logfile="output.log"))
	m[:solver_name] = "Xpress"
	m[:timeout_seconds] = MAXTIME
	m[:optimality_tolerance] = MIPRELSTOP
	return m
end


text:
module MKL

const libmkl = Libdl.find_library(["mkl_rt"])
isempty(libmkl) && throw("MKL library cannot be found.")

include("lib/libmkl_dnn.jl")
include("lib/libmkl_dnn_types.jl")

include("convolution.jl")
include("pooling.jl")
include("relu.jl")

function check_dnnerror(e)
    e == E_SUCCESS && return
    e == E_INCORRECT_INPUT_PARAMETER && throw("INCORRECT_INPUT_PARAMETER.")
    e == E_UNEXPECTED_NULL_POINTER && throw("UNEXPECTED_NULL_POINTER.")
    e == E_MEMORY_ERROR && throw("MEMORY_ERROR.")
    e == E_UNSUPPORTED_DIMENSION && throw("UNSUPPORTED_DIMENSION.")
    e == E_UNIMPLEMENTED && throw("UNIMPLEMENTED.")
end

function dnnLayoutCreate{T}(x::Array{T})
    p = Ptr{Cvoid}[0]
    csize = Csize_t[size(x)...]
    cstrides = Csize_t[strides(x)...]
    dnnLayoutCreate_F32(p, ndims(x), csize, cstrides)
    p[1]
end

function dnnPrimitiveAttributesCreate()
    p = Ptr{Cvoid}[0]
    dnnPrimitiveAttributesCreate_F32(p)
    p[1]
end

function dnnExecute(primitive, resources)
    dnnExecute_F32(primitive, resources)
end

function dnnPrimitiveAttributesDestroy(attributes)
    dnnPrimitiveAttributesDestroy_F32(attributes)
end

function dnnDelete(primitive)
    dnnDelete_F32(primitive)
end

end


text:
module SDDiP

export setSDDiPsolver!,
Pattern,
@binarystate, @binarystates

# Our Lagrangian solver
include(joinpath(dirname(@__FILE__),"Lagrangian","Lagrangian.jl"))

using SDDP, JuMP, Compat, Reexport
@reexport using .Lagrangian

include("solver.jl")
include("binary_expansion.jl")
include("binarystate.jl")
end


text:
# Generate some simple random data just to make sure the
# models work; note that we don't check the validity of
# the models, we assume they work fine, we just check
# they get called properly and have appropriate metadata

Random.seed!(1551)
n, p = 100, 5
X  = randn(n, p)
X1 = hcat(X, ones(100))
θ  = randn(p+1)
y  = X1 * θ .+ 0.1 .* randn(n)

θ_ls = X1 \ y

@testset "ARD" begin
    m = ARDRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.998859, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y),  0.032693, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :alpha, :lambda, :sigma, :scores)
end

@testset "BayesianRidge" begin
    m = BayesianRidgeRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.998859, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.0326918, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :alpha, :lambda, :sigma, :scores)
end

@testset "ElasticNet" begin
    m = ElasticNetRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.769795, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.464447, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
end

@testset "ElasticNetCV" begin
    m = ElasticNetCVRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.99884, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.0328523, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :l1_ratio, :mse_path, :alphas)
end

@testset "Huber" begin
    m = HuberRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.99884, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.032935, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :scale, :outliers)
end

@testset "Lars" begin
    m = LarsRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.998859, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.0326918, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :alphas, :active, :coef_path)
end

@testset "LarsCV" begin
    m = LarsCVRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.998859, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.0326918, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :alpha, :alphas, :cv_alphas, :mse_path, :coef_path)
end

@testset "Lasso" begin
    m = LassoRegressor(alpha = 0.1)
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.9946343, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.0709070, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
end

@testset "LassoCV" begin
    m = LassoCVRegressor(random_state=0, cv=6)
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.998857, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.03272594, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :alpha, :alphas, :mse_path, :dual_gap)
end

@testset "LassoLars" begin
    m = LassoLarsRegressor(alpha = 0.01)
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.994317, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.072973, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :alphas, :active, :coef_path)
end

@testset "LassoLarsCV" begin
    m = LassoLarsCVRegressor(cv=4)
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.998859, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.0326918, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :coef_path, :alpha, :alphas, :cv_alphas, :mse_path)
end

@testset "LassoLarsIC" begin
    m = LassoLarsICRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.9920703, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.08619978, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :alpha)
end

@testset "LinReg" begin
    m = LinearRegressor()
    f, _, _ = fit(m, 1, X, y)
    # test if things work
    @test isapprox(f[1].score(X, y), 0.998859, rtol=1e-4)
    @test isapprox(norm(predict(m, f, X) .- y) / norm(y), 0.032691, rtol=1e-4)
    fp = fitted_params(m, f)
    @test fp.coef      ≈ θ_ls[1:p]
    @test fp.intercept ≈ θ_ls[end]
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
end

@testset "OMP" begin
    m = OrthogonalMatchingPursuitRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.4593868, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.711741, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
end

@testset "OMPCV" begin
    m = OrthogonalMatchingPursuitCVRegressor(cv = 5)
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.998859, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.0326918, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :n_nonzero_coefs)
end

@testset "OMPCV" begin
    m = OrthogonalMatchingPursuitCVRegressor(cv = 5)
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.998859, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.0326918, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :n_nonzero_coefs)
end

@testset "PassAggr" begin
    m = PassiveAggressiveRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test 0.980 ≤ f[1].score(X, y) ≤ 0.999
    @test norm(predict(m, f, X) .- y)/norm(y) ≤ 0.1
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
end

@testset "RANSAC" begin
    m = RANSACRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test 0.99 ≤ f[1].score(X, y) ≤ 0.999
    @test 0.02 ≤ norm(predict(m, f, X) .- y)/norm(y) ≤ 0.05
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:estimator, :n_trials, :inlier_mask, :n_skips_no_inliers, :n_skips_invalid_data, :n_skips_invalid_model)
end

@testset "Ridge" begin
    m = RidgeRegressor(alpha = 1.0)
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.998793, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.0336254, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
end

@testset "RidgeCV" begin
    m = RidgeCVRegressor(cv=nothing, store_cv_values=true)
    f, _, _ = fit(m, 1, X, y)
    @test isapprox(f[1].score(X, y), 0.998858, rtol=1e-5)
    @test isapprox(norm(predict(m, f, X) .- y)/norm(y), 0.0327014, rtol=1e-5)
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :alpha, :cv_values)
    @test fp[:cv_values] isa Matrix
end

@testset "SGDReg" begin
    m = SGDRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test 0.980 ≤ f[1].score(X, y) ≤ 0.999
    @test norm(predict(m, f, X) .- y)/norm(y) ≤ 0.1
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :average_coef, :average_intercept)
end

@testset "TheilSen" begin
    m = TheilSenRegressor()
    f, _, _ = fit(m, 1, X, y)
    @test 0.980 ≤ f[1].score(X, y) ≤ 0.999
    @test norm(predict(m, f, X) .- y)/norm(y) ≤ 0.1
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :breakdown, :n_subpopulation)
end

##############
# MULTI TASK #
##############

y2 = (t1=y, t2=y)

@testset "MTLassoCV" begin
    m = MultiTaskLassoRegressor()
    f, _, _ = fit(m, 1, X, y2)
    @test f[1].coef_ isa Matrix
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
    pred = predict(m, f, X)
    @test pred isa Tables.MatrixTable
    @test Tables.schema(pred).names == (:t1, :t2)
end

@testset "MTLassoCV" begin
    m = MultiTaskLassoCVRegressor(cv = 5, random_state = 0)
    f, _, _ = fit(m, 1, X, y2)
    @test f[1].coef_ isa Matrix
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :alpha, :mse_path, :alphas)
    @test predict(m, f, X) isa Tables.MatrixTable
end

X = MLJBase.table([0 0; 1 1; 2 2])
y = MLJBase.table([0 0; 1 1; 2 2])

@testset "MTElNet" begin
    m = MultiTaskElasticNetRegressor(alpha = 0.1, l1_ratio = 0.5)
    f, _, _ = fit(m, 1, X, y)
    @test f[1].coef_ isa Matrix
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept)
    @test predict(m, f, X) isa Tables.MatrixTable
end

@testset "MTElNetCV" begin
    m = MultiTaskElasticNetCVRegressor(cv = 3, n_alphas = 100, l1_ratio = 0.5)
    f, _, _ = fit(m, 1, X, y)
    @test f[1].coef_ isa Matrix
    # testing that the fitted params is proper
    fp = fitted_params(m, f)
    @test keys(fp) == (:coef, :intercept, :alpha, :mse_path, :l1_ratio)
    @test predict(m, f, X) isa Tables.MatrixTable
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Read

import ...LibGit2, ..Cache, ..Reqs, ...Pkg.PkgError, ..Dir
using ..Types

readstrip(path...) = strip(readstring(joinpath(path...)))

url(pkg::AbstractString) = readstrip(Dir.path("METADATA"), pkg, "url")
sha1(pkg::AbstractString, ver::VersionNumber) =
    readstrip(Dir.path("METADATA"), pkg, "versions", string(ver), "sha1")

function available(names=readdir("METADATA"))
    pkgs = Dict{String,Dict{VersionNumber,Available}}()
    for pkg in names
        isfile("METADATA", pkg, "url") || continue
        versdir = joinpath("METADATA", pkg, "versions")
        isdir(versdir) || continue
        for ver in readdir(versdir)
            ismatch(Base.VERSION_REGEX, ver) || continue
            isfile(versdir, ver, "sha1") || continue
            haskey(pkgs,pkg) || (pkgs[pkg] = Dict{VersionNumber,Available}())
            pkgs[pkg][convert(VersionNumber,ver)] = Available(
                readchomp(joinpath(versdir,ver,"sha1")),
                Reqs.parse(joinpath(versdir,ver,"requires"))
            )
        end
    end
    return pkgs
end
available(pkg::AbstractString) = get(available([pkg]),pkg,Dict{VersionNumber,Available}())

function latest(names=readdir("METADATA"))
    pkgs = Dict{String,Available}()
    for pkg in names
        isfile("METADATA", pkg, "url") || continue
        versdir = joinpath("METADATA", pkg, "versions")
        isdir(versdir) || continue
        pkgversions = VersionNumber[]
        for ver in readdir(versdir)
            ismatch(Base.VERSION_REGEX, ver) || continue
            isfile(versdir, ver, "sha1") || continue
            push!(pkgversions, convert(VersionNumber,ver))
        end
        isempty(pkgversions) && continue
        ver = string(maximum(pkgversions))
        pkgs[pkg] = Available(
                readchomp(joinpath(versdir,ver,"sha1")),
                Reqs.parse(joinpath(versdir,ver,"requires"))
            )
    end
    return pkgs
end

isinstalled(pkg::AbstractString) =
    pkg != "METADATA" && pkg != "REQUIRE" && pkg[1] != '.' && isdir(pkg)

function isfixed(pkg::AbstractString, prepo::LibGit2.GitRepo, avail::Dict=available(pkg))
    isinstalled(pkg) || throw(PkgError("$pkg is not an installed package."))
    isfile("METADATA", pkg, "url") || return true
    ispath(pkg, ".git") || return true

    LibGit2.isdirty(prepo) && return true
    LibGit2.isattached(prepo) && return true
    LibGit2.need_update(prepo)
    if isnull(find("REQUIRE", LibGit2.GitIndex(prepo)))
        isfile(pkg,"REQUIRE") && return true
    end
    head = string(LibGit2.head_oid(prepo))
    for (ver,info) in avail
        head == info.sha1 && return false
    end

    cache = Cache.path(pkg)
    cache_has_head = if isdir(cache)
        crepo = LibGit2.GitRepo(cache)
        LibGit2.iscommit(head, crepo)
    else
        false
    end
    res = true
    try
        for (ver,info) in avail
            if cache_has_head && LibGit2.iscommit(info.sha1, crepo)
                if LibGit2.is_ancestor_of(head, info.sha1, crepo)
                    res = false
                    break
                end
            elseif LibGit2.iscommit(info.sha1, prepo)
                if LibGit2.is_ancestor_of(head, info.sha1, prepo)
                    res = false
                    break
                end
            else
                Base.warn_once("unknown $pkg commit $(info.sha1[1:8]), metadata may be ahead of package cache")
            end
        end
    finally
        cache_has_head && LibGit2.close(crepo)
    end
    return res
end

function ispinned(pkg::AbstractString)
    ispath(pkg,".git") || return false
    LibGit2.with(LibGit2.GitRepo, pkg) do repo
        return ispinned(repo)
    end
end

function ispinned(prepo::LibGit2.GitRepo)
    LibGit2.isattached(prepo) || return false
    br = LibGit2.branch(prepo)
    # note: regex is based on the naming scheme used in Entry.pin()
    return ismatch(r"^pinned\.[0-9a-f]{8}\.tmp$", br)
end

function installed_version(pkg::AbstractString, prepo::LibGit2.GitRepo, avail::Dict=available(pkg))
    ispath(pkg,".git") || return typemin(VersionNumber)

    # get package repo head hash
    local head
    try
        head = string(LibGit2.head_oid(prepo))
    catch ex
        # refs/heads/master does not exist
        if isa(ex,LibGit2.GitError) &&
            ex.code == LibGit2.Error.EUNBORNBRANCH
            head = ""
        else
            rethrow(ex)
        end
    end
    isempty(head) && return typemin(VersionNumber)

    vers = collect(keys(filter((ver,info)->info.sha1==head, avail)))
    !isempty(vers) && return maximum(vers)

    cache = Cache.path(pkg)
    cache_has_head = if isdir(cache)
        crepo = LibGit2.GitRepo(cache)
        LibGit2.iscommit(head, crepo)
    else
        false
    end
    ancestors = VersionNumber[]
    descendants = VersionNumber[]
    try
        for (ver,info) in avail
            sha1 = info.sha1
            base = if cache_has_head && LibGit2.iscommit(sha1, crepo)
                LibGit2.merge_base(crepo, head, sha1)
            elseif LibGit2.iscommit(sha1, prepo)
                LibGit2.merge_base(prepo, head, sha1)
            else
                Base.warn_once("unknown $pkg commit $(sha1[1:8]), metadata may be ahead of package cache")
                continue
            end
            string(base) == sha1 && push!(ancestors,ver)
            string(base) == head && push!(descendants,ver)
        end
    finally
        cache_has_head && LibGit2.close(crepo)
    end
    both = sort!(intersect(ancestors,descendants))
    isempty(both) || warn("$pkg: some versions are both ancestors and descendants of head: $both")
    if !isempty(descendants)
        v = minimum(descendants)
        return VersionNumber(v.major, v.minor, v.patch, ("",), ())
    elseif !isempty(ancestors)
        v = maximum(ancestors)
        return VersionNumber(v.major, v.minor, v.patch, (), ("",))
    else
        return typemin(VersionNumber)
    end
end

function requires_path(pkg::AbstractString, avail::Dict=available(pkg))
    pkgreq = joinpath(pkg,"REQUIRE")
    ispath(pkg,".git") || return pkgreq
    repo = LibGit2.GitRepo(pkg)
    head = LibGit2.with(LibGit2.GitRepo, pkg) do repo
        LibGit2.isdirty(repo, "REQUIRE") && return pkgreq
        LibGit2.need_update(repo)
        if isnull(find("REQUIRE", LibGit2.GitIndex(repo)))
            isfile(pkgreq) && return pkgreq
        end
        string(LibGit2.head_oid(repo))
    end
    for (ver,info) in avail
        if head == info.sha1
            return joinpath("METADATA", pkg, "versions", string(ver), "requires")
        end
    end
    return pkgreq
end

requires_list(pkg::AbstractString, avail::Dict=available(pkg)) =
    collect(keys(Reqs.parse(requires_path(pkg,avail))))

requires_dict(pkg::AbstractString, avail::Dict=available(pkg)) =
    Reqs.parse(requires_path(pkg,avail))

function installed(avail::Dict=available())
    pkgs = Dict{String,Tuple{VersionNumber,Bool}}()
    for pkg in readdir()
        isinstalled(pkg) || continue
        ap = get(avail,pkg,Dict{VersionNumber,Available}())
        if ispath(pkg,".git")
            LibGit2.with(LibGit2.GitRepo, pkg) do repo
                ver = installed_version(pkg, repo, ap)
                fixed = isfixed(pkg, repo, ap)
                pkgs[pkg] = (ver, fixed)
            end
        else
            pkgs[pkg] = (typemin(VersionNumber), true)
        end
    end
    return pkgs
end

function fixed(avail::Dict=available(), inst::Dict=installed(avail), dont_update::Set{String}=Set{String}(),
    julia_version::VersionNumber=VERSION)
    pkgs = Dict{String,Fixed}()
    for (pkg,(ver,fix)) in inst
        (fix || pkg in dont_update) || continue
        ap = get(avail,pkg,Dict{VersionNumber,Available}())
        pkgs[pkg] = Fixed(ver,requires_dict(pkg,ap))
    end
    pkgs["julia"] = Fixed(julia_version)
    return pkgs
end

function free(inst::Dict=installed(), dont_update::Set{String}=Set{String}())
    pkgs = Dict{String,VersionNumber}()
    for (pkg,(ver,fix)) in inst
        (fix || pkg in dont_update) && continue
        pkgs[pkg] = ver
    end
    return pkgs
end

function issue_url(pkg::AbstractString)
    ispath(pkg,".git") || return ""
    m = match(LibGit2.GITHUB_REGEX, url(pkg))
    m === nothing && return ""
    return "https://github.com/" * m.captures[1] * "/issues"
end

end # module

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
module CPLEX

const _DEPS_FILE = joinpath(dirname(@__DIR__), "deps", "deps.jl")
if isfile(_DEPS_FILE)
    include(_DEPS_FILE)
else
    error("CPLEX not properly installed. Please run Pkg.build(\"CPLEX\")")
end

using CEnum
import SparseArrays

function _get_version_number()
    if libcplex == "julia_registryci_automerge"
        return VersionNumber(12, 10, 0)  # Fake a valid version for AutoMerge.
    end
    stat = Ref{Cint}()
    env = ccall((:CPXopenCPLEX, libcplex), Ptr{Cvoid}, (Ptr{Cint},), stat)
    p = ccall((:CPXversion, libcplex), Cstring, (Ptr{Cvoid},), env)
    version_string = unsafe_string(p)
    ccall((:CPXcloseCPLEX, libcplex), Cint, (Ptr{Ptr{Cvoid}},), Ref(env))
    return VersionNumber(parse.(Int, split(version_string, ".")[1:3])...)
end

const _CPLEX_VERSION = _get_version_number()

if _CPLEX_VERSION == v"12.10.0"
    include("gen1210/ctypes.jl")
    include("gen1210/libcpx_common.jl")
    include("gen1210/libcpx_api.jl")
elseif _CPLEX_VERSION == v"20.1.0"
    include("gen2010/ctypes.jl")
    include("gen2010/libcpx_common.jl")
    include("gen2010/libcpx_api.jl")
else
    error("""
You have installed version $_CPLEX_VERSION of CPLEX, which is not supported
by CPLEX.jl. We require CPLEX version 12.10 or 20.1.

After installing CPLEX, run:

    import Pkg
    Pkg.rm("CPLEX")
    Pkg.add("CPLEX")

Make sure you set the environment variable `CPLEX_STUDIO_BINARIES` following
the instructions in the CPLEX.jl README, which is available at
https://github.com/jump-dev/CPLEX.jl.

If you have a newer version of CPLEX installed, changes may need to be made
to the Julia code. Please open an issue at
https://github.com/jump-dev/CPLEX.jl.""")
end

include("MOI/MOI_wrapper.jl")
include("MOI/MOI_callbacks.jl")
include("MOI/conflicts.jl")
include("MOI/indicator_constraint.jl")

# CPLEX exports all `CPXxxx` symbols. If you don't want all of these symbols in
# your environment, then use `import CPLEX` instead of `using CPLEX`.

for sym in names(@__MODULE__, all = true)
    sym_string = string(sym)
    if startswith(sym_string, "CPX")
        @eval export $sym
    end
end

include("deprecated_functions.jl")

# Special overload to deprecate the `model.inner` field access.
function Base.getproperty(opt::Optimizer, key::Symbol)
    if key == :inner
        error(_DEPRECATED_ERROR_MESSAGE)
    end
    return getfield(opt, key)
end

end


text:
function rotation_matrix(q0::Float64, q1::Float64, q2::Float64, q3::Float64)

	a11::Float64 = q0^2 + q1^2 - q2^2 - q3^2
	a12::Float64 = 2 * q1 * q2 + 2 * q0 * q3
	a13::Float64 = 2 * q1 * q3 - 2 * q0 * q2
	a21::Float64 = 2 * q1 * q2 - 2 * q0 * q3
	a22::Float64 = q0^2 - q1^2 + q2^2 - q3^2
	a23::Float64 = 2 * q2 * q3 + 2 * q0 * q1
	a31::Float64 = 2 * q1 * q3 + 2 * q0 * q2
	a32::Float64 = 2 * q2 * q3 - 2 * q0 * q1
	a33::Float64 = q0^2 - q1^2 - q2^2 + q3^2

	return (a11, a12, a13, a21, a22, a23, a31, a32, a33)

end


text:
using LightGraphs
include("util.jl")

struct Point
    x::Real
    y::Real
    z::Real
    t::Real
end

function mdist(x::Point,y::Point)::Integer
    return abs(x.x-y.x)+abs(x.y-y.y)+abs(x.z-y.z)+abs(x.t-y.t)
end

function parse_points()::Dict{Integer,Point}
    retval::Dict{Integer,Point} = Dict()
    pointnum = 1
    for line in getlines("inputs/25.txt")
        parts = map(x->parse(Int64,x), split(line, ","))
        retval[pointnum] = Point(parts...)
        pointnum += 1
    end
    return retval
end

function make_graph(points::Dict{Integer,Point})::SimpleGraph
    g = SimpleGraph(length(points))
    for i in 1:length(points)
        p = points[i]
        for j in 1:length(points)
            q = points[j]
            if mdist(p,q) <= 3
                add_edge!(g, i, j)
            end
        end
    end
    return g
end

points = parse_points()
g = make_graph(points)
part_1 = length(connected_components(g))


text:
module LibInt2Module
using Libdl

  """
  The maximal LQuantumNumber to which libint2 has been compiled. Compare :/deps/usr/src/Makefile: MAX_AM.
  Creating shells with higher LQuantumNumber causes segfaults within the libint2.
  """
  const MAX_AM = 7

export libInt2Initialize, libInt2Finalize
export LibInt2Shell
export LibInt2Engine, LibInt2EngineCoulomb, LibInt2EngineOverlap, LibInt2EngineKinetic, LibInt2EngineNuclearAttraction
export destroy!, lqn, nprims
export computeBasisShellsLibInt2
import ..IntegralsModule.computeMatrixBlockOverlap
import ..IntegralsModule.computeMatrixBlockKinetic
import ..IntegralsModule.computeMatrixBlockNuclearAttraction
import ..IntegralsModule.computeTensorBlockElectronRepulsionIntegrals
import ..ShellModule.computeFactorsNonAxial

import ...QuantumLab.libint2_available

if (libint2_available) # the normal case

  import ...QuantumLab.lib_path
  libint2jl = Ptr{Nothing}(0)
  function __init__()
    push!(Libdl.DL_LOAD_PATH,lib_path)
    try
      Libdl.dlopen("libint2-QuantumLab.so")
      global libint2jl = Libdl.dlopen("libint2jl.so")
      LibInt2Module.libInt2Initialize()
    catch e
      @warn "QuantumLab was compiled with libint2_QuantumLab.so. Maybe it has disappeared?"
      error(e)
    end
  end

  import Base.convert
  import Base.show
  import ..ShellModule.nbf
  using TensorOperations
  using ..BaseModule
  using ..GeometryModule
  using ..BasisSetModule
  using ..ShellModule
  using ..BasisFunctionsModule
  using ..IntegralsModule
  using ..AtomModule

  ## Bitstypes

  ## LibInt2Shell
  #  Type Declaration
  primitive type LibInt2Shell64 <: AbstractShell 64 end
  primitive type LibInt2Shell32 <: AbstractShell 32 end
  if Int === Int32
    const LibInt2Shell = LibInt2Shell32
  else
    const LibInt2Shell = LibInt2Shell64
  end

  #  Constructors
  """
      LibInt2Shell(origin::Vector{Float64}, lqn::Int, nprim::Int, exponents::Vector{Float64}, coefficients::Vector{Float64}; renorm::Bool=false)
  The libint2 library expects the coefficients to be input the same way as with basis set definition files. It then renormalizes the coefficients accordingly.
  This behavior is the default for efficiency reasons. To directly enter the coefficients set the renorm flag to false.
  When converting a Shell into a LibInt2Shell (e.g. with convert()), this is taken care of automatically.
  Note, that coefficients are generally only specified up to a global scaling factor - only relative factors are handled by renormalization.
  """
  function LibInt2Shell(origin::Vector{Float64},lqn::Int,nprim::Int,exponents::Vector{Float64},coefficients::Vector{Float64}; renorm::Bool=true)
    if lqn > MAX_AM
      error("libint2 has only been compiled to MAX_AM $MAX_AM.
      Can't create LibInt2Shell with lqn larger than that ($lqn).
      Consider recompiling to higher MAX_AM, or use `Shell`s instead of `LibInt2Shell`s.")
    end
    if renorm==true
      return reinterpret(LibInt2Shell,ccall(dlsym(libint2jl,:_Z11createShellPdiiS_S_),Ptr{Nothing},(Ptr{Float64},Int,Int,Ptr{Float64},Ptr{Float64}),origin,lqn,nprim,exponents,coefficients))
    else
      return convert(LibInt2Shell,Shell(Position(origin...),LQuantumNumber(lqn),exponents,coefficients;renorm=false))
    end
  end

  function LibInt2Shell(origin::Position,lqn::LQuantumNumber,exponents::Vector{Float64},coefficients::Vector{Float64}; renorm::Bool=true)
    orig = convert(Vector{Float64},origin)
    nprim = length(coefficients)
    LibInt2Shell(orig,lqn.exponent,nprim,exponents,coefficients;renorm=renorm)
  end

  function convert(::Type{LibInt2Shell},sh::Shell)
    scaledcoefficients = Float64[]
    for (coeff,expon) in zip(sh.coefficients,sh.exponents)
      pgb = PrimitiveGaussianBasisFunction(origin,expon,MQuantumNumber(sh.lqn.exponent,0,0))
      push!(scaledcoefficients,coeff * sqrt(computeIntegralOverlap(pgb,pgb)))
    end

    LibInt2Shell([sh.center.x,sh.center.y,sh.center.z],sh.lqn.exponent,length(sh.exponents),sh.exponents,scaledcoefficients)
  end

  #  Destructor
  function destroy!(shell::LibInt2Shell)
    ccall(dlsym(libint2jl,:_Z12destroyShellPN7libint25ShellE),Nothing,(LibInt2Shell,),shell)
  end

  #  Further Functions
  function show(io::IO,::MIME"text/plain",sh::LibInt2Shell)
    ccall(dlsym(libint2jl,:_Z10printShellPN7libint25ShellE),Nothing,(LibInt2Shell,),sh)
  end

  function lqn(l2sh::LibInt2Shell)
    # cmp. convert(::Type{Shell},::LibInt2Shell)
    sh_ptrptrptr = reinterpret(Ptr{Ptr{Ptr{Cdouble}}},l2sh)
    contract_ptrptr = unsafe_load(sh_ptrptrptr,4)
    contract_ptr = reinterpret(Ptr{Cint},contract_ptrptr)
    return Int(unsafe_load(contract_ptr,1))
  end

  function nprims(l2sh::LibInt2Shell)
    # cmp. convert(::Type{Shell},::LibInt2Shell)
    sh_ptrptr    = reinterpret(Ptr{Ptr{Cdouble}},l2sh)
    alpha_ptr = unsafe_load(sh_ptrptr,1)
    alphaEnd_ptr = unsafe_load(sh_ptrptr,2)
    return Int(div(alphaEnd_ptr-alpha_ptr,8))
  end

  nbf(shell::LibInt2Shell) = BaseModule.numberMQNsCartesian(LQuantumNumber(lqn(shell)))

  function computeDimensions(shells::Vector{LibInt2Shell})
    totaldim = 0
    maxprims = 0
    maxlqn = 0
    for sh in shells
      totaldim += nbf(sh)
      maxprims = maxprims < nprims(sh) ? nprims(sh) : maxprims
      maxlqn = maxlqn < lqn(sh) ? lqn(sh) : maxlqn
    end
    return (totaldim, maxprims, LQuantumNumber(maxlqn))
  end

  function convert(::Type{Shell},l2sh::LibInt2Shell)
    # LibInt2Shell objects reside in memory as
    # class Shell {
    #   std::vector<real_t> alpha;
    #   std::vector<Contraction> contr;
    #   std::array<real_t, 3> O;        //!< this is inlined in memory
    #   std::vector<real_t> max_ln_coeff;
    # }
    # where
    # struct Contraction {
    #   int l;
    #   bool pure;
    #   std::vector<real_t> coeff;
    # }
    # A std::vector is of the same size as 3 ptrs = 24 bytes on 64bit with gnu compiler.
    # The first pointer gives the start of the allocated space, the second pointer gives the pointer
    # beyond the last vector element and the third pointer beyond the allocated space of the vector
    # (cmp. std::vector.reserve()).

    # I assume one contraction per shell -> no shared-sp

    sh_ptrptrptr = reinterpret(Ptr{Ptr{Ptr{Cdouble}}},l2sh)
    sh_ptrptr    = reinterpret(Ptr{Ptr{Cdouble}},l2sh)
    sh_ptr       = reinterpret(Ptr{Cdouble},l2sh)

    alpha_ptr = unsafe_load(sh_ptrptr,1)
    alphaEnd_ptr = unsafe_load(sh_ptrptr,2)
    primitivecount = div(alphaEnd_ptr-alpha_ptr,8)
    expons = unsafe_wrap(Array,alpha_ptr,primitivecount)

    center_ptr = sh_ptr + 24*2
    center = Position(unsafe_wrap(Array,center_ptr,3)...)

    contract_ptrptr = unsafe_load(sh_ptrptrptr,4)
    coeff_ptr = unsafe_load(contract_ptrptr,2) # the int and bool together take the space of one pointer
    coeffs = unsafe_wrap(Array,coeff_ptr,primitivecount)

    contract_ptr = reinterpret(Ptr{Cint},contract_ptrptr)
    lqn = LQuantumNumber(Int(unsafe_load(contract_ptr,1)))

    return Shell(center,lqn,expons,coeffs;renorm=false)
  end

  ## LibInt2Engine
  #  Type Declaration
  primitive type LibInt2Engine64 64 end
  primitive type LibInt2Engine32 32 end
  if Int === Int32
    const LibInt2Engine = LibInt2Engine32
  else
    const LibInt2Engine = LibInt2Engine64
  end

  #  Constructors
  function LibInt2EngineCoulomb(maxNumberPrimitives::Int,maxLQN::LQuantumNumber)
    reinterpret(LibInt2Engine,ccall(dlsym(libint2jl,:_Z19createEngineCoulombii),Ptr{Nothing},(Cint,Cint),maxNumberPrimitives,maxLQN.exponent))
  end
  function LibInt2EngineKinetic(maxNumberPrimitives::Int,maxLQN::LQuantumNumber)
    reinterpret(LibInt2Engine,ccall(dlsym(libint2jl,:_Z19createEngineKineticii),Ptr{Nothing},(Cint,Cint),maxNumberPrimitives,maxLQN.exponent))
  end
  function LibInt2EngineNuclearAttraction(maxNumberPrimitives::Int,maxLQN::LQuantumNumber,atom::Atom)
    l2atoms = [(atom.element.atomicNumber,atom.position.x,atom.position.y,atom.position.z)]
    reinterpret(LibInt2Engine,ccall(dlsym(libint2jl,:_Z29createEngineNuclearAttractioniiPN7libint24AtomEi),Ptr{Nothing},(Cint,Cint,Ptr{Nothing},Cint),maxNumberPrimitives,maxLQN.exponent,l2atoms,length(l2atoms)))
  end
  function LibInt2EngineNuclearAttraction(maxNumberPrimitives::Int,maxLQN::LQuantumNumber,geo::Geometry)
    l2atoms = Vector{Tuple{Int,Float64,Float64,Float64}}()
    for atom in geo.atoms
      push!(l2atoms,Tuple{Int,Float64,Float64,Float64}((atom.element.atomicNumber,atom.position.x,atom.position.y,atom.position.z)))
    end
    reinterpret(LibInt2Engine,ccall(dlsym(libint2jl,:_Z29createEngineNuclearAttractioniiPN7libint24AtomEi),Ptr{Nothing},(Cint,Cint,Ptr{Nothing},Cint),maxNumberPrimitives,maxLQN.exponent,l2atoms,length(l2atoms)))
  end
  function LibInt2EngineOverlap(maxNumberPrimitives::Int,maxLQN::LQuantumNumber)
    reinterpret(LibInt2Engine,ccall(dlsym(libint2jl,:_Z19createEngineOverlapii),Ptr{Nothing},(Cint,Cint),maxNumberPrimitives,maxLQN.exponent))
  end

  #  Destructor
  function destroy!(engine::LibInt2Engine)
    ccall(dlsym(libint2jl,:_Z13destroyEnginePN7libint26EngineE),Nothing,(LibInt2Engine,),engine)
  end

  ## library functionality
  function libInt2Initialize()
    ccall(dlsym(libint2jl,:_Z12libint2startv),Nothing,())
  end

  function libInt2Finalize()
    ccall(dlsym(libint2jl,:_Z11libint2stopv),Nothing,())
  end

  function computeFactorsNonAxial(shells::Vector{LibInt2Shell})
    naxfacs = Vector{Float64}()
    for sh in shells
      push!(naxfacs,computeFactorsNonAxial(LQuantumNumber(lqn(sh)))...)
    end
    return naxfacs
  end

  """
  If a LibInt2Engine is specified, it is used without assertion if it is of the correct type.
  """
  function IntegralsModule.computeMatrixBlockOverlap(engine::LibInt2Engine, μ::LibInt2Shell, ν::LibInt2Shell)
    buf = ccall(dlsym(libint2jl,:_Z13compute2cIntsPN7libint26EngineEPNS_5ShellES3_),Ptr{Cdouble},(LibInt2Engine,LibInt2Shell,LibInt2Shell), engine, μ,ν)
    block = reshape(unsafe_wrap(Array,buf,nbf(μ)*nbf(ν)),(nbf(ν),nbf(μ)))' # the inner loop within LibInt2 is over the second shell
    ShellModule.renormNonAxial!(block,[μ],[ν])
    return block
  end

  function IntegralsModule.computeMatrixBlockOverlap(μlib::LibInt2Shell,νlib::LibInt2Shell)
    totaldim, maxprims, maxlqn = computeDimensions([μlib,νlib])
    engine = LibInt2EngineOverlap(maxprims,maxlqn)

    result = copy(computeMatrixBlockOverlap(engine,μlib,νlib))

    destroy!(engine)
    return result
  end

  """
  If a LibInt2Engine is specified, it is used without assertion if it is of the correct type.
  """
  function computeMatrixBlockKinetic(engine::LibInt2Engine, μ::LibInt2Shell, ν::LibInt2Shell)
    buf = ccall(dlsym(libint2jl,:_Z13compute2cIntsPN7libint26EngineEPNS_5ShellES3_),Ptr{Cdouble},(LibInt2Engine,LibInt2Shell,LibInt2Shell), engine, μ,ν)
    block = reshape(unsafe_wrap(Array,buf,nbf(μ)*nbf(ν)),(nbf(ν),nbf(μ)))' # the inner loop within LibInt2 is over the second shell
    ShellModule.renormNonAxial!(block,[μ],[ν])
    return block
  end

  function computeMatrixBlockKinetic(μlib::LibInt2Shell,νlib::LibInt2Shell)
    totaldim, maxprims, maxlqn = computeDimensions([μlib,νlib])
    engine = LibInt2EngineKinetic(maxprims,maxlqn)

    result = copy(computeMatrixBlockKinetic(engine,μlib,νlib))

    destroy!(engine)
    return result
  end

  """
  If a LibInt2Engine is specified, it is used without assertion if it is of the correct type.
  """
  function computeMatrixBlockNuclearAttraction(engine::LibInt2Engine, μ::LibInt2Shell, ν::LibInt2Shell)
    buf = ccall(dlsym(libint2jl,:_Z13compute2cIntsPN7libint26EngineEPNS_5ShellES3_),Ptr{Cdouble},(LibInt2Engine,LibInt2Shell,LibInt2Shell), engine, μ,ν)
    block = reshape(unsafe_wrap(Array,buf,nbf(μ)*nbf(ν)),(nbf(ν),nbf(μ)))' # the inner loop within LibInt2 is over the second shell
    ShellModule.renormNonAxial!(block,[μ],[ν])
    return block
  end

  function computeMatrixBlockNuclearAttraction(μlib::LibInt2Shell,νlib::LibInt2Shell, atom::Atom)
    totaldim, maxprims, maxlqn = computeDimensions([μlib,νlib])
    engine = LibInt2EngineNuclearAttraction(maxprims,maxlqn,atom)

    result = copy(computeMatrixBlockNuclearAttraction(engine,μlib,νlib))

    destroy!(engine)
    return result
  end

  function computeMatrixBlockNuclearAttraction(μlib::LibInt2Shell,νlib::LibInt2Shell, geo::Geometry)
    totaldim, maxprims, maxlqn = computeDimensions([μlib,νlib])
    engine = LibInt2EngineNuclearAttraction(maxprims,maxlqn,geo)

    result = copy(computeMatrixBlockNuclearAttraction(engine,μlib,νlib))

    destroy!(engine)
    return result
  end

  """
  If a LibInt2Engine is specified, it is used without assertion if it is of the correct type.
  """
  function computeTensorBlockElectronRepulsionIntegrals(engine::LibInt2Engine, μ::LibInt2Shell, ν::LibInt2Shell, λ::LibInt2Shell, σ::LibInt2Shell)
    buf = ccall(dlsym(libint2jl,:_Z13compute4cIntsPN7libint26EngineEPNS_5ShellES3_S3_S3_),Ptr{Cdouble},(LibInt2Engine,LibInt2Shell,LibInt2Shell,LibInt2Shell,LibInt2Shell), engine, σ,λ,ν,μ)
    block = reshape(unsafe_wrap(Array,buf,nbf(μ)*nbf(ν)*nbf(λ)*nbf(σ)),(nbf(μ),nbf(ν),nbf(λ),nbf(σ))) # Apparently for ERIs the order is correct already (in contrast to the matrix functions
    ShellModule.renormNonAxial!(block,[μ],[ν],[λ],[σ])
    return block
  end

  function computeTensorBlockElectronRepulsionIntegrals(μlib::LibInt2Shell, νlib::LibInt2Shell, λlib::LibInt2Shell, σlib::LibInt2Shell)
    totaldim, maxprims, maxlqn = computeDimensions([μlib,νlib,λlib,σlib])
    engine = LibInt2EngineCoulomb(maxprims,maxlqn)

    result = copy(computeTensorBlockElectronRepulsionIntegrals(engine, μlib,νlib,λlib,σlib))

    destroy!(engine)
    return result
  end

  function computeBasisShellsLibInt2(basSet::BasisSet,geo::Geometry)
    shells = Vector{LibInt2Shell}()
    for atom in geo.atoms
      for contractedDefinition in basSet.definitions[atom.element]
        exponents = [prim.exponent for prim in contractedDefinition.primitives]
        coefficients = [prim.prefactor for prim in contractedDefinition.primitives]
        sh = LibInt2Shell(atom.position,contractedDefinition.lQuantumNumber,exponents,coefficients)
        push!(shells,sh)
      end
    end
    return shells
  end

else # if the libint2 is not available we have to stub the relevant functions
  include("LibInt2Module.stubs")
end

end # Module


text:
Base.map(fn, field::Field) = Base.broadcast(fn, field)

# useful operations
weighted_jacobian(field) = weighted_jacobian(axes(field))
weighted_jacobian(space::Spaces.AbstractSpace) =
    Spaces.local_geometry_data(space).WJ

# sum will give the integral over the field
function Base.sum(field::Union{Field, Base.Broadcast.Broadcasted{<:FieldStyle}})
    Base.reduce(
        RecursiveApply.radd,
        Base.Broadcast.broadcasted(
            RecursiveApply.rmul,
            weighted_jacobian(field),
            todata(field),
        ),
    )
end

function Base.sum(fn, field::Field)
    # Can't just call mapreduce as we need to weight _after_ applying the function
    Base.sum(Base.Broadcast.broadcasted(fn, field))
end

function LinearAlgebra.norm(field::Field, p::Real = 2)
    if p == 2
        # currently only one which supports structured types
        sqrt(sum(LinearAlgebra.norm_sqr, field))
    elseif p == 1
        sum(abs, field)
    elseif p == Inf
        error("Inf norm not yet supported")
    else
        sum(x -> x^p, field)^(1 / p)
    end
end

function Base.isapprox(
    x::Field,
    y::Field;
    atol::Real = 0,
    rtol::Real = Base.rtoldefault(eltype(parent(x)), eltype(parent(y)), atol),
    nans::Bool = false,
    norm::Function = LinearAlgebra.norm,
)
    d = norm(x .- y)
    return isfinite(d) && d <= max(atol, rtol * max(norm(x), norm(y)))
end


text:
using Random

"""
    random_candidate(problem::Problem)

Returns a randomized `Solution` based on specification in `problem.var_types`.
"""
function random_candidate(problem::Problem)
    return Solution(
        problem,
        [rand(t) for t in problem.var_types],
        false,
        zeros(problem.nobjs),
        0.0,
        0
    )
end

"""
    candidate_from_genome(problem::Problem, genome)

Returns a Solution object that packs genome as the solution's decision variable.
"""
function candidate_from_genome(problem::Problem, genome)
    return Solution(
        problem,
        genome,
        false,
        zeros(problem.nobjs),
        0.0,
        0
    )
end

function tournament_selector(population::Vector{Solution}, tournament_size::Int;
                             dominance::Function=compare_pareto_dominance)
    popsize = length(population)

    winner = rand(population)

    for i in 1:tournament_size-1
        challenger = rand(population)
        if dominance(winner, challenger) > 0
            winner = challenger
        end
    end

    return winner
end

function PM(parent::Solution, probability=1, distribution_index=20.0)
    child = copy(parent)
    problem = child.problem

    if typeof(probability) <: Int
        probability /= sum([typeof(t) == MOGA_Real for t in problem.var_types])
    end

    for i in 1:problem.nvars
        if typeof(problem.var_types[i]) == MOGA_Real
            if rand() < probability
                child.x[i] = pm_mutation(child.x[i],
                                         problem.var_types[1].min_value,
                                         problem.var_types[1].max_value,
                                         distribution_index)
                child.evaluated = false
            end
        end
    end
    return child
end

function PM!(s::Solution, probability=1, distribution_index=20.0)
    problem = s.problem
    if typeof(probability) <: Int
        probability /= sum([typeof(t) == MOGA_Real for t in problem.var_types])
    end

    for i in 1:problem.nvars
        if typeof(problem.var_types[i]) == MOGA_Real
            if rand() < probability
                s.x[i] = pm_mutation(s.x[i],
                                     problem.var_types[1].min_value,
                                     problem.var_types[1].max_value,
                                     distribution_index)
                s.evaluated = false
            end
        end
    end
end

function pm_mutation(x::Float64, lb::Float64, ub::Float64, di::Float64)
    u = rand()
    dx = ub - lb

    if u < 0.5
        bl = (x - lb) / dx
        b = 2*u + (1-2*u)*(1-bl)^(di+1)
        delta = b^(1/(di + 1)) - 1
    else
        bu = (ub-x)/dx
        b = 2*(1-u) + 2(u-0.5)*(1-bu)^(di+1)
        delta = 1 - b^(1/(di+1))
    end
    x = x + delta*dx
    x = clip(x, lb, ub)

    return x
end

function SBX(p1::Solution, p2::Solution, probability=1.0, distribution_index=15.0)
    c1 = copy(p1)
    c2 = copy(p2)

    if rand() <= probability
        problem = c1.problem
        nvars = problem.nvars
        for i = 1:nvars
            if typeof(problem.var_types[i]) == MOGA_Real
                if rand() <= 0.5
                    x1 = c1.x[i]
                    x2 = c2.x[i]
                    lb = problem.var_types[i].min_value
                    ub = problem.var_types[i].max_value

                    x1, x2 = sbx_crossover(x1, x2, lb, ub, distribution_index)

                    c1.x[i] = x1
                    c2.x[i] = x2
                    c1.evaluated = false
                    c2.evaluated = false
                end
            end
        end
    end
    return c1, c2
end

"Note that this function modifies s1 and s2"
function SBX!(s1::Solution, s2::Solution, probability=1.0, distribution_index=15.0)
    if rand() <= probability
        problem = s1.problem
        nvars = problem.nvars
        for i = 1:nvars
            if typeof(problem.var_types[i]) == MOGA_Real
                if rand() <= 0.5
                    x1 = s1.x[i]
                    x2 = s2.x[i]
                    lb = problem.var_types[i].min_value
                    ub = problem.var_types[i].max_value

                    x1, x2 = sbx_crossover(x1, x2, lb, ub, distribution_index)

                    s1.x[i] = x1
                    s2.x[i] = x2
                    s1.evaluated = false
                    s2.evaluated = false
                end
            end
        end
    end
end

function sbx_crossover(x1::Float64, x2::Float64, lb::Float64, ub::Float64, di::Float64)
    dx = abs(x2 - x1)

    if dx > eps()
        if x2 > x1
            y1 = x1
            y2 = x2
        else
            y1 = x2
            y2 = x1
        end

        rx = rand()

        # calc x1
        b = 1.0 / (1.0 + (2.0 * (y1 - lb) / (y2 - y1)))
        alpha = 2.0 - b ^ (di + 1.0)
        if rx <= 1.0 / alpha
            alpha *= rx
            betaq = alpha ^ (1.0/(di+1.0))
        else
            alpha = 1.0 / (2.0 - alpha*rx)
            betaq = alpha ^ (1.0/(di+1.0))
        end
        x1 = 0.5 * ((y1 + y2) - betaq * (y2 - y1))

        # calc x2
        b = 1.0 / (1.0 + (2.0 * (ub - y2) / (y2 - y1)))
        alpha = 2.0 - b^(di + 1.0)
        if rx <= 1.0 / alpha
            alpha = alpha * rx
            betaq = alpha ^ (1.0 / (di + 1.0))
        else
            alpha = alpha * rx
            alpha = 1.0 / (2.0 - alpha)
            betaq = alpha ^ (1.0 / (di + 1.0))
        end
        x2 = 0.5 * ((y1 + y2) + betaq * (y2 - y1))

        # randomly swap the values
        if Random.bitrand()[1]
            x1, x2 = x2, x1
        end

        x1 = clip(x1, lb, ub)
        x2 = clip(x2, lb, ub)

    end
    return x1, x2
end

function clip(x::Real, lb::Real, ub::Real)
    return max(lb, min(x, ub))
end

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna