text:
using Test, InMemoryDatasets, Random, PooledArrays, CategoricalArrays, DataAPI,
    Combinatorics, Unitful

import DataFrames as DF
const ≅ = isequal
const ≇ = !isequal

function _levels!(x::PooledArray, levels::AbstractVector)
    res = similar(x)
    copyto!(res, levels)
    copyto!(res, x)
end
_levels!(x::CategoricalArray, levels::AbstractVector) = levels!(x, levels)

function get_groups_for_groupby(gds)
    res = Vector{Int}(undef, nrow(parent(gds)))
    cnt = 1
    cnt2 = 1
    for i in 1:IMD._ngroups(gds)
        i == IMD._ngroups(gds) ? hi = nrow(parent(gds)) : hi = gds.starts[i+1]-1
        lo = gds.starts[i]
        for j in lo:hi
            res[cnt2] = cnt
            cnt2 += 1
        end
        cnt += 1
    end
    res[invperm(gds.perm)]
end

function isequal_unordered(gd,
                            gdsparts::AbstractVector)
    n = gd.lastvalid
    @assert n == length(gdsparts)
    remaining = Set(1:n)
    for i in 1:n
        for j in remaining
            if parent(gd)[IMD._get_perms(gd)[IMD.getindex_group(gd, i)], :] == gdsparts[j]
                pop!(remaining, j)
                break
            end
        end
    end
    isempty(remaining) || error("gd is not equal to provided groups")
end

@testset "parent" begin
    ds = Dataset(a = [1, 1, 2, 2], b = [5, 6, 7, 8])
    gds = gatherby(ds, :a)
    @test parent(gds) === ds
    df = DF.DataFrame(ds)
    gdf = DF.groupby(df, :a, sort = false)
    @test gds.groups == gdf.groups
    ds = Dataset(x1 = rand([missing, 2.2, 1.1], 1000), x2 = rand(1000))
    gds = gatherby(ds, :x1)
    @test parent(gds) === ds
    df = DF.DataFrame(ds)
    gdf = DF.groupby(df, :x1, sort = false)
    @test gds.groups == gdf.groups
    ds = Dataset(a = [1, 1, 2, 2], b = [5, 6, 7, 8])
    gds = groupby(ds, :a)
    @test parent(gds) === ds
    df = DF.DataFrame(ds)
    gdf = DF.groupby(df, :a, sort = true)
    @test get_groups_for_groupby(gds) == gdf.groups
    ds = Dataset(x1 = rand([missing, 2.2, 1.1], 1000), x2 = rand(1000))
    gds = groupby(ds, :x1)
    @test parent(gds) === ds
    df = DF.DataFrame(ds)
    gdf = DF.groupby(df, :x1, sort = true)
    @test get_groups_for_groupby(gds) == gdf.groups
    ds = Dataset(x1 = rand([missing; collect(1:1000000)], 1000), x2 = rand(1000))
    gds = groupby(ds, :x1)
    @test parent(gds) === ds
    df = DF.DataFrame(ds)
    gdf = DF.groupby(df, :x1, sort = true)
    @test get_groups_for_groupby(gds) == gdf.groups
    ds = Dataset(x1 = categorical(rand([missing; collect(1:1000000)], 1000)), x2 = rand(1000))
    gds = groupby(ds, :x1)
    @test parent(gds) === ds
    df = DF.DataFrame(ds)
    gdf = DF.groupby(df, :x1, sort = true)
    @test get_groups_for_groupby(gds) == gdf.groups
    ds = Dataset(x1 = PooledArray(rand([missing; collect(1:1000000)], 1000)), x2 = rand(1000))
    gds = groupby(ds, :x1)
    @test parent(gds) === ds
    df = DF.DataFrame(ds)
    gdf = DF.groupby(df, :x1, sort = true)
    @test get_groups_for_groupby(gds) == gdf.groups
    ds = Dataset(x1 = (rand([missing; collect(1:1000000)], 1000)) .* 1.1, x2 = rand(1000))
    gds = groupby(ds, :x1)
    @test parent(gds) === ds
    df = DF.DataFrame(ds)
    gdf = DF.groupby(df, :x1, sort = true)
    @test get_groups_for_groupby(gds) == gdf.groups
    ds = Dataset(x1 = (rand([missing; collect(1:1000000)], 1000)) .* 1.1, x2 = rand(1000))
    gds = gatherby(ds, :x1)
    @test parent(gds) === ds
    df = DF.DataFrame(ds)
    gdf = DF.groupby(df, :x1, sort = false)
    @test gds.groups == gdf.groups
    ds = Dataset(x1 = categorical(rand([missing; collect(1:1000000)], 1000)), x2 = rand(1000))
    gds = gatherby(ds, :x1)
    @test parent(gds) === ds
    df = DF.DataFrame(ds)
    gdf = DF.groupby(df, :x1, sort = false)
    @test gds.groups == gdf.groups
    ds = Dataset(rand(1:1000000, 1000, 10), :auto)
    gds = gatherby(ds, :)
    @test parent(gds) === ds
    df = DF.DataFrame(ds)
    gdf = DF.groupby(df, :, sort = false)
    @test gds.groups == gdf.groups
    ds = Dataset(rand(1:1000000, 1000, 10), :auto)
    gds = groupby(ds, :)
    @test parent(gds) === ds
    df = DF.DataFrame(ds)
    gdf = DF.groupby(df, :, sort = true)
    @test get_groups_for_groupby(gds) == gdf.groups
end

@testset "consistency" begin
    ds = Dataset(a = [1, 1, 2, 2], b = [5, 6, 7, 8], c = 1:4)
    push!(ds.c.val, 5)
    @test_throws AssertionError gatherby(ds, :a)

    ds = Dataset(a = [1, 1, 2, 2], b = [5, 6, 7, 8], c = 1:4)
    push!(IMD._columns(ds), ds[:, :a])
    @test_throws AssertionError gatherby(ds, :a)

    ds = Dataset(a = [1, 1, 2, 2], b = [5, 6, 7, 8], c = 1:4)
    push!(ds.c.val, 5)
    @test_throws AssertionError groupby(ds, :a)

    ds = Dataset(a = [1, 1, 2, 2], b = [5, 6, 7, 8], c = 1:4)
    push!(IMD._columns(ds), ds[:, :a])
    @test_throws AssertionError groupby(ds, :a)
end

@testset "general usage" begin
    Random.seed!(1)
    ds = Dataset(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)))

   df = DF.DataFrame(ds)
   @test ds[IMD._get_perms(groupby(ds, 1:2)), :] == Dataset(df[DF.groupby(df, 1:2, sort = true).idx, :])
   @test ds[IMD._get_perms(gatherby(ds, 1:2)), :] == Dataset(df[DF.groupby(df, 1:2, sort = false).idx, :])
   @test ds[IMD._get_perms(groupby(ds, [2, 3, 1])), :] == Dataset(df[DF.groupby(df, [2, 3, 1], sort = true).idx, :])
   @test ds[IMD._get_perms(gatherby(ds, [2, 3, 1])), :] == Dataset(df[DF.groupby(df, [2, 3, 1], sort = false).idx, :])
   @test ds[IMD._get_perms(groupby(ds, 1:3)), :] == Dataset(df[DF.groupby(df, 1:3, sort = true).idx, :])
   @test ds[IMD._get_perms(gatherby(ds, 1:3)), :] == Dataset(df[DF.groupby(df, 1:3, sort = false).idx, :])
    # test number of potential combinations higher than typemax(Int32)
    N = 2000
    ds2 = Dataset(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)))
    df2 = DF.DataFrame(ds2)
    ds2b = mapcols(ds2, Vector{Int}, :)
    @test ds2[IMD._get_perms(groupby(ds2, 1:2)), :] == Dataset(df2[DF.groupby(df2, 1:2, sort = true).idx, :])
    @test ds2[IMD._get_perms(gatherby(ds2, 1:2)), :] == Dataset(df2[DF.groupby(df2, 1:2, sort = false).idx, :])
    @test ds2[IMD._get_perms(groupby(ds2, [2, 3, 1])), :] == Dataset(df2[DF.groupby(df2, [2, 3, 1], sort = true).idx, :])
    @test ds2[IMD._get_perms(gatherby(ds2, [2, 3, 1])), :] == Dataset(df2[DF.groupby(df2, [2, 3, 1], sort = false).idx, :])
    @test ds2[IMD._get_perms(groupby(ds2, 1:3)), :] == Dataset(df2[DF.groupby(df2, 1:3, sort = true).idx, :])
    @test ds2[IMD._get_perms(gatherby(ds2, 1:3)), :] == Dataset(df2[DF.groupby(df2, 1:3, sort = false).idx, :])
    @test ds2b[IMD._get_perms(groupby(ds2b, 1:2)), :] == Dataset(df2[DF.groupby(df2, 1:2, sort = true).idx, :])
    @test ds2b[IMD._get_perms(gatherby(ds2b, 1:2)), :] == Dataset(df2[DF.groupby(df2, 1:2, sort = false).idx, :])
    @test ds2b[IMD._get_perms(groupby(ds2b, [2, 3, 1])), :] == Dataset(df2[DF.groupby(df2, [2, 3, 1], sort = true).idx, :])
    @test ds2b[IMD._get_perms(gatherby(ds2b, [2, 3, 1])), :] == Dataset(df2[DF.groupby(df2, [2, 3, 1], sort = false).idx, :])
    @test ds2b[IMD._get_perms(groupby(ds2b, 1:3)), :] == Dataset(df2[DF.groupby(df2, 1:3, sort = true).idx, :])
    @test ds2b[IMD._get_perms(gatherby(ds2b, 1:3)), :] == Dataset(df2[DF.groupby(df2, 1:3, sort = false).idx, :])
    x = CategoricalArray(collect(1:20))
    ds = Dataset(v1=x, v2=x)
    df = DF.DataFrame(v1 = x, v2 = x)
    @test groupby!(ds, 1, rev = true) == Dataset(sort!(df, DF.order(:v1, rev = true)))
end

@testset "grouping for PAs - From DataFrames.jl" begin
    xv = ["A", missing, "B", "B", "A", "B", "A", "A"]
    yv = ["B", "A", "A", missing, "A", missing, "A", "A"]
    xvars = (xv,
             categorical(xv),
             levels!(categorical(xv), ["A", "B", "X"]),
             levels!(categorical(xv), ["X", "B", "A"]),
             _levels!(PooledArray(xv), ["A", "B", missing]),
             _levels!(PooledArray(xv), ["B", "A", missing, "X"]),
             _levels!(PooledArray(xv), [missing, "X", "A", "B"]))
    yvars = (yv,
             categorical(yv),
             levels!(categorical(yv), ["A", "B", "X"]),
             levels!(categorical(yv), ["B", "X", "A"]),
             _levels!(PooledArray(yv), ["A", "B", missing]),
             _levels!(PooledArray(yv), [missing, "A", "B", "X"]),
             _levels!(PooledArray(yv), ["B", "A", "X", missing]))
    for x in xvars, y in yvars
        ds = Dataset(Key1=x, Key2=y, Value=1:8)

        @testset "gatherby" begin
            gd = gatherby(ds, :Key1)
            @test gd.lastvalid == 3
            @test isequal_unordered(gd, [Dataset(Key1="A", Key2=["B", "A", "A", "A"], Value=[1, 5, 7, 8]),
                                        Dataset(Key1="B", Key2=["A", missing, missing], Value=[3, 4, 6]),
                                        Dataset(Key1=missing, Key2="A", Value=2)])

            gd = gatherby(ds, [:Key1, :Key2])
            @test gd.lastvalid == 5
            @test isequal_unordered(gd, [Dataset(Key1="A", Key2="A", Value=[5, 7, 8]),
                                Dataset(Key1="A", Key2="B", Value=1),
                                Dataset(Key1="B", Key2="A", Value=3),
                                Dataset(Key1="B", Key2=missing, Value=[4, 6]),
                                Dataset(Key1=missing, Key2="A", Value=2)])
        end

        @testset "gatherby with dropmissing" begin
            gd = gatherby(dropmissing(ds, :Key1, view=true), :Key1)
            @test gd.lastvalid == 2
            @test isequal_unordered(gd, [Dataset(Key1="A", Key2=["B", "A", "A", "A"], Value=[1, 5, 7, 8]),
            Dataset(Key1="B", Key2=["A", missing, missing], Value=[3, 4, 6])])

            gd = gatherby(dropmissing(ds, [:Key1, :Key2], view=true), [:Key1, :Key2])
            @test gd.lastvalid == 3
            @test isequal_unordered(gd, [Dataset(Key1="A", Key2="A", Value=[5, 7, 8]),
            Dataset(Key1="A", Key2="B", Value=1),
            Dataset(Key1="B", Key2="A", Value=3)])

        end

        @testset "groupby" begin
            gd = groupby(ds, :Key1)
            @test gd.lastvalid == 3
            @test isequal_unordered(gd, [
                Dataset(Key1="A", Key2=["B", "A", "A", "A"], Value=[1, 5, 7, 8]),
                Dataset(Key1="B", Key2=["A", missing, missing], Value=[3, 4, 6]),
                Dataset(Key1=missing, Key2="A", Value=2)
            ])

            gd = groupby(ds, [:Key1, :Key2])
            @test gd.lastvalid == 5
            @test isequal_unordered(gd, [
                Dataset(Key1="A", Key2="A", Value=[5, 7, 8]),
                Dataset(Key1="A", Key2="B", Value=1),
                Dataset(Key1="B", Key2="A", Value=3),
                Dataset(Key1="B", Key2=missing, Value=[4, 6]),
                Dataset(Key1=missing, Key2="A", Value=2)
            ])
        end

        @testset "groupby with dropmissing" begin
            gd = groupby(dropmissing(ds, :Key1, view=true), :Key1)
            @test gd.lastvalid == 2
            @test isequal_unordered(gd, [
                Dataset(Key1="A", Key2=["B", "A", "A", "A"], Value=[1, 5, 7, 8]),
                Dataset(Key1="B", Key2=["A", missing, missing], Value=[3, 4, 6])
            ])

            gd = groupby(dropmissing(ds, [:Key1, :Key2], view=true), [:Key1, :Key2])
            @test gd.lastvalid == 3
            @test isequal_unordered(gd, [
                Dataset(Key1="A", Key2="A", Value=[5, 7, 8]),
                Dataset(Key1="A", Key2="B", Value=1),
                Dataset(Key1="B", Key2="A", Value=3)
            ])
        end
    end

    @test gatherby(Dataset(x=[missing]), :x).groups ==
        gatherby(Dataset(x=Union{Int, Missing}[missing]), :x).groups ==
        gatherby(Dataset(x=Union{String, Missing}[missing]), :x).groups ==
        gatherby(Dataset(x=Any[missing]), :x).groups == [1]
    @test isempty(gatherby(dropmissing(Dataset(x=[missing])), :x).groups)

end

@testset "combine - set 1" begin
    Random.seed!(1)
    # 5 is there to ensure we test a single-row group
    ds = Dataset(a=[rand([1:4;missing], 19); 5],
                   x1=rand(1:100, 20),
                   x2=rand(1:3, 20) + im*rand(1:3, 20),
                   x3=repeat(1:2, 10) .* u"m")
    combine_out1 = combine(groupby(ds, 1), :x1=>sum=>:x1_sum, :x1=>(x->sum(x))=>:x1_sum2)
    combine_out2 = combine(gatherby(ds, 1), :x1=>sum=>:x1_sum, :x1=>(x->sum(x))=>:x1_sum2)
    combine_out3 = combine(gatherby(dropmissing(ds, 1, view=true), 1), :x1=>sum=>:x1_sum, :x1=>(x->sum(x))=>:x1_sum2)
    combine_out4 = combine(gatherby(ds, 1), :x1=>maximum, :x1=>minimum, 2:3=>byrow(-))
    combine_out5 = combine(gatherby(dropmissing(ds, 1, view=true), 1), :x1=>maximum, :x1=>minimum, 2:3=>byrow(-))
    combine_out6 = combine(gatherby(ds, 1), :x1=>IMD.n)

    @test combine_out1 == combine(groupby(ds, 1), :x1=>sum=>:x1_sum, :x1=>(x->sum(x))=>:x1_sum2, threads = false)
    @test combine_out2 == combine(gatherby(ds, 1), :x1=>sum=>:x1_sum, :x1=>(x->sum(x))=>:x1_sum2, threads = false)
    @test combine_out3 == combine(gatherby(dropmissing(ds, 1, view=true), 1), :x1=>sum=>:x1_sum, :x1=>(x->sum(x))=>:x1_sum2, threads = false)
    @test combine_out4 == combine(gatherby(ds, 1), :x1=>maximum, :x1=>minimum, 2:3=>byrow(-), threads = false)
    @test combine_out5 == combine(gatherby(dropmissing(ds, 1, view=true), 1), :x1=>maximum, :x1=>minimum, 2:3=>byrow(-), threads = false)
    @test combine_out6 == combine(gatherby(ds, 1), :x1=>IMD.n, threads = false)

    @test combine_out1[!, 2] == combine_out1[!, 3]
    @test combine_out2[!, 2] == combine_out2[!, 3]
    @test combine_out3[!, 2] == combine_out3[!, 3]
    @test all(x->x>=0, combine_out4[!, "row_-"])
    @test all(x->x>=0, combine_out5[!, "row_-"])
    @test dropmissing(combine_out4, 1) == combine_out5
    @test combine_out6[nrow(combine_out6), 1] == 5
    combine_out7 = combine(groupby(ds, :x3), :x1=>mean)
    combine_out8 = combine(gatherby(ds, :x3), :x2=>IMD.n)
    combine_out9 = combine(groupby(ds, :x3), :x3=>maximum, :x3=>minimum)
    combine_out7_v = combine(groupby(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), :x3), :x1=>mean)
    combine_out8_v = combine(gatherby(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), :x3), :x2=>IMD.n)
    combine_out8_v = combine_out8_v[nrow(combine_out8_v):-1:1, :]
    combine_out9_v = combine(groupby(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), :x3), :x3=>maximum, :x3=>minimum)

    @test combine_out7 == combine(groupby(ds, :x3), :x1=>mean, threads = false)
    @test combine_out8 == combine(gatherby(ds, :x3), :x2=>IMD.n, threads = false)
    @test combine_out9 == combine(groupby(ds, :x3), :x3=>maximum, :x3=>minimum, threads = false)
    @test combine_out7_v == combine(groupby(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), :x3), :x1=>mean, threads = false)
    @test combine_out9_v == combine(groupby(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), :x3), :x3=>maximum, :x3=>minimum, threads = false)

    @test combine_out7[!, 1] == [1 * u"m", 2 * u"m"]
    @test combine_out8[!, 2] == [10, 10]
    @test combine_out9[!, 2] == combine_out9[!, 1]
    @test combine_out9[!, 3] == combine_out9[!, 1]

    @test all(byrow(compare(combine_out7, combine_out7_v, on = names(combine_out7)), all, :))
    @test all(byrow(compare(combine_out8, combine_out8_v, on = names(combine_out8)), all, :))
    @test all(byrow(compare(combine_out9, combine_out9_v, on = names(combine_out9)), all, :))

    ds = Dataset(rand(Int, 100, 1), :auto)
    @test sort(combine(gatherby(ds, 1), 1=>length), :) == combine(groupby(ds, 1), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false), 1=>length), :) == combine(groupby(ds, 1, stable = false), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false), 1=>(x->length(x))=>:length_x1), :) == combine(groupby(ds, 1, stable = false, alg = QuickSort), 1=>length)

    @test sort(combine(gatherby(ds, 1, threads = false), 1=>length, threads = false), :) == combine(groupby(ds, 1), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false, threads = false), 1=>length, threads = false), :) == combine(groupby(ds, 1, stable = false), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false, threads = false), 1=>(x->length(x))=>:length_x1, threads = false), :) == combine(groupby(ds, 1, stable = false, alg = QuickSort), 1=>length)

    ds = Dataset(rand(Int128, 100, 1), :auto)
    @test sort(combine(gatherby(ds, 1), 1=>length), :) == combine(groupby(ds, 1), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false), 1=>length), :) == combine(groupby(ds, 1, stable = false), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false), 1=>(x->length(x))=>:length_x1), :) == combine(groupby(ds, 1, stable = false, alg = QuickSort), 1=>length)

    @test sort(combine(gatherby(ds, 1, threads = false), 1=>length, threads = false), :) == combine(groupby(ds, 1), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false, threads = false), 1=>length, threads = false), :) == combine(groupby(ds, 1, stable = false), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false, threads = false), 1=>(x->length(x))=>:length_x1, threads = false), :) == combine(groupby(ds, 1, stable = false, alg = QuickSort), 1=>length)

    ds = Dataset(rand(UInt64, 100, 1), :auto)
    @test sort(combine(gatherby(ds, 1), 1=>length), :) == combine(groupby(ds, 1), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false), 1=>length), :) == combine(groupby(ds, 1, stable = false), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false), 1=>(x->length(x))=>:length_x1), :) == combine(groupby(ds, 1, stable = false, alg = QuickSort), 1=>length)

    @test sort(combine(gatherby(ds, 1, threads = false), 1=>length, threads = false), :) == combine(groupby(ds, 1), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false, threads = false), 1=>length, threads = false), :) == combine(groupby(ds, 1, stable = false), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false, threads = false), 1=>(x->length(x))=>:length_x1, threads = false), :) == combine(groupby(ds, 1, stable = false, alg = QuickSort), 1=>length)

    ds = Dataset(rand(Int8, 100, 1), :auto)
    @test sort(combine(gatherby(ds, 1), 1=>length), :) == combine(groupby(ds, 1), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false), 1=>length), :) == combine(groupby(ds, 1, stable = false), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false), 1=>(x->length(x))=>:length_x1), :) == combine(groupby(ds, 1, stable = false, alg = QuickSort), 1=>length)

    @test sort(combine(gatherby(ds, 1, threads = false), 1=>length, threads = false), :) == combine(groupby(ds, 1), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false, threads = false), 1=>length, threads = false), :) == combine(groupby(ds, 1, stable = false), 1=>length)
    @test sort(combine(gatherby(ds, 1, stable = false, threads = false), 1=>(x->length(x))=>:length_x1, threads = false), :) == combine(groupby(ds, 1, stable = false, alg = QuickSort), 1=>length)

end

# TODO windows on CI stuck here ??
if !Base.Sys.iswindows()
    @testset "combine - set 2" begin
        ds = Dataset(a=[rand([1:4;missing], 19); 5],
                       x1=rand(1:100, 20),
                       x2=rand(1:3, 20) + im*rand(1:3, 20),
                       x3=repeat(1:2, 10) .* u"m")
        c1_IMD = combine(gatherby(ds, [:a, :x3]), :x1=>sum=>:o1, :x1=>argmax=>:o2, :x1=>sort=>:o3)
        @test c1_IMD == combine(gatherby(ds, [:a, :x3], threads = false), :x1=>sum=>:o1, :x1=>argmax=>:o2, :x1=>sort=>:o3, threads = false)
        c1_DF = DF.combine(DF.groupby(DF.DataFrame(ds), [:a, :x3], sort = false), :x1=>sum=>:o1, :x1=>argmax=>:o2, :x1=>sort=>:o3)

        @test all(byrow(compare(c1_IMD, Dataset(c1_DF), on = names(c1_IMD)), all, :))
        c1_IMD = combine(groupby(ds, [:a, :x3]), :x1=>sum=>:o1, :x1=>argmax=>:o2, :x1=>sort=>:o3)
        @test c1_IMD == combine(groupby(ds, [:a, :x3], threads = false), :x1=>sum=>:o1, :x1=>argmax=>:o2, :x1=>sort=>:o3, threads = false)

        c1_DF = DF.combine(DF.groupby(DF.DataFrame(ds), [:a, :x3], sort = true), :x1=>sum=>:o1, :x1=>argmax=>:o2, :x1=>sort=>:o3)
        @test all(byrow(compare(c1_IMD, Dataset(c1_DF), on = names(c1_IMD)), all, :))

        for r in 0:.1:1
            ds = Dataset(rand(1:100, 100, 10), :auto)
            modify!(ds, [1,5,3] => x->Float32.(x))
            map!(ds, x->rand()<r ? missing : x, :)
            c1 = combine(groupby(ds, [:x1,:x3]), :x10=>sort, :x4=>sum)
            c2 = combine(groupby(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), [:x1, :x3]), :x10=>sort, :x4=>sum)

            @test c1 == combine(groupby(ds, [:x1,:x3], threads = false), :x10=>sort, :x4=>sum, threads = false)
            @test c2 == combine(groupby(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), [:x1, :x3], threads = false), :x10=>sort, :x4=>sum, threads = false)
            @test c1 == c2
            c1 = combine(gatherby(ds, [:x1,:x3]), :x10=>sort, :x4=>sum)
            c2 = combine(gatherby(view(ds, nrow(ds):-1:1, [10, 1, 3, 4]), [:x1, :x3]), :x10=>sort, :x4=>sum)
            @test c1 == combine(gatherby(ds, [:x1,:x3], threads = false), :x10=>sort, :x4=>sum, threads = false)
            @test c2 == combine(gatherby(view(ds, nrow(ds):-1:1, [10, 1, 3, 4]), [:x1, :x3], threads = false), :x10=>sort, :x4=>sum, threads = false)
            @test sort(c1, :) == sort(c2, :)
        end
        ds = Dataset(rand(1000, 2), :auto)
        insertcols!(ds, 1, :g=>rand(1:10, nrow(ds)))
        c1_IMD = combine(groupby(ds, 1), (2, 3)=>cor=>:cor, 3=>sum=>:sum, 2 => mean => :mean)
        @test c1_IMD == combine(groupby(ds, 1, threads = false), (2, 3)=>cor=>:cor, 3=>sum=>:sum, 2 => mean => :mean, threads = false)
        df = DF.DataFrame(ds)
        c1_DF = DF.combine(DF.groupby(df, 1, sort = true), 2:3=>cor=>:cor, 3=>sum=>:sum, 2=>mean=>:mean)
        @test all(byrow(compare(c1_IMD, Dataset(c1_DF), eq = isapprox), all, :))
        c1_IMD = combine(groupby(view(ds, nrow(ds):-1:1, ncol(ds):-1:1), 3), (2,1)=>cor=>:cor, 1=>sum=>:sum, 2 => mean => :mean)
        df = DF.DataFrame(ds)
        c1_DF = DF.combine(DF.groupby(df, 1, sort = true), 2:3=>cor=>:cor, 3=>sum=>:sum, 2=>mean=>:mean)
        @test all(byrow(compare(c1_IMD, Dataset(c1_DF), eq = isapprox), all, :))
    end
end

@testset "combine - set 3" begin
    ds = Dataset(x1 = [1,2,1,2,2,1], x2=[1.0,missing,1.1,1.1,1.1,1.1],y=100:100:600.0)
    c1 = combine(gatherby(ds, 1), :x2=>IMD.sum)
    c2 = combine(groupby(ds, 1, stable = false), :y=>mean)
    c3 = combine(gatherby(ds,[:x2], isgathered=true), :y=>IMD.maximum)
    c4 = combine(groupby(ds, r"x", stable = false), (:y, :x2)=>(x,y)->y[IMD.argmax(x)])

    @test c1 == Dataset(x1=[1,2], sum_x2=[3.2,2.2])
    @test c2 == Dataset(x1=[1,2], mean_y=[1000/3, 1100/3])
    @test c3 == Dataset(x2=[1,missing,1.1], maximum_y=[100.0, 200, 600])
    @test c4 == Dataset(x1=[1,1,2,2], x2=[1.0,1.1,1.1,missing], function_y_x2=[1,1.1,1.1,missing])
    combinefmt(x)=ismissing(x) ? missing : 0
    setformat!(ds, :x2=>combinefmt)
    c1 = combine(groupby(ds, :x2), :y=>IMD.maximum)
    c2 = combine(groupby(ds, :x2, mapformats = false), :y=>IMD.maximum)
    c3 = combine(gatherby(ds, :x2), (:x1,:y)=>(x,y)->(IMD.maximum(x), IMD.minimum(y)))

    @test c1 == combine(groupby(ds, :x2, threads = false), :y=>IMD.maximum, threads = false)
    @test c2 == combine(groupby(ds, :x2, mapformats = false, threads = false), :y=>IMD.maximum, threads = false)
    @test c3 == combine(gatherby(ds, :x2, threads = false), (:x1,:y)=>(x,y)->(IMD.maximum(x), IMD.minimum(y)), threads = false)

    @test all(byrow(compare(c1, Dataset(x2=[0, missing], maximum_y=[600.0, 200]), mapformats=true), all))
    @test all(byrow(compare(c2, Dataset(x2=[0, 0, missing], maximum_y=[100, 600.0, 200]), mapformats=true), all))
    @test all(byrow(compare(c3, Dataset(x2=[0,missing], function_x1_y=[(2, 100.0), (2, 200.0)]), mapformats=true), all))
    ds = Dataset(x1 = [1,2,1,2,2,1], x2=[1.0,missing,1.1,1.1,1.1,1.1],y=100:100:600.0)
    @test byrow(compare(combine(gatherby(ds, 1), :x2=>[IMD.sum, IMD.maximum],2:3=>byrow(+)=>:row), Dataset(x1=[1,2], sum_x2=[3.2,2.2], maximum_x2=[1.1,1.1], row = [4.3, 3.3] ), eq = isapprox), all)|>all
    @test byrow(compare(combine(gatherby(ds, 1), :x2=>[x->IMD.sum(x), IMD.maximum],2:3=>byrow(+)=>:row), Dataset(x1=[1,2], function_x2=[3.2,2.2], maximum_x2=[1.1,1.1], row = [4.3, 3.3] ), eq = isapprox), all)|>all
    @test byrow(compare(combine(groupby(ds, 1), :x2=>[IMD.sum, IMD.maximum],2:3=>byrow(+)=>:row), Dataset(x1=[1,2], sum_x2=[3.2,2.2], maximum_x2=[1.1,1.1], row = [4.3, 3.3] ), eq = isapprox), all)|>all

    @test byrow(compare(combine(gatherby(view(ds, :, [2,1]), :x1), :x2=>[IMD.sum, IMD.maximum],2:3=>byrow(+)=>:row), Dataset(x1=[1,2], sum_x2=[3.2,2.2], maximum_x2=[1.1,1.1], row = [4.3, 3.3] ), eq = isapprox), all)|>all
    @test byrow(compare(combine(gatherby(view(ds, :, [2,1]), :x1), :x2=>[x->IMD.sum(x), IMD.maximum],2:3=>byrow(+)=>:row), Dataset(x1=[1,2], function_x2=[3.2,2.2], maximum_x2=[1.1,1.1], row = [4.3, 3.3] ), eq = isapprox), all)|>all
    @test byrow(compare(combine(groupby(view(ds, :, [2,1]), :x1), :x2=>[IMD.sum, IMD.maximum],2:3=>byrow(+)=>:row), Dataset(x1=[1,2], sum_x2=[3.2,2.2], maximum_x2=[1.1,1.1], row = [4.3, 3.3] ), eq = isapprox), all)|>all
end

@testset "combine - set4 - checking normalize_combine" begin
    ds = Dataset(g1 = [1,2,1,2,2,1], g2 = [2,3,1,1,2,3], x2=[1.0,missing,1.1,1.1,1.1,1.1],y=100:100:600.0)
    @test combine(groupby(ds, 2), [1,3] => [IMD.maximum, IMD.minimum])== Dataset(g2 = [1,2,3], maximum_g1 = [2,2,2], minimum_g1 = [1,1,1], maximum_x2 = [1.1,1.1,1.1], minimum_x2 = [1.1,1,1.1])
    @test combine(gatherby(ds, 2), [1,3] => [IMD.maximum, IMD.minimum])== Dataset(g2 = [2,3,1], maximum_g1 = [2,2,2], minimum_g1 = [1,1,1], maximum_x2 = [1.1,1.1,1.1], minimum_x2 = [1.,1.1,1.1])
    @test combine(groupby(ds, [2,1]), r"x"=>IMD.sum) == Dataset(g2 = [1,1,2,2,3,3], g1=[1,2,1,2,1,2], sum_x2 = [1.1,1.1,1,1.1,1.1, missing])
    @test combine(gatherby(ds, [2,1]), r"x"=>IMD.sum) == Dataset(g2 = [2,3,1,1,2,3], g1=[1,2,1,2,2,1], sum_x2 = [1, missing, 1.1,1.1,1.1,1.1])
    groupby!(ds, [2,1])
    @test combine(ds, r"x"=>IMD.sum) == Dataset(g2 = [1,1,2,2,3,3], g1=[1,2,1,2,1,2], sum_x2 = [1.1,1.1,1,1.1,1.1, missing])
    @test combine(groupby(ds, 2), [1,3] => [IMD.maximum, IMD.minimum])== Dataset(g2 = [1,2,3], maximum_g1 = [2,2,2], minimum_g1 = [1,1,1], maximum_x2 = [1.1,1.1,1.1], minimum_x2 = [1.1,1,1.1])
    @test combine(groupby(ds, [2,1]), r"x"=>IMD.sum) == Dataset(g2 = [1,1,2,2,3,3], g1=[1,2,1,2,1,2], sum_x2 = [1.1,1.1,1,1.1,1.1, missing])

    @test combine(groupby(ds, 1), 3:4 => IMD.sum,  2:3=> byrow.([maximum, minimum])) == Dataset(g1 = [1,2], sum_x2 = [3.2, 2.2], sum_y = [1000, 1100], row_maximum = [1000,1100], row_minimum=[3.2,2.2])
    @test combine(groupby(ds, 1), 3:4 => IMD.sum,  2:3=> byrow.([maximum, minimum]), 4:5=>byrow(sum)) == Dataset(g1 = [1,2], sum_x2 = [3.2, 2.2], sum_y = [1000, 1100], row_maximum = [1000,1100], row_minimum=[3.2,2.2], row_sum = [1003.2, 1102.2])

    @test combine(groupby(view(ds, :, [4,3,2,1]), 4), [2,1] => IMD.sum,  2:3=> byrow.([maximum, minimum])) == Dataset(g1 = [1,2], sum_x2 = [3.2, 2.2], sum_y = [1000, 1100], row_maximum = [1000,1100], row_minimum=[3.2,2.2])
    @test combine(groupby(view(ds,:, [4,3,2,1]), 4), [:x2, :y] => IMD.sum,  2:3=> byrow.([maximum, minimum]), 4:5=>byrow(sum)) == Dataset(g1 = [1,2], sum_x2 = [3.2, 2.2], sum_y = [1000, 1100], row_maximum = [1000,1100], row_minimum=[3.2,2.2], row_sum = [1003.2, 1102.2])

end

@testset "modifying and combining views" begin
    ds = Dataset(x = [3,1,2,2,missing,3,3], y = [1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], z = [11,15,7,-11,12,0,0])
    sds1 = dropmissing(ds, 2, view = true)
    sds2 = view(ds, [2,1,1,3,5,6,7,4,4], [:z, :x])

    @test combine(groupby(sds1, :x), :z => IMD.sum) == Dataset(x = [2, 3], sum_z = [-4, 11])
    @test combine(gatherby(sds1, :x), :z => IMD.sum) == Dataset(x = [3, 2], sum_z = [11, -4])
    @test combine(groupby(sds2, :x), :z => IMD.sum) == Dataset(x = [1,2, 3, missing], sum_z = [15, -15, 22, 12])
    @test combine(gatherby(sds2, :x), :z => IMD.sum) == Dataset(x = [1, 3, 2, missing], sum_z = [15, 22, -15, 12])

    @test combine(groupby(sds1, :x), :z => (x->IMD.sum(x))=>:sum_z) == Dataset(x = [2, 3], sum_z = [-4, 11])
    @test combine(gatherby(sds1, :x), :z => (x->IMD.sum(x))=>:sum_z) == Dataset(x = [3, 2], sum_z = [11, -4])
    @test combine(groupby(sds2, :x), :z => (x->IMD.sum(x))=>:sum_z) == Dataset(x = [1,2, 3, missing], sum_z = [15, -15, 22, 12])
    @test combine(gatherby(sds2, :x), :z => (x->IMD.sum(x))=>:sum_z) == Dataset(x = [1, 3, 2, missing], sum_z = [15, 22, -15, 12])

    @test combine(groupby(sds1, :x), :y=>(sort!)=>:s_y) == Dataset(x=[2,2,3,3,3], s_y=[-3.0,-1.0,1.1,4.0,5.0])
    @test combine(gatherby(sds1, :x), :y=>(sort!)=>:s_y) == Dataset(x=[3,3,3,2,2], s_y=[1.1,4.0,5.0,-3.0,-1.0])

    @test combine(groupby(sds2, :x), :z=>(sort!)=>:s_z) == Dataset(x=[1,2,2,2,3,3,3,3,missing], s_z=[15,-11,-11,7,0,0,11,11,12])
    @test combine(gatherby(sds2, :x), :z=>(sort!)=>:s_z) == Dataset(x=[1,3,3,3,3,2,2,2,missing], s_z=[15,0,0,11,11,-11,-11,7,12])

    @test combine(groupby(sds1, :x), (:y,:z)=>cor) == Dataset(x=[2,3],cor_y_z=[1.0, -0.9690582663799521])
    @test combine(gatherby(sds1, :x), (:y,:z)=>cor) == Dataset(x=[3,2],cor_y_z=[-0.9690582663799521, 1.0])

    @test combine(groupby(sds2, :x), (1,2)=>((x,y)->IMD.maximum(x)/length(y))=>:q) == Dataset(x=[1,2,3,missing], q=[15.0, 7/3,11/4,12.0])
    @test combine(gatherby(sds2, :x), (1,2)=>((x,y)->IMD.maximum(x)/length(y))=>:q) == Dataset(x=[1,3,2,missing], q=[15.0, 11/4,7/3,12.0])

    @test combine(groupby(sds1, :x), :y=>IMD.maximum, :z=>IMD.maximum, 2:3=>byrow(-)=>:q) == Dataset(x=[2,3], maximum_y=[-1.0,5],maximum_z=[7, 11], q=[-8,-6.0])
    @test combine(gatherby(sds1, :x), :y=>IMD.maximum, :z=>IMD.maximum, 2:3=>byrow(-)=>:q) == Dataset(x=[3,2], maximum_y=reverse([-1.0,5]),maximum_z=reverse([7, 11]), q=reverse([-8,-6.0]))

    @test combine(groupby(sds2, :x), :z=>(sort!)=>:s_z, :x=>IMD.maximum, :z=>IMD.minimum, 3:4=>byrow(-)=>:q) == Dataset([Union{Missing, Int64}[1, 2, 2, 2, 3, 3, 3, 3, missing], Union{Missing, Int64}[15, -11, -11, 7, 0, 0, 11, 11, 12], Union{Missing, Int64}[1, 2, 2, 2, 3, 3, 3, 3, missing], Union{Missing, Int64}[15, -11, -11, -11, 0, 0, 0, 0, 12], Union{Missing, Int64}[-14, 13, 13, 13, 3, 3, 3, 3, missing]], [:x, :s_z, :maximum_x, :minimum_z, :q])

    @test combine(gatherby(sds2, :x), :z=>(sort!)=>:s_z, :x=>IMD.maximum, :z=>IMD.minimum, 3:4=>byrow(-)=>:q) == Dataset([Union{Missing, Int64}[1, 3, 3, 3, 3, 2, 2, 2, missing], Union{Missing, Int64}[15, 0, 0, 11, 11, -11, -11, 7, 12], Union{Missing, Int64}[1, 3, 3, 3, 3, 2, 2, 2, missing], Union{Missing, Int64}[15, 0, 0, 0, 0, -11, -11, -11, 12], Union{Missing, Int64}[-14, 3, 3, 3, 3, 13, 13, 13, missing]], [:x, :s_z, :maximum_x, :minimum_z, :q])

    @test ds == Dataset(x = [3,1,2,2,missing,3,3], y = [1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], z = [11,15,7,-11,12,0,0])

    @test_throws ArgumentError combine(gatherby(sds2, :x), :q=>IMD.sum)
    @test_throws ArgumentError combine(gatherby(sds2, :x), :y=>IMD.sum)

    modify!(groupby(sds1, :x), :y=>IMD.maximum=>:q)
    @test ds == Dataset(x = [3,1,2,2,missing,3,3], y = [1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], z = [11,15,7,-11,12,0,0], q=[5.0, missing, -1,-1,missing,5,5])
    modify!(gatherby(sds1, :x), :y=>IMD.maximum=>:q2)
    @test ds == Dataset(x = [3,1,2,2,missing,3,3], y = [1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], z = [11,15,7,-11,12,0,0], q=[5.0, missing, -1,-1,missing,5,5], q2=[5.0, missing, -1,-1,missing,5,5])

    # original data set
    ds = Dataset(x = [3,1,2,2,missing,3,3], y = [1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], z = [11,15,7,-11,12,0,0])
    sds1 = dropmissing(ds, 2, view = true)
    sds2 = view(ds, [2,1,1,3,5,6,7,4,4], [:z, :x])

    modify!(groupby(sds2, :x), :z=>IMD.maximum=>:q)
    @test ds == Dataset([Union{Missing, Int64}[3, 1, 2, 2, missing, 3, 3], Union{Missing, Float64}[1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], Union{Missing, Int64}[11, 15, 7, -11, 12, 0, 0], Union{Missing, Int64}[11, 15, 7, 7, 12, 11, 11]], ["x", "y", "z", "q"])
    modify!(gatherby(sds2, :x), :z=>IMD.maximum=>:q2)
    @test ds == Dataset([Union{Missing, Int64}[3, 1, 2, 2, missing, 3, 3], Union{Missing, Float64}[1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], Union{Missing, Int64}[11, 15, 7, -11, 12, 0, 0], Union{Missing, Int64}[11, 15, 7, 7, 12, 11, 11], [11, 15, 7, 7, 12, 11, 11]], ["x", "y", "z", "q","q2"])

    # original data set
    ds = Dataset(x = [3,1,2,2,missing,3,3], y = [1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], z = [11,15,7,-11,12,0,0])
    sds1 = dropmissing(ds, 2, view = true)
    sds2 = view(ds, [2,1,1,3,5,6,7,4,4], [:z, :x])

    modify!(groupby(sds1, :x), (2,3)=>cor)
    @test ds == Dataset([Union{Missing, Int64}[3, 1, 2, 2, missing, 3, 3], Union{Missing, Float64}[1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], Union{Missing, Int64}[11, 15, 7, -11, 12, 0, 0], Union{Missing, Float64}[-0.9690582663799521, missing, 1.0, 1.0, missing, -0.9690582663799521, -0.9690582663799521]], ["x", "y", "z", "cor_y_z"])
    modify!(gatherby(sds1, :x), (2,3)=>cor)
    @test ds == Dataset([Union{Missing, Int64}[3, 1, 2, 2, missing, 3, 3], Union{Missing, Float64}[1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], Union{Missing, Int64}[11, 15, 7, -11, 12, 0, 0], Union{Missing, Float64}[-0.9690582663799521, missing, 1.0, 1.0, missing, -0.9690582663799521, -0.9690582663799521]], ["x", "y", "z", "cor_y_z"])
    # sds2 is not well defined anymore, since the original data has been changed
    sds2 = view(ds, [2,1,1,3,5,6,7,4,4], [:z, :x])

    modify!(groupby(sds2, :x), (:x, :z)=>((x,y)->length(x)/IMD.maximum(y))=>:q, [:z, :q]=>byrow(maximum)=>:q2)
    @test ds == Dataset([Union{Missing, Int64}[3, 1, 2, 2, missing, 3, 3], Union{Missing, Float64}[1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], Union{Missing, Int64}[11, 15, 7, -11, 12, 0, 0], Union{Missing, Float64}[-0.9690582663799521, missing, 1.0, 1.0, missing, -0.9690582663799521, -0.9690582663799521], Union{Missing, Float64}[0.36363636363636365, 0.06666666666666667, 0.42857142857142855, 0.42857142857142855, 0.08333333333333333, 0.36363636363636365, 0.36363636363636365], Union{Missing, Float64}[11.0, 15.0, 7.0, 0.42857142857142855, 12.0, 0.36363636363636365, 0.36363636363636365]], ["x", "y", "z", "cor_y_z", "q", "q2"])

    @test IMD.index(view(ds,[2,1,1,3,5,6,7,4,4],[:z,:x,:q,:q2])) == IMD.index(sds2)

    ds = Dataset(x = [3,1,2,2,missing,3,3], y = [1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], z = [11,15,7,-11,12,0,0])
    groupby!(ds, :y)
    sds1 = dropmissing(ds, 2, view = true)
    sds2 = view(ds, [2,1,1,3,5,6,7,4,4], [:z, :x])

    modify!(groupby(sds1, :x), (2,3)=>cor)
    @test ds == Dataset([Union{Missing, Int64}[2, 2, 3, 3, 3, 1, missing], Union{Missing, Float64}[-3.0, -1.0, 1.1, 4.0, 5.0, missing, missing], Union{Missing, Int64}[-11, 7, 11, 0, 0, 15, 12], Union{Missing, Float64}[1.0, 1.0, -0.9690582663799521, -0.9690582663799521, -0.9690582663799521, missing, missing]], ["x", "y", "z", "cor_y_z"])
    modify!(gatherby(sds1, :x), (2,3)=>cor)
    @test ds == Dataset([Union{Missing, Int64}[2, 2, 3, 3, 3, 1, missing], Union{Missing, Float64}[-3.0, -1.0, 1.1, 4.0, 5.0, missing, missing], Union{Missing, Int64}[-11, 7, 11, 0, 0, 15, 12], Union{Missing, Float64}[1.0, 1.0, -0.9690582663799521, -0.9690582663799521, -0.9690582663799521, missing, missing]], ["x", "y", "z", "cor_y_z"])
    # sds2 is not well defined anymore, since the original data has been changed
    sds2 = view(ds, [2,1,1,3,5,6,7,4,4], [:z, :x])

    modify!(groupby(sds2, :x), (:x, :z)=>((x,y)->length(x)/IMD.maximum(y))=>:q, [:z, :q]=>byrow(maximum)=>:q2)
    @test ds == Dataset([Union{Missing, Int64}[2, 2, 3, 3, 3, 1, missing], Union{Missing, Float64}[-3.0, -1.0, 1.1, 4.0, 5.0, missing, missing], Union{Missing, Int64}[-11, 7, 11, 0, 0, 15, 12], Union{Missing, Float64}[1.0, 1.0, -0.9690582663799521, -0.9690582663799521, -0.9690582663799521, missing, missing], Union{Missing, Float64}[0.42857142857142855, 0.42857142857142855, 0.36363636363636365, 0.36363636363636365, 0.36363636363636365, 0.06666666666666667, 0.08333333333333333], Union{Missing, Float64}[0.42857142857142855, 7.0, 11.0, 0.36363636363636365, 0.36363636363636365, 15.0, 12.0]], ["x", "y", "z", "cor_y_z", "q", "q2"])

    @test IMD.index(view(ds,[2,1,1,3,5,6,7,4,4],[:z,:x,:q,:q2])) == IMD.index(sds2)

    ds = Dataset(x = [3,1,2,2,missing,3,3], y = [1.1, missing, -1.0, -3.0, missing, 4.0, 5.0], z = [11,15,7,-11,12,0,0])
    sds2 = view(ds, [2,1,1,3,5,6,7,4,4], 1:2)
    @test combine(gatherby(sds2, 1), :y=>IMD.sum) == Dataset(x=[1,3,2, missing],  sum_y=[missing, 11.2,-7.0,missing])
    @test combine(gatherby(sds2, 1), :y=>(x->IMD.sum(x))=>:sum_y) == Dataset(x=[1,3,2, missing],  sum_y=[missing, 11.2,-7.0,missing])
    @test combine(groupby(sds2, 1), :y=>IMD.sum) == Dataset(x=[1,2,3, missing],  sum_y=[missing, -7.0, 11.2,missing])
    @test combine(gatherby(sds2, 1, isgathered = true), :y=>IMD.sum) == Dataset(x = [1,3,2,missing,3,2], sum_y=[missing, 2.2, -1,missing, 9,-6])

    ds = Dataset(x = [1,2,1,2,3], y1 = Union{Int8, Missing}[1,2,missing,4,missing], y2 = Union{Int32, Missing}[1,2,3,4,missing], y3=Union{Int16, Missing}[100,20,3000,4,missing], y4=Float16.(rand(5)), y5=rand(BigFloat, 5), y6=[missing, missing, missing, missing, missing])
    sds = view(ds, [1,2,3,4,5], [1,2,3,4,5,6,7])

    @test combine(gatherby(sds, 1), 2:4 .=>Ref([IMD.sum, mean, IMD.maximum, IMD.minimum, IMD.n, IMD.nmissing])) == Dataset([Union{Missing, Int64}[1, 2, 3], Union{Missing, Int64}[1, 6, missing], Union{Missing, Float64}[1.0, 3.0, missing], Union{Missing, Int8}[1, 4, missing], Union{Missing, Int8}[1, 2, missing], Union{Missing, Int64}[1, 2, 0], Union{Missing, Int64}[1, 0, 1], Union{Missing, Int64}[4, 6, missing], Union{Missing, Float64}[2.0, 3.0, missing], Union{Missing, Int32}[3, 4, missing], Union{Missing, Int32}[1, 2, missing], Union{Missing, Int64}[2, 2, 0], Union{Missing, Int64}[0, 0, 1], Union{Missing, Int64}[3100, 24, missing], Union{Missing, Float64}[1550.0, 12.0, missing], Union{Missing, Int16}[3000, 20, missing], Union{Missing, Int16}[100, 4, missing], Union{Missing, Int64}[2, 2, 0], Union{Missing, Int64}[0, 0, 1]], ["x", "sum_y1", "mean_y1", "maximum_y1", "minimum_y1", "n_y1", "nmissing_y1", "sum_y2", "mean_y2", "maximum_y2", "minimum_y2", "n_y2", "nmissing_y2", "sum_y3", "mean_y3", "maximum_y3", "minimum_y3", "n_y3", "nmissing_y3"])

    @test combine(gatherby(sds,1), :y6=>IMD.sum=>:mm) == combine(gatherby(sds,1), :y6=>IMD.minimum=>:mm) == combine(gatherby(sds,1), :y6=>var=>:mm) == combine(gatherby(sds,1), :y6=>(x->IMD.sum(x))=>:mm) == combine(gatherby(sds,1), :y6=>(x->var(x))=>:mm) == combine(gatherby(sds,1), :y6=>(x->IMD.minimum(x))=>:mm) == Dataset(x = [1,2,3] , mm = [missing, missing, missing])

    var1(x) = var(x)
    std1(x) = std(x)
    median1(x) = median(x)
    c1 =combine(gatherby(sds, 1), 2:4 .=>Ref([var, std, median]))
    c2 =  combine(gatherby(copy(sds), 1), 2:4 .=> Ref([var1, std1, median1]))
    @test byrow(compare(c1, c2, on = names(c1) .=> names(c2)) , all)|>all

    c3 = combine(gatherby(sds,1), :y4=>IMD.sum)
    @test eltype(c3.sum_y4) == Union{Missing, Float16}
    c3 = combine(gatherby(sds,1), :y5=>IMD.sum)
    @test eltype(c3.sum_y5) == Union{Missing, BigFloat}

    ds = Dataset(rand(1:10, 300_000,3), :auto)
    insertcols!(ds, 1, :g => repeat(1:150_000, inner = 2))
    map!(ds, x->rand()<.7 ? missing : x, r"x")
    sds = view(ds, :, :)

    @test combine(gatherby(sds, 1), r"x"=>IMD.sum) == combine(groupby(sds, 1), r"x"=>IMD.sum)
    @test combine(gatherby(sds, 1), r"x"=>IMD.maximum) == combine(groupby(sds, 1), r"x"=>IMD.maximum)
    @test combine(gatherby(sds, 1), r"x"=>IMD.minimum) == combine(groupby(sds, 1), r"x"=>IMD.minimum)
    @test combine(gatherby(sds, 1), r"x"=>var) == combine(groupby(sds, 1), r"x"=>var)

end


text:
module FluxPrune

export prune, iterativeprune,
       LevelPrune, ThresholdPrune,
       ChannelPrune

using MaskedArrays
using MaskedArrays: mask
using Functors: isleaf, functor
using LinearAlgebra
using Flux

include("functor.jl")
include("unstructured.jl")
include("structured.jl")
include("prune.jl")

end

ext: jl


text:
## Determine the number of flights from each airport, include a list of any airports not used.
NMR.clean_split_dir()

jobid = 1
NMR.runjob(NMR.Nmr(jobid, ["AComp_Passenger_data.csv"], NMR.mapper_parserecordacomp, "acomp.csv"))
NMR.runjob(NMR.Nmr(jobid, ["acomp.csv"], NMR.reducer_numberofflights, "acomp_flights_reduced.csv"))
c = NMR.runcombiner("acomp_flights_reduced.csv", NMR.combiner_parsejson)

# store output
header = "airportid,flights"
output = [string(JSON.parse(l)[1]) * "," * string(JSON.parse(l)[2]) for l in c]
NMR.save_job_output(jobid, header, output)


text:
using PyPlot, Seismic, ElasticWave

nz = 200; nx = 200; # model size
ext= 20; # thickness of PML layers
iflag = 1; #when iflag=1, four side obsorbing boundary, iflag=2, free surface
dx = 5. ; dz = 5. ; dt = 5e-4; f0=30.; tmax=0.8; # recording time and dominant frequency of source wavelet
rho=2.5e-4*ones(nz,nx);
# vp=2500.*ones(nz,nx);  vs=2500./sqrt(3)*ones(nz,nx);
# vp[20:end,:]=3000.; vs[20:end,:] = 3000./sqrt(3);
vp=3000.*ones(nz,nx);  vs=3000./sqrt(3)*ones(nz,nx);
vp[150:end,:]=4000.;  vs[150:end,:]=4000./sqrt(3); #physical model
fidMtx = CreateFidMtx(nz, nx, ext, iflag, vp, vs, rho, dz, dx, dt, f0); #finite difference matrix

# single source
isz = [100, 100]; isx = [85, 115]; ot = [0.0, 0.1]; #source location and starting time
flags = vec([false false true true false]); # add source to normal stress, simulate explosive sources
srcs = InitMultiSources(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);

woz = 51; wox=51; wnz = 100; wnx = 100;
wsc = Srcs2Wsc(woz, wox, wnz, wnx, srcs);
wnt = 500; tmp = zeros(wnt,2);
tmp[1:346,1] = wsc.p[50,35,:];
tmp[1:346,2] = wsc.p[50,65,:];

fig = figure("1", figsize=(12,6)); s=25;
subplot(1,2,1); Seismic.SeisPlot(fignum="1", vp, cmap="jet", vmax=4000, vmin=3000, dx=5, dy=5, xlabel="X (m)", ylabel="Y (m)", ticksize=s, labelsize=s);
ax = gca(); ax[:set_xticks]([200, 400, 600, 800]); ax[:set_yticks]([200, 400, 600, 800]);
ax[:text](-100,60, "a)", fontsize=s, fontweight="bold")
subplot(1,2,2); Seismic.SeisPlot(fignum="1", tmp, style="wiggles", xcur=1.0, dy=0.0005, xlabel="Trace number", ylabel="Time (s)", ticksize=s, labelsize=s);
ax = gca(); ax[:set_xticks]([0, 1]); ax[:set_yticks]([0.05, 0.10, 0.15, 0.20]);
ax[:text](-1.3,0.015, "b)", fontsize=s, fontweight="bold")
tight_layout()

irx = collect(1:4:nx); irz = ones(Int64,length(irx));
tmp = collect(1:4:nz); irz = vcat(irz, tmp); irx = vcat(irx, 5*ones(Int64,length(tmp)));
shot = MultiStepForward(nz, nx, irz, irx, wsc, fidMtx,  tmax=tmax);
tmpx = shot.Vx; tmpz = shot.Vz;
shot.Vz = SeisAddNoise(shot.Vz, 0.5, L=5);
shot.Vx = SeisAddNoise(shot.Vx, 0.5, L=5);

fig = figure("2", figsize=(12,12)); s=15;
subplot(2,2,1); SeisPlot(fignum="2", tmpx, cmap="gray", pclip=95, dy=0.0005, ylabel="Time (s)", ticksize=s, labelsize=s)
ax = gca(); ax[:set_xticks]([20, 40, 60, 80]); ax[:set_yticks]([0.2, 0.4, 0.6]);
ax[:text](-8,0.04, "a)", fontsize=20, fontweight="bold")

subplot(2,2,2); SeisPlot(fignum="2", tmpz, cmap="gray", pclip=95, dy=0.0005, ticksize=s, labelsize=s)
ax = gca(); ax[:set_xticks]([20, 40, 60, 80]); ax[:set_yticks]([0.2, 0.4, 0.6]);
ax[:text](-8,0.04, "b)", fontsize=20, fontweight="bold")

subplot(2,2,3); SeisPlot(fignum="2", shot.Vx, cmap="gray", pclip=95, dy=0.0005, ylabel="Time (s)", xlabel="Trace number", ticksize=s, labelsize=s)
ax = gca(); ax[:set_xticks]([20, 40, 60, 80]); ax[:set_yticks]([0.2, 0.4, 0.6]);
ax[:text](-8,0.04, "c)", fontsize=20, fontweight="bold")

subplot(2,2,4); SeisPlot(fignum="2", shot.Vz, cmap="gray", pclip=95, dy=0.0005, xlabel="Trace number", ticksize=s, labelsize=s)
ax = gca(); ax[:set_xticks]([20, 40, 60, 80]); ax[:set_yticks]([0.2, 0.4, 0.6]);
ax[:text](-8,0.04, "d)", fontsize=20, fontweight="bold")

# wrand = RandWsc(woz, wox, wnz, wnx, ext, iflag, dt, wnt, f0);
# lambda = power_wsc(nz, nx, irz, irx, fidMtx, wrand, tmax; niter=30);

lambda = 1415.; mu=0.1;
(J, winv) = FISTA(woz, wox, wnz, wnx, wnt, shot, fidMtx, mu, lambda; niter=50);
lambda = 1415.; mu=0.01;
(J1, winv1) = FISTA(woz, wox, wnz, wnx, wnt, shot, fidMtx, mu, lambda; niter=50);
wadj = MultiStepAdjoint(1, 1, 200, 200, wnt, shot, fidMtx);
a = EngDisWsc(wadj);
p = EngDisWsc(winv); p1 = zeros(nz,nx); p1[woz:woz+wnz-1, wox:wox+wnx-1] = p[:,:];
v = maximum(p1)
SeisPlot(p1, vmax=v, vmin=0, cmap="jet", pclip=100)
v1 = maximum(a)
SeisPlot(a, vmax=v1, vmin=0, cmap="jet", pclip=100)

max_stars_repo_name: WneleiGao/ElasticWave


text:
using JuMP, ReSHOP, EMP
using Test

@testset "simple mopec" begin

ATmat = [1., -1., -1.]
s = [.9, .1, 0]
b = [0, 5, 3]

mopec = EMPmaster()

ag = MathPrgm(mopec)
mkt = MathPrgm(mopec)

EquilibriumProblem(mopec, [ag, mkt])

n = 3

@variable(mkt, y >= 0)
@variable(ag, x[1:n] >= 0)
@variable(mkt, p[1:n] >= 0)

JuMP.fix(p[2], 1.; force=true)
JuMP.set_start_value.(x, 1.)

@NLobjective(ag, Max, sum(s[i] * log(x[i]) for i=1:n))
@constraint(ag, sum(p[i]*x[i] for i=1:n) <= sum(p[i]*b[i] for i=1:n) )

@constraint(mkt, [b[i] + ATmat[i]*y - x[i] for i=1:n] ⟂ p)
@constraint(mkt, sum(-ATmat[i]*p[i] for i=1:n) ⟂ y)

solveEMP(mopec)

@test isapprox(value.(x)[:], [3,2,0], atol=1e-4)
@test isapprox(value.(p)[:], [6,1,5], atol=1e-4)
@test isapprox(value(y), 3, atol=1e-4)

#@test isapprox(getobjectivevalue(m),5326.851310161077, atol=1e-5)

end


text:
#
#	Quantum State struct
#
struct QState{T<:Number}
	Coef::AbstractVector{T}
	QBasis::Union{<:AbstractQBasis}
end

#
#	basic function on Quantum States
#

# shits basis of periodic state by one to the right 11000 -> 10001 useful to obtain entanglement
# at different cut (for instance intra unit-cell cut instead of inter for a 2 site unite-cell system
function shiftStateBasis(state::QState)
	new_coef = zeros(eltype(state.Coef),size(state.Coef))
	@inbounds for x in state.QBasis.OccRep
		add = false
		if getBitPos(x[1])[1] == 1
			add = true
		end

		new_n = x[1]>>1
		if add == true
			new_n += 2^(state.QBasis.N-1)
		end

		new_idx = fermionicHashfunc(new_n)+1

		new_coef[new_idx] = state.Coef[x[2]]

	end

	state.Coef[:] .= new_coef[:]
	return state
end

function getCijMatrix(Coef::AbstractVector{<:Number}, N::Int, Np::Int; Basis::Union{AbstractQBasis, Nothing} = nothing)
	CijMatrix = Array{ComplexF64,2}(undef, N, N)

	for i in 1:N
		cr_i = nothing
		try
			cr_i = loadFermionicCr(i, N, Np)
		catch
			cr_i = fermionicCrMatrix(i, N, Np; preQBasis=Basis)
			saveFermionicCr(cr_i, i, N, Np)
		end
		for j in i:N
			cr_j = nothing
			try
				cr_j = loadFermionicCr(j, N, Np)
			catch
				cr_j = fermionicCrMatrix(j, N, Np; preQBasis=Basis)
				saveFermionicCr(cr_j, j, N, Np)
			end

			#cr_j = fermionicCrMatrix(j, N, Np; preQBasis = Basis)
			C = cr_i*transpose(cr_j)
			CijMatrix[i,j] = transpose(conj(Coef))*C*Coef
			if i != j
				CijMatrix[j,i] = conj(CijMatrix[i,j])
			end
		end
	end

	return CijMatrix
end

function getCijMatrix(Coef::AbstractVector{<:Number}, N::Int, Np::Int, symSec::Vector{Int}; Basis::Union{AbstractQBasis, Nothing} = nothing)
	CijMatrix = Array{ComplexF64,2}(undef, N, N)
	Coef = Coef[symSec]

	for i in 1:N
		cr_i = nothing
		try
			cr_i = loadFermionicCr(i, N, Np)
		catch
			cr_i = fermionicCrMatrix(i, N, Np; preQBasis=Basis)
			saveFermionicCr(cr_i, i, N, Np)
		end

		cr_i = cr_[symSec,:]

		for j in i:N
			cr_j = nothing
			try
				cr_j = loadFermionicCr(j, N, Np)
			catch
				cr_j = fermionicCrMatrix(j, N, Np; preQBasis=Basis)
				saveFermionicCr(cr_j, j, N, Np)
			end

			cr_j = cr_j[:,symSec]
			#cr_j = fermionicCrMatrix(j, N, Np; preQBasis = Basis)
			C = cr_i*transpose(cr_j)
			CijMatrix[i,j] = transpose(conj(Coef))*C*Coef
			if i != j
				CijMatrix[j,i] = conj(CijMatrix[i,j])
			end
		end
	end

	return CijMatrix
end

function getCijMatrix(Coef::AbstractVector{<:Number}, i_idx::Vector{Int}, j_idx::Vector{Int}, N::Int, Np::Int; Basis::Union{AbstractQBasis, Nothing} = nothing)
	CijMatrix = Array{ComplexF64,2}(undef, size(i_idx)[1], size(j_idx)[1])

	for (n,i) in enumerate(i_idx)
		cr_i = nothing
		try
			cr_i = loadFermionicCr(i, N, Np)
		catch
			cr_i = fermionicCrMatrix(i, N, Np; preQBasis=Basis)
			saveFermionicCr(cr_i, i, N, Np)
		end
		for (m,j) in enumerate(j_idx)
			cr_j = nothing
			try
				cr_j = loadFermionicCr(j, N, Np)
			catch
				cr_j = fermionicCrMatrix(j, N, Np; preQBasis=Basis)
				saveFermionicCr(cr_j, j, N, Np)
			end

			#cr_j = fermionicCrMatrix(j, N, Np; preQBasis = Basis)
			C = cr_i*transpose(cr_j)
			CijMatrix[n,m] = transpose(conj(Coef))*C*Coef
			if i != j
				CijMatrix[m,n] = conj(CijMatrix[n,m])
			end
		end
	end

	return CijMatrix
end

getCijMatrix(State::QState) = getCijMatrix(State.Coef, State.QBasis.N, State.QBasis.Np-1)
getCijMatrix(State::QState, i_idx::Vector{Int}, j_idx::Vector{Int}) = getCijMatrix(State.Coef, i_idx, j_idx, State.QBasis.N, State.QBasis.Np-1)
getCijMatrix(State::QState, symSec::Vector{Int}) = getCijMatrix(State.Coef, State.QBasis.N, State.QBasis.Np-1, symSec)

function projectOntoStates(States::AbstractVector{<:Any}, CijMatrix::AbstractArray{<:Number})

	occ = 0
	for x in States
		occ += (transpose(x)*CijMatrix*conj(x))[1,1]
	end

	return occ

end

function entanglementHalfCut(Coef::AbstractVector{<:Number}, N::Int, Np::Int; Basis::Union{AbstractQBasis, Nothing} = nothing)
	if Basis == nothing

		Basis = QFermionicBasis(N, Np)
	end

	dim = Int(sum(2.0.^[i for i = 0:Int(N/2-1)]))+1
	CoefMatrix = spzeros(ComplexF64,dim,dim)

	for idx in Basis.OccRep

		nleft = idx[1] >> Int(N/2)
		nright = idx[1] - (nleft << Int(N/2))+1
		nleft += 1
		#@show idx[2], idx[1], nleft, nright, Coef[idx[2]]
		#if I want to do symsec I have to count bits here

		CoefMatrix[nleft,nright] = Coef[idx[2]]

	end

	#@show CoefMatrix

	quality = 0.0
	k = 10
	svdHandle = []
	svdHandle2 = []
	s = []
	while quality < 0.99
		try
			svdHandle = svds(CoefMatrix, nsv = k)
			svdHandle2 = svdsolve(CoefMatrix, k)
		catch
			# replace by more meaningful analysis of the issue
			# add logging for the package
			@show "something went wrong"
		end
		quality = sum(svdHandle[1].S.^2.0)
		k += 10
		s = svdHandle[1].S[:]

	end

	@show s
	@show svdHandle2[1]
	@show max(real.((CoefMatrix - sparse(svdHandle[1].U*diagm(0 => s)*transpose(conj(svdHandle[1].V)))).nzval)...)

	@show max(imag.((CoefMatrix - sparse(svdHandle[1].U*diagm(0 => s)*transpose(conj(svdHandle[1].V)))).nzval)...)
	return s, svdHandle[1].U, svdHandle[1].V, quality
end

entanglementHalfCut(State::QState) = entanglementHalfCut(State.Coef, State.QBasis.N, State.QBasis.Np, Basis = State.QBasis)

max_stars_repo_name: Gawatz/exactQmodule

max_issues_repo_name: Gawatz/exactQmodule

max_forks_repo_name: Gawatz/exactQmodule


text:
using ThArrays
using Test

@testset "Tensor Creation" begin

    @testset "Creation with Array" begin
        ary = rand(2, 3)
        @test Tensor(ary) == Tensor(ary)
        @test convert(Array, Tensor(ary)) == ary
    end

    @testset "Creation with Array (sharing data)" begin
        ary = rand(2, 3)
        t = Tensor(ary)
        ThArrays.ThC.sin!(t)
        @test t == Tensor(ary)
    end

    @testset "Creation with Array (copying data)" begin
        ary = rand(2, 3)
        t = Tensor(ary, detach=true)
        ThArrays.ThC.sin!(t)
        @test isapprox(convert(Array, t), sin.(ary), atol=0.001)
    end

    @testset "Create with Number (0-dim Tensor)" begin
        t = Tensor(1.0)
        @test size(t) == ()
        @test length(t) == 1
        @test t[] == 1.0
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

# test Pyr5 elasticity with point load

fn = dirname(@__DIR__) * "/geometry/3d_pyr/Pyr5.med"
mesh = aster_read_mesh(fn)
element_sets = join(keys(mesh.element_sets), ", ")
@debug("element sets: $element_sets")

element1 = create_elements(mesh,"Pyr5")
baseQuad = create_elements(mesh,"baseQuad")
update!(element1, "youngs modulus", 288.0)
update!(element1, "poissons ratio", 1/3)

tip_node_id = first(mesh.node_sets[:tipPoint])
tipPoint = Element(Poi1, [tip_node_id])
update!(tipPoint, "geometry", mesh.nodes)
update!(tipPoint, "displacement traction force 1",  5.0)
update!(tipPoint, "displacement traction force 2", -7.0)
update!(tipPoint, "displacement traction force 3",  3.0)

problem = Problem(Elasticity, "solve continuum block", 3)
problem.properties.finite_strain = false
add_elements!(problem, element1, tipPoint)

update!(first(baseQuad), "displacement 1", 0.0)
update!(first(baseQuad), "displacement 2", 0.0)
update!(first(baseQuad), "displacement 3", 0.0)
bc = Problem(Dirichlet, "Boundary conditions", 3, "displacement")
add_elements!(bc, baseQuad)

analysis = Analysis(Linear)
add_problems!(analysis, problem, bc)
run!(analysis)

xi, time = (0.0, 0.0, 1.0), 0.0
u = first(element1)("displacement", xi, time)
# Code_Aster Result in verification/2017-05-27-pyramids/Pyr5_displacement.txt
u_expected = [6.9444444444427100E-02,-9.7222222222197952E-02,1.0416666666679683E-02]
@debug("displacement at tip", u, u_expected)
@test isapprox(u, u_expected)


text:
using GadgetronExamples
using Test

@testset "GadgetronExamples.jl" begin
    # Write your tests here.
end


text:
ForwardDiff_gradient(x...) = ForwardDiff.gradient(x...)

function kr(psi, alpha, N)
	if psi < 0
		m = (N - 1) / N
		denom = 1 + abs(alpha * psi) ^ N
		numer = 1 - abs(alpha * psi) ^ (N - 1) * denom ^ (-m)
		return numer ^ 2 / denom ^ (m / 2)
	else
		return one(psi)
	end
end

kr(x::AbstractArray) = kr(x[1], x[2], x[3])

function Calculus.differentiate(x::Calculus.SymbolParameter{:kr}, args, wrt)
	chain_part = map(x->Calculus.simplify(Calculus.differentiate(x, wrt)), args)
	if chain_part == [0, 0, 0]
		return :(0)
	else
		return :(sum(ForwardDiff_gradient(kr, StaticArrays.SA[$(args...)]) .* StaticArrays.SA[$(chain_part...)]))
	end
end

function hm(x, y)
	return 2 / (1 / x + 1 / y)
end

hm(x::AbstractArray) = hm(x[1], x[2])

function Calculus.differentiate(x::Calculus.SymbolParameter{:hm}, args, wrt)
	chain_part = map(x->Calculus.simplify(Calculus.differentiate(x, wrt)), args)
	if chain_part == [0, 0]
		return :(0)
	else
		return :(sum(ForwardDiff_gradient(hm, StaticArrays.SA[$(args...)]) .* StaticArrays.SA[$(chain_part...)]))
	end
end

function Calculus.differentiate(x::Calculus.SymbolParameter{:abs}, args, wrt)
	if length(args) > 1
		error("Too many arguments passed to abs()")
	end
	arg = args[1]
	return :(ifelse($arg > 0, 1, -1) * $(Calculus.differentiate(arg, wrt)))
end

function effective_saturation(alpha::Number, psi::Number, N::Number)
	m = (N - 1) / N
	if psi < 0
		return (1 + abs(alpha * psi) ^ N) ^ (-m)
	else
		return 1
	end
end

@NonlinearEquations.equations exclude=(coords, dirichletnodes, neighbors, areasoverlengths) function richards(psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs, specificstorage, volumes)
	NonlinearEquations.setnumequations(length(psi))
	dirichletnodes_set = Set(dirichletnodes)
	for i = 1:length(psi)
		if i in dirichletnodes_set
			NonlinearEquations.addterm(i, psi[i] - dirichletpsis[i])
		else
			NonlinearEquations.addterm(i, Qs[i] / (specificstorage[i] * volumes[i]))
		end
	end
	for (i, (node_a, node_b)) in enumerate(neighbors)
		for (node1, node2) in [(node_a, node_b), (node_b, node_a)]
			if !(node1 in dirichletnodes_set) && !(node2 in dirichletnodes_set)
				NonlinearEquations.addterm(node1, hm(kr(psi[node1], alphas[i], Ns[i]), kr(psi[node2], alphas[i], Ns[i])) * Ks[i] * (psi[node2] + coords[3, node2] - psi[node1] - coords[3, node1]) * areasoverlengths[i] / (specificstorage[node1] * volumes[node1]))
			elseif !(node1 in dirichletnodes_set) && node2 in dirichletnodes_set
				NonlinearEquations.addterm(node1, hm(kr(psi[node1], alphas[i], Ns[i]), kr(dirichletpsis[node2], alphas[i], Ns[i])) * Ks[i] * (dirichletpsis[node2] + coords[3, node2] - psi[node1] - coords[3, node1]) * areasoverlengths[i] / (specificstorage[node1] * volumes[node1]))
			end
		end
	end
end

function richards_steadystate(psi0, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs; kwargs...)
	args = (psi0, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs, ones(length(Qs)), ones(length(Qs)))
	function residuals!(residuals, psi)
		myargs = (psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs, ones(length(Qs)), ones(length(Qs)))
		copy!(residuals, DPFEHM.richards_residuals(myargs...))
	end
	function jacobian!(J, psi)
		myargs = (psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs, ones(length(Qs)), ones(length(Qs)))
		DPFEHM.richards_psi!(J, myargs...)
	end
	J0 = DPFEHM.richards_psi(args...)
	residuals0 = DPFEHM.richards_residuals(args...)
	df = NLsolve.OnceDifferentiable(residuals!, jacobian!, psi0, residuals0, J0)
	soln = NLsolve.nlsolve(df, psi0; kwargs...)
	if !NLsolve.converged(soln)
		display(soln)
		error("solution did not converge")
	end
	return soln.zero
end

function ChainRulesCore.rrule(::typeof(richards_steadystate), psi0, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs; kwargs...)
	psi = richards_steadystate(psi0, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs; kwargs...)
	function pullback(delta)
		args = (psi, Ks, neighbors, areasoverlengths, dirichletnodes, dirichletpsis, coords, alphas, Ns, Qs, ones(length(Qs)), ones(length(Qs)))
		lambda = richards_psi(args...)' \ delta
		req_Ks = richards_Ks(args...)
		req_dirichletpsis = richards_dirichletpsis(args...)
		req_alphas = richards_alphas(args...)
		req_Ns = richards_Ns(args...)
		req_Qs = richards_Qs(args...)
		return (ChainRulesCore.NoTangent(),#step function
				@ChainRulesCore.thunk(zeros(length(psi0))),#psi0 -- should be insensitive to psi0
				@ChainRulesCore.thunk(-(req_Ks' * lambda)),#Ks
				@ChainRulesCore.thunk(ChainRulesCore.NoTangent()),#neighbors
				@ChainRulesCore.thunk(ChainRulesCore.NoTangent()),#areasoverlengths
				@ChainRulesCore.thunk(ChainRulesCore.NoTangent()),#dirichletnodes
				@ChainRulesCore.thunk(-(req_dirichletpsis' * lambda)),#dirichletpsis
				@ChainRulesCore.thunk(ChainRulesCore.NoTangent()),#coords
				@ChainRulesCore.thunk(-(req_alphas' * lambda)),#alphas
				@ChainRulesCore.thunk(-(req_Ns' * lambda)),#Ns
				@ChainRulesCore.thunk(-(req_Qs' * lambda)))#Qs
	end
	return psi, pullback
end


text:
# --------------------------------------------------------------------------
# ACE.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

@testset "Fourier  Basis" begin

@info("Testset FourierBasis")

##

using StaticArrays, Test
using LinearAlgebra
using JuLIP.Testing: print_tf

using ACE
using JuLIP: evaluate!, evaluate, evaluate_d!
using ACE: alloc_B, alloc_dB
using ACE.Bonds: CylindricalCoordinateSystem, cylindrical,
                  cartesian,
                  FourierBasis,
                  CylindricalCoordinates

##

deg = 10
fB = FourierBasis(10, Float64)

##

@info("test cyl_l2i ∘ cyl_i2l = id")
is = sort([ ACE.Bonds.cyl_l2i(k) for k = -deg:deg ])
println(@test is == 1:(2*deg+1))
for i = 1:(2*deg+1)
   k = ACE.Bonds.cyl_i2l(i)
   print_tf(@test ACE.Bonds.cyl_l2i(k) == i)
end

##
# choose a few random angles, then confirm that the fourier basis
# does the right thing against an explicit expression using a
# comprehension

@info("Test correctness of Fourier / Trigonometric Polynomial basis")
for ntest = 1:30
   θ = 2 * π * rand()
   cosθ, sinθ = cos(θ), sin(θ)
   c = CylindricalCoordinates(cosθ, sinθ, 0.0, 0.0)

   P = evaluate!(alloc_B(fB), nothing, fB, c)
   Px = zeros(ComplexF64, 2*deg+1)
   for l = -deg:deg
      Px[ACE.Bonds.cyl_l2i(l)] = exp(im * l * θ)
   end
   print_tf(@test P ≈ Px)
end

##

@info("Test correctness of gradient of Fourier Basis")
P = alloc_B(fB)
Ph = alloc_B(fB)
dP = alloc_dB(fB)

for ntest = 1:30
   θ = 2 * π * rand()
   cosθ, sinθ = cos(θ), sin(θ)
   c = CylindricalCoordinates(cosθ, sinθ, 0.0, 0.0)

   evaluate_d!(P, dP, nothing, fB, c)
   P1 = evaluate!(alloc_B(fB), nothing, fB, c)
   print_tf(@test P1 ≈ P)
   dP_dθ = [ dp[1] * (-sinθ) + dp[2] * cosθ  for dp in  dP ]

   errs = Float64[]
   for p = 2:10
      h = 0.1^p
      θh = θ + h
      ch = CylindricalCoordinates(cos(θh), sin(θh), 0.0, 0.0)
      evaluate!(Ph, nothing, fB, ch)
      dhP_dθ = (Ph - P) / h
      push!(errs, norm(dP_dθ - dhP_dθ, Inf))
   end
   print_tf(@test minimum(errs[2:end]) < 1e-3 * maximum(errs[1:3]))
end

##

end


text:
# high-level interface and shorthands

export current_window, current_buffer, current_tabpage, current_line, set_current, set_current_line, cursor, cursor!

current_window(c::NvimClient) = get_current_window(c)
current_buffer(c::NvimClient) = get_current_buffer(c)
current_tabpage(c::NvimClient) = get_current_tabpage(c)
current_line(c::NvimClient) = get_current_line(c)

set_current(o::Buffer) = set_current_buffer(o.client, o)
set_current(o::Window) = set_current_window(o.client, o)
set_current(o::Tabpage) = set_current_tabpage(o.client, o)

cursor(o::Window) = ((r,c) = get_cursor(o); (r, c + 1))
cursor!(o::Window, r::Integer, c::Integer) = set_cursor(o, (r, c - 1))

# array interface for buffer

# index relative end ( -1 is last element)
immutable EndRelIndex
    i::Int
end

vimindex(a::Integer) = a > 0 ? a - 1 : throw(BoundsError())
vimindex(a::EndRelIndex) = a.i < 0 ? a.i : throw(BoundsError())

# range where indices already are converted to api indicies
immutable CappedRange
    start::Int
    stop::Int
    CappedRange(a,b) = new(vimindex(a),vimindex(b))
end

Base.length(b::Buffer) = line_count(b)
Base.endof(b::Buffer) = EndRelIndex(-1)

import Base.-, Base.+

-(a::EndRelIndex, b::Integer) = EndRelIndex(a.i - b)
+(a::EndRelIndex, b::Integer) = EndRelIndex(a.i + b)
+(a::Integer, b::EndRelIndex) = EndRelIndex(a + b.i)

Base.colon(a::Integer, b::EndRelIndex) = CappedRange(a, b)
Base.colon(a::EndRelIndex, b::Union{Integer, EndRelIndex}) = CappedRange(a, b)

Base.getindex(b::Buffer, r::CappedRange) = get_line_slice(b, r.start, r.stop, true, true)
Base.getindex(b::Buffer, ::Colon) = b[1:end]
function Base.getindex(b::Buffer, i::Union{Integer, EndRelIndex})
    line = b[i:i]
    length(line) > 0 ? line[1] : ""
end
function Base.getindex{T<:Integer}(b::Buffer, r::UnitRange{T})
    if (r.start > r.stop)
        Array(ByteString, 0)
    else
        b[CappedRange(r.start, r.stop)]
    end
end

Base.setindex!(b::Buffer, lines::Array, r::CappedRange) =
    set_line_slice(b, r.start, r.stop, true, true, lines)
Base.setindex!(b::Buffer, lines::Array, ::Colon) = b[1:end] = lines
Base.setindex!(b::Buffer, s::AbstractString, i::Integer) = b[i:i] = [s]
Base.setindex!(b::Buffer, s::AbstractString, i::EndRelIndex) = b[i:i] = [s]
Base.setindex!(b::Buffer, s::AbstractString, r::CappedRange) = b[r] = [s]
function Base.setindex!{T<:Integer}(b::Buffer, s::AbstractString, r::UnitRange{T})
    ndests = r.stop - r.start + 1
    ndests <= 0 ? s : b[r] = fill(s, ndests)
end
function Base.setindex!{T<:Integer}(b::Buffer, lines::Array, r::UnitRange{T})
    b[CappedRange(r.start, r.stop)] = lines
end

Base.push!(b::Buffer, items...) = (insert(b, -1, [items...]); b)
Base.append!(b::Buffer, items) = (insert(b, -1, items); b)
Base.unshift!(b::Buffer, items...) = (set_line_slice(b, 0, 0, true, false, [items...]); b)

Base.deleteat!(b::Buffer, i::Integer) = deleteat!(b, i:i)
Base.deleteat!{T<:Integer}(b::Buffer, rng::UnitRange{T}) = (b[rng] = []; b)


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: braket

using AWS.Compat
using AWS.UUIDs
"""
    CancelQuantumTask()

Cancels the specified task.

# Required Parameters
- `clientToken`: The client token associated with the request.
- `quantumTaskArn`: The ARN of the task to cancel.

"""

cancel_quantum_task(clientToken, quantumTaskArn; aws_config::AWSConfig=global_aws_config()) = braket("PUT", "/quantum-task/$(quantumTaskArn)/cancel", Dict{String, Any}("clientToken"=>clientToken); aws_config=aws_config)
cancel_quantum_task(clientToken, quantumTaskArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("PUT", "/quantum-task/$(quantumTaskArn)/cancel", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>clientToken), args)); aws_config=aws_config)

"""
    CreateQuantumTask()

Creates a quantum task.

# Required Parameters
- `action`: The action associated with the task.
- `clientToken`: The client token associated with the request.
- `deviceArn`: The ARN of the device to run the task on.
- `outputS3Bucket`: The S3 bucket to store task result files in.
- `outputS3KeyPrefix`: The key prefix for the location in the S3 bucket to store task results in.
- `shots`: The number of shots to use for the task.

# Optional Parameters
- `deviceParameters`: The parameters for the device to run the task on.
"""

create_quantum_task(action, clientToken, deviceArn, outputS3Bucket, outputS3KeyPrefix, shots; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/quantum-task", Dict{String, Any}("action"=>action, "clientToken"=>clientToken, "deviceArn"=>deviceArn, "outputS3Bucket"=>outputS3Bucket, "outputS3KeyPrefix"=>outputS3KeyPrefix, "shots"=>shots); aws_config=aws_config)
create_quantum_task(action, clientToken, deviceArn, outputS3Bucket, outputS3KeyPrefix, shots, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/quantum-task", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("action"=>action, "clientToken"=>clientToken, "deviceArn"=>deviceArn, "outputS3Bucket"=>outputS3Bucket, "outputS3KeyPrefix"=>outputS3KeyPrefix, "shots"=>shots), args)); aws_config=aws_config)

"""
    GetDevice()

Retrieves the devices available in Amazon Braket.

# Required Parameters
- `deviceArn`: The ARN of the device to retrieve.

"""

get_device(deviceArn; aws_config::AWSConfig=global_aws_config()) = braket("GET", "/device/$(deviceArn)"; aws_config=aws_config)
get_device(deviceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("GET", "/device/$(deviceArn)", args; aws_config=aws_config)

"""
    GetQuantumTask()

Retrieves the specified quantum task.

# Required Parameters
- `quantumTaskArn`: the ARN of the task to retrieve.

"""

get_quantum_task(quantumTaskArn; aws_config::AWSConfig=global_aws_config()) = braket("GET", "/quantum-task/$(quantumTaskArn)"; aws_config=aws_config)
get_quantum_task(quantumTaskArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("GET", "/quantum-task/$(quantumTaskArn)", args; aws_config=aws_config)

"""
    SearchDevices()

Searches for devices using the specified filters.

# Required Parameters
- `filters`: The filter values to use to search for a device.

# Optional Parameters
- `maxResults`: The maximum number of results to return in the response.
- `nextToken`: A token used for pagination of results returned in the response. Use the token returned from the previous request continue results where the previous request ended.
"""

search_devices(filters; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/devices", Dict{String, Any}("filters"=>filters); aws_config=aws_config)
search_devices(filters, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/devices", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("filters"=>filters), args)); aws_config=aws_config)

"""
    SearchQuantumTasks()

Searches for tasks that match the specified filter values.

# Required Parameters
- `filters`: Array of SearchQuantumTasksFilter objects.

# Optional Parameters
- `maxResults`: Maximum number of results to return in the response.
- `nextToken`: A token used for pagination of results returned in the response. Use the token returned from the previous request continue results where the previous request ended.
"""

search_quantum_tasks(filters; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/quantum-tasks", Dict{String, Any}("filters"=>filters); aws_config=aws_config)
search_quantum_tasks(filters, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = braket("POST", "/quantum-tasks", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("filters"=>filters), args)); aws_config=aws_config)


text:
# Independent Component Analysis

"""
This type contains ICA model parameters: mean and component matrix ``W``.

**Note:** Each column of the component matrix ``W`` corresponds to an independent component.
"""
struct ICA{T<:Real} <: LinearDimensionalityReduction
    mean::Vector{T}   # mean vector, of length m (or empty to indicate zero mean)
    W::Matrix{T}      # component coefficient matrix, of size (m, k)
end
function show(io::IO, M::ICA)
    indim, outdim = size(M)
    print("ICA(indim=$indim, outdim=$outdim)")
end

"""
    size(M::ICA)

Returns a tuple with the input dimension, *i.e* the number of observed mixtures, and
the output dimension, *i.e* the number of independent components.
"""
size(M::ICA) = size(M.W)

"""
    mean(M::ICA)

Returns the mean vector.
"""
mean(M::ICA) = fullmean(size(M.W,1), M.mean)

"""
    predict(M::ICA, x)

Transform `x` to the output space to extract independent components,
as ``\\mathbf{W}^T (\\mathbf{x} - \\boldsymbol{\\mu})``, given the model `M`.
"""
predict(M::ICA, x::AbstractVecOrMat{<:Real}) = transpose(M.W) * centralize(x, M.mean)

#### core algorithm

"""
The abstract type for all `g` (derivative) functions.

Let `g` be an instance of such type, then `update!(g, U, E)` given

- `U = w'x`

returns updated in-place `U` and `E`, s.t.

- `g(w'x) --> U` and `E{g'(w'x)} --> E`

"""
abstract type ICAGDeriv end

"""
Derivative for ``(1/a_1)\\log\\cosh a_1 u``
"""
struct Tanh{T} <: ICAGDeriv
    a::T
end
function update!(f::Tanh{T}, U::AbstractMatrix{T}, E::AbstractVector{T}) where {T}
    n,k = size(U)
    a = f.a
    @inbounds for j in 1:k
        _s = zero(T)
        @fastmath for i in 1:n
            t = tanh(a * U[i,j])
            U[i,j] = t
            _s += a * (1 - t^2)
        end
        E[j] = _s / n
    end
end

"""
Derivative for ``-e^{\\frac{-u^2}{2}}``
"""
struct Gaus <: ICAGDeriv end
function update!(f::Gaus, U::AbstractMatrix{T}, E::AbstractVector{T}) where {T}
    n,k = size(U)
    @inbounds for j in 1:k
        _s = zero(T)
        for i in 1:n
            u = U[i,j]
            u2 = u^2
            e = exp(-u2/2)
            U[i,j] = u * e
            _s += (1 - u2) * e
        end
        E[j] = _s / n
    end
end

# Fast ICA

"""
    fastica!(W, X, fun, maxiter, tol, verbose)

Invoke the Fast ICA algorithm[^1].

**Parameters:**
- `W`: The initial un-mixing matrix, of size ``(m, k)``. The function updates this matrix inplace.
- `X`: The data matrix, of size ``(m, n)``. This matrix is input only, and won't be modified.
- `fun`: The approximate neg-entropy functor of type [`ICAGDeriv`](@ref).
- `maxiter`: Maximum number of iterations.
- `tol`: Tolerable change of `W` at convergence.

Returns the updated `W`.

**Note:** The number of components is inferred from `W` as `size(W, 2)`.
"""
function fastica!(W::DenseMatrix{T},         # initialized component matrix, size (m, k)
                  X::DenseMatrix{T},         # (whitened) observation sample matrix, size(m, n)
                  fun::ICAGDeriv,            # approximate neg-entropy functor
                  maxiter::Int,              # maximum number of iterations
                  tol::Real) where {T<:Real} # convergence tolerance

    # argument checking
    m = size(W, 1)
    k = size(W, 2)
    size(X, 1) == m || throw(DimensionMismatch("Sizes of W and X mismatch."))
    n = size(X, 2)
    k <= min(m, n) || throw(DimensionMismatch("k must not exceed min(m, n)."))

    @debug "FastICA Algorithm" m=m n=n k=k

    # pre-allocated storage
    Wp = similar(W)                # to store previous version of W
    U  = Matrix{T}(undef, n, k)    # to store w'x & g(w'x)
    Y  = Matrix{T}(undef, m, k)    # to store E{x g(w'x)} for components
    E1 = Vector{T}(undef, k)       # store E{g'(w'x)} for components

    # normalize each column
    for j = 1:k
        w = view(W,:,j)
        rmul!(w, 1.0 / sqrt(sum(abs2, w)))
    end

    # main loop
    chg = T(NaN)
    t = 0
    converged = false
    while !converged && t < maxiter
        t += 1
        copyto!(Wp, W)

        # apply W of previous step
        mul!(U, transpose(X), W) # u <- w'x

        # compute g(w'x) --> U and E{g'(w'x)} --> E1
        update!(fun, U, E1)

        # compute E{x g(w'x)} --> Y
        rmul!(mul!(Y, X, U), 1 / n)

        # update w: E{x g(w'x)} - E{g'(w'x)}w, i.e. w := y - e1 * w
        for j = 1:k
            w = view(W,:,j)
            y = view(Y,:,j)
            e1 = E1[j]
            @. w = y - e1 * w
        end

        # symmetric decorrelation: W <- W * (W'W)^{-1/2}
        copyto!(W, W * _invsqrtm!(W'W))

        # compare with Wp to evaluate a conversion change
        chg = maximum(abs.(abs.(diag(W*Wp')) .- 1))
        converged = (chg < tol)

        @debug "Iteration $t" change=chg tolerance=tol
    end
    converged || throw(ConvergenceException(maxiter, chg, oftype(chg, tol)))
    return W
end

#### interface function
"""
    fit(ICA, X, k; ...)

Perform ICA over the data set given in `X`.

**Parameters:**
-`X`: The data matrix, of size ``(m, n)``. Each row corresponds to a mixed signal,
while each column corresponds to an observation (*e.g* all signal value at a particular time step).
-`k`: The number of independent components to recover.

**Keyword Arguments:**
- `alg`: The choice of algorithm (*default* `:fastica`)
- `fun`: The approx neg-entropy functor (*default* [`Tanh`](@ref))
- `do_whiten`: Whether to perform pre-whitening (*default* `true`)
- `maxiter`: Maximum number of iterations (*default* `100`)
- `tol`: Tolerable change of ``W`` at convergence (*default* `1.0e-6`)
- `mean`: The mean vector, which can be either of:
    - `0`: the input data has already been centralized
    - `nothing`: this function will compute the mean (*default*)
    - a pre-computed mean vector
- `winit`: Initial guess of ``W``, which should be either of:
    - empty matrix: the function will perform random initialization (*default*)
    - a matrix of size ``(k, k)`` (when `do_whiten`)
    - a matrix of size ``(m, k)`` (when `!do_whiten`)

Returns the resultant ICA model, an instance of type [`ICA`](@ref).

**Note:** If `do_whiten` is `true`, the return `W` satisfies
``\\mathbf{W}^T \\mathbf{C} \\mathbf{W} = \\mathbf{I}``,
otherwise ``W`` is orthonormal, *i.e* ``\\mathbf{W}^T \\mathbf{W} = \\mathbf{I}``.
"""
function fit(::Type{ICA}, X::AbstractMatrix{T},# sample matrix, size (m, n)
             k::Int;                           # number of independent components
             alg::Symbol=:fastica,             # choice of algorithm
             fun::ICAGDeriv=Tanh(one(T)),      # approx neg-entropy functor
             do_whiten::Bool=true,             # whether to perform pre-whitening
             maxiter::Integer=100,             # maximum number of iterations
             tol::Real=1.0e-6,                 # convergence tolerance
             mean=nothing,                     # pre-computed mean
             winit::Matrix{T}=zeros(T,0,0)     # init guess of W, size (m, k)
            ) where {T<:Real}

    # check input arguments
    m, n = size(X)
    n > 1 || error("There must be more than one samples, i.e. n > 1.")
    k <= min(m, n) || error("k must not exceed min(m, n).")

    alg == :fastica || error("alg must be :fastica")
    maxiter > 1 || error("maxiter must be greater than 1.")
    tol > 0 || error("tol must be positive.")

    # preprocess data
    mv = preprocess_mean(X, mean)
    Z::Matrix{T} = centralize(X, mv)

    W0= zeros(T, 0, 0)  # whitening matrix
    if do_whiten
        C = rmul!(Z * transpose(Z), 1.0 / (n - 1))
        Efac = eigen(C)
        ord = sortperm(Efac.values; rev=true)
        (v, P) = extract_kv(Efac, ord, k)
        W0 = rmul!(P, Diagonal(1 ./ sqrt.(v)))
        Z = W0'Z
    end

    # initialize
    W = (isempty(winit) ? randn(T, size(Z,1), k) : copy(winit))

    # invoke core algorithm
    fastica!(W, Z, fun, maxiter, tol)

    # construct model
    if do_whiten
        W = W0 * W
    end
    return ICA(mv, W)
end


text:
using Jive

ignores = Set()

push!(ignores, joinpath("adapters", "hive"))

if Sys.iswindows()
    push!(ignores, joinpath("adapters", "mysql"))
    push!(ignores, joinpath("adapters", "odbc"))
    push!(ignores, joinpath("adapters", "jdbc"))
end

if haskey(ENV, "TRAVIS") && ENV["TRAVIS_OS_NAME"] == "osx"
    push!(ignores, joinpath("adapters", "mysql"))
    push!(ignores, joinpath("adapters", "odbc"))
end

using Jive.Distributed: nprocs
if nprocs() > 1 || !(get(ENV, "JIVE_PROCS", "") in ["", "0"]) || VERSION >= v"1.3.0-alpha"
    push!(ignores, joinpath("adapters", "jdbc"))
end

# LibPQ v0.11.1
#     ERROR: LoadError: ArgumentError: cannot convert NULL to string
push!(ignores, "adapters/postgresql/copy_test")

runtests(@__DIR__, skip=["adapters/mysql/options.jl", ignores...])


text:
using Random, Statistics, Distributions, Plots; pyplot()
Random.seed!(0)

mu, sig = 50, 5
dist = Normal(mu, sig)
alpha = 0.01
nMax = 40

observations = rand(dist,1)
piLarray, piUarray = [], []

for _ in 2:nMax
    xNew = rand(dist)
    push!(observations,xNew)

    xbar, sd = mean(observations), std(observations)
    n = length(observations)
    tVal = quantile(TDist(n-1),1-alpha/2)
    delta = tVal * sd * sqrt(1+1/n)
    piL, piU = xbar - delta, xbar + delta

    push!(piLarray,piL); push!(piUarray,piU)
end

scatter(1:nMax, observations,
	c=:blue, msw=0, label="Observations")
plot!(2:nMax, piUarray,
	c=:red, shape=:xcross, msw=0, label="Prediction Interval")
plot!(2:nMax, piLarray,
	c=:red, shape=:xcross, msw=0, label="",
	ylims=(0,100), xlabel="Number of observations", ylabel="Value")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
import DataFrames
import Dates
import ODBC

"""
```julia
read(configuration::Configuration,
        range::Union{Tuple{Union{Date, Nothing}, Union{Date, Nothing}}, Nothing},
        fields::Union{Vector{Symbol}, Nothing})::DataFrame
```

# Arguments

- `configuration::BiTWeather.Configuration`: A [`BiTWeather.Configuration`](@ref)
    providing the location and structure of the weather data.
- `range::Union{Tuple{Union{Dates.Date, Nothing}, Union{Dates.Date, Nothing}}, Nothing}`:
    The date range of the records to retrieve. A value of nothing for the
    lower, upper or both date bounds means there is no lower, upper or any
    date limit respectively. By selecting a date range here, you
    can do things such limit the records to the dates you need for your
    downstream processing.
- `fields::Union{Vector{Symbol}, Nothing}`: The fields to include from each
    record. A value of nothing means all fields included in
    `configure.fieldMappings` will be included. The `:dateTime` field will
    be included reguardless of the value of `fields`. By providing a list
    of fields here, you can limit the records to the fields you need for
    your downstream processing.

# Return

-  A `DataFrames.DataFrame` containing the requested weather data in SI units
    stored as Float64 and indexed by year/month/day/hour/minute/second.

Fetches and returns the weather data described by `configuration`. It
assumes the Julia environment is configured so `configure.dsn` is a DSN that is
able to connect to the SQL database containing the weather data and read the
table `configure.table`. It renames the `configure.table` fields using the
`configure.fieldMappings`. It parses the `:dateTime` field into separate
columns of integers containing the year, month, day, hour, minute and second.
All other fields are converted to SI units equivalents and stored as Float64
values. It drops any records with missing field values.

# Example

```julia
import DataFrames
import BiTWeather

configuration = BiTWeather.Configuration(
    dsn = "meteobridge",
    table = "backinthirty",
    fieldMappings = Dict(
        :dateTime =>    BiTWeather.FieldMapping(name = "DateTime",   units = nothing),
        :temperature => BiTWeather.FieldMapping(name = "TempOutNow", units = "F")
    )
)

weatherData::DataFrames.DataFrame = BiTWeather.read(configuration, nothing, nothing)
```
"""
function read(configuration::Configuration, range::Union{Tuple{Union{Dates.Date, Nothing}, Union{Dates.Date, Nothing}}, Nothing}, fields::Union{Vector{Symbol}, Nothing})::DataFrames.DataFrame

    fieldMappings::Dict{Symbol, FieldMapping} = configuration.fieldMappings

    #=
    Build the SQL query to retrieve the data columns.
    =#
    query::String = "SELECT"
    # Add the desired columns to the query.
    columns::Vector{Symbol} = []
    if isnothing(fields)
        columns = keys(fieldMappings)
    else
        for field in fields
            if field in keys(fieldMappings)
                push!(columns, field)
            end
        end
        if !(:dateTime in columns)
            push!(columns, :dateTime)
        end
    end
    filter!(x -> x != :dateTime, columns)
    query = "$(query) $(fieldMappings[:dateTime].name)"
    for column in columns
        query = "$(query), $(fieldMappings[column].name)"
    end
    query = "$(query) FROM"
    query = "$(query) $(configuration.table)"
    if !isnothing(range)
        if     !isnothing(range[1]) && !isnothing(range[2])
            query = "$(query) WHERE $(fieldMappings[:dateTime].name) BETWEEN '$(Dates.convert(Dates.DateTime, range[1]))' AND '$(Dates.convert(Dates.DateTime, range[2]))'"
        elseif !isnothing(range[1]) &&  isnothing(range[2])
            query = "$(query) WHERE $(fieldMappings[:dateTime].name) >= '$(Dates.convert(Dates.DateTime, range[1]))'"
        elseif  isnothing(range[1]) && !isnothing(range[2])
            query = "$(query) WHERE $(fieldMappings[:dateTime].name) < '$(Dates.convert(Dates.DateTime, range[2] + Dates.Day(1)))'"
        end
    end

    # Connect to the weather database.

    connection::ODBC.Connection = ODBC.Connection(configuration.dsn)
    # Read the requested fields.
    weatherData::DataFrames.DataFrame = ODBC.DBInterface.execute(connection, query) |> DataFrames.DataFrame

    # Drop records with missing values.
    DataFrames.dropmissing!(weatherData)

    # Set the field names.
    nameMappings::Dict{String, Symbol} = Dict()
    nameMappings[fieldMappings[:dateTime].name] = :dateTime
    for column in columns
        nameMappings[fieldMappings[column].name] = column
    end
    DataFrames.rename!(weatherData, nameMappings)

    # Set the field types.
    for fieldName in DataFrames.propertynames(weatherData)
        if fieldName != :dateTime
            weatherData[!, fieldName] = convert.(Float64, weatherData[!, fieldName])
        else
            weatherData[!, fieldName] = Dates.convert.(Dates.DateTime, weatherData[!, fieldName])
        end
    end

    # Set the field units.
    for fieldName in DataFrames.propertynames(weatherData)
        if fieldName != :dateTime
            weatherData[!, fieldName] = unitMappings[fieldName][fieldMappings[fieldName].units].(weatherData[!, fieldName])
        end
    end

    weatherData[!, :year] = Dates.year.(weatherData.dateTime)
    weatherData[!, :month] = Dates.month.(weatherData.dateTime)
    weatherData[!, :day] = Dates.day.(weatherData.dateTime)
    weatherData[!, :hour] = Dates.hour.(weatherData.dateTime)
    weatherData[!, :minute] = Dates.minute.(weatherData.dateTime)
    weatherData[!, :second] = Dates.second.(weatherData.dateTime)

    DataFrames.sort!(weatherData, [:year, :month, :day, :hour, :minute, :second])

    return weatherData
end

max_stars_repo_name: pebender/BiTWeather

max_issues_repo_name: pebender/BiTWeather

max_forks_repo_name: pebender/BiTWeather


text:
### This implementation is roughly in the order defined in
### `Uniform Discrete Curvelet Transform' by Truong T. Nguyen and Hervé Chauris
### Equations will be referenced out of this paper.

### Common Lisp comment style will be used here for clarity
module Curvelet
    export curveletTransform, inverseCurveletTransform
    using FFTW
    using DSP

    Curvelet_verbose = false

    ### Eqn. 13 defines β(t) s.t. it smoothly maps the range -1...1 to 0...1
    ### This function is used to specify the transition regions of the windows
    function _beta(t::AbstractFloat)
        ## While the absolute value is not included in any of the equations, it seems
        ## to be required to properly calculate values with t~ -1.0
        if(abs(t) <= 1)
            sqrt(abs(sum([-5/32 21/32 -35/32 35/32 1/2].*(t.^[7 5 3 1 0]))))
        elseif(1 <= t)
            1.0
        else
            0.0
        end
    end

    ### ηa & ηb are used to define the respective radial and angular falloff rates
    ### Within the paper the value of 0.15 was used for both, though recovery
    ### guarentees hold for more values.
    const ETA_A = 0.15
    const ETA_B = 0.15

    #### Windowing functions
    ###
    _w1p(t) = _beta((pi-abs(t))/(pi*ETA_A))
    _w0p(t) = _w1p(2*t*(1+ETA_A))
    _w0(w::Tuple{Float64,Float64}) = _w0p(w[1]).*_w0p(w[2])

    _w1(w::Tuple{Float64, Float64}) = ((1-_w0(w)^2)^0.5)*_w1p(w[1])*_w1p(w[2])

    _v1p(t::Float64,N) = _beta((2/N-1-t)/(2*ETA_B/N))*_beta((t+1)/(2ETA_B/N))

    _vp(t,N,l) = _v1p(t-(2*(l-1))/N,N)

    ## Eqn. 29 defines T(θ) as a warping function that warps the 1D function ṽ(t,l)
    ## to the angular function v(θ,l)
    function T(theta)
        if(-pi/4 <= theta && theta <= pi/4)
            tan(theta)
        elseif(pi/4 < theta && theta < pi/2)
            2-tan(pi/2-theta)
        elseif(-pi/2 < theta && theta < -pi/4)
            -2-tan(pi/2-theta)
        elseif(theta < -pi/2)
            -2
        else#  theta > pi/2
            2
        end
    end

    ### Warped version of ṽ(t,l)
    function _vl(theta,N,l)
        if(l<N)
            _vp(float(T(theta)), N, l)
        else
            _vp(T(pi/2-theta),N,2*N+1-l)
        end
    end

    ### Intermediate in defining the windowing function for l!=0
    function _ulp(w::Tuple{AbstractFloat, AbstractFloat}, N, l)
        _vl(atan2(w[2],w[1]),N,l)*_w1(w)
    end

    function _ul(w::Tuple{AbstractFloat, AbstractFloat}, N::Int, l::Int)
        if(l==0)
            _w0(w)
        else
            result = 0
            ## These extra loops should not be needed, but something is missed
            ## if these are ommited
            for a=[-2pi 0.0 2pi], b=[-2pi 0.0 2pi]
                result += _ulp((w[1]+a,w[2]+b), N, l)
            end
            result
        end
    end

    grid(x,size) = (-pi+2.0*pi/size*(x-1))

    ### Defines the window (without scaling term) exactly as described in the paper.
    ### This approach is quite expensive though, resulting ~100x the work really
    ### needed to get the output as many points outside the ROS are dealt with
    function direct_window_definition(size::Int, N::Int, l::Int)
        fn = (x,y)->_ul((x,y),N,l)

        result = zeros(size, size)

        ## Calculate the window, allowing it to have its full image expressed by
        ## wraping around the edge
        for x=1:size, y=1:size,
            result[y,x] = fn(grid(x,size), grid(y,size))
        end
        return result
    end

    ###Generate scaled window, see fast_window_definition in the optimization
    ### section
    function window(size::Int, N::Int, l::Int)
        if(l==0)
            fast_window_definition(size, N, l)*2
        else
            fast_window_definition(size, N, l)*2*sqrt(2)
        end
    end

    ### Performs the downsample by a factor of two in both directions
    function downsample(X)
        DSP.fftshift(DSP.fft(DSP.ifft(DSP.ifftshift(X))[1:2:end,1:2:end]))
    end

    function downsample2(X)
        Xx = DSP.fftshift(X)
        y  = X[1:Int(end/2),1:Int(end/2)]
        y += X[1:Int(end/2),Int(end/2)+1:end]
        y += X[Int(end/2)+1:end,Int(end/2)+1:end]
        y += X[Int(end/2)+1:end,1:Int(end/2)]
        DSP.fftshift(y)/4
    end

    function super_downsample(X)
        Xx = DSP.fftshift(X)
        y  = X[1:end/2,1:end/2]
        y += X[1:end/2,end/2+1:end]
        y += X[end/2+1:end,end/2+1:end]
        y += X[end/2+1:end,1:end/2]
        yy = y[1:end/2,:] + y[end/2+1:end,:]
        yyy = yy[1:end/2,:] + yy[end/2+1:end,:]
        DSP.fftshift(yyy)
    end

    function super_upsample(X)
        Xx = DSP.fftshift(X)
        XX = [Xx Xx; Xx Xx]
        XXX = [XX; XX]
        XXXX = [XXX; XXX]
    end

    function upsample(X)
        N = size(X,1)
        xx = zeros(2*N,2*N)+0im
        xx[1:2:end,1:2:end] = DSP.ifft(DSP.ifftshift(X))
        DSP.fftshift(DSP.fft(xx))
    end

    function upsample2(X)
        Xx = DSP.fftshift(X)
        [Xx Xx; Xx Xx]
    end

    C1 = 7
    C2 = 2*C1+1

    function _gen_windows(Wsize,N)
        win = Array{Matrix{Float64}}(2*N+1)
        for i=0:N
            if(Curvelet_verbose)
                print('.')
            end
            win[i+1] = window(Wsize,N,i)
        end
        for i=N+1:2N
            win[i+1] = flipdim(win[i-N+1],1)'
            win[i+1][:,[2:end..., 1]] = win[i+1][:, [1:end...]]
        end
        win
    end

    struct CurveletPlan
        dim::Int
        angles::Int
        windows::Vector{Matrix{Float64}}
        subplan::Union{CurveletPlan,Void}
    end

    ## no need to print out a ton of matricies...
    import Base.show
    show(io::IO, cp::Curvelet.CurveletPlan) = print(io, "CurveletPlan")

    function planCurveletTransform(dim::Int, N::Int)
        if(Curvelet_verbose)
            print("<planning>")
        end
        @assert mod(dim,2)==0
        @assert dim > 0
        @assert N > 0

        if(dim <= 16)
            return nothing
        end

        CurveletPlan(dim,2*N+1,_gen_windows(dim,N),planCurveletTransform(Int(dim/2), N))
    end

    struct CurveletCoeff
        HP::Vector{Matrix{Complex{Float64}}}
        LP::Union{CurveletCoeff,Matrix{Complex{Float64}}}
    end

    curveletTransform(x) = curveletTransform(x, planCurveletTransform(size(x,1), C1))

    function curveletTransform(x::Matrix{Float64}, plan::CurveletPlan)
        if(Curvelet_verbose)
            println("FFT")
        end
        S = DSP.fftshift(DSP.fft(x))

        tf = Array{Matrix{Complex{Float64}}}(plan.angles)
        if(Curvelet_verbose)
            print("Transforming")
        end
        for i=1:plan.angles
            if(Curvelet_verbose)
                print(".")
            end
            tf[i] = S.*plan.windows[i]
        end
        if(Curvelet_verbose)
            println()
        end

        ds = Array{Matrix{Complex{Float64}}}(plan.angles)
        if(Curvelet_verbose)
            print("Downsampling")
        end
        for i=1:plan.angles
            ds[i] = downsample2(tf[i])
            if(Curvelet_verbose)
                print(".")
            end
        end
        if(Curvelet_verbose)
            println()
        end

        coeff = Array{Matrix{Complex{Float64}}}(plan.angles)
        if(Curvelet_verbose)
            print("creating coeff")
        end
        for i=1:plan.angles
            coeff[i] = DSP.ifft(DSP.fftshift(ds[i]))
            if(Curvelet_verbose)
                print(".")
            end
        end
        if(Curvelet_verbose)
            println()
        end

        CurveletCoeff(coeff[2:end], curveletTransform(ds[1], plan.subplan))
    end

    function curveletTransform(S::Matrix{Complex{Float64}}, plan::Void)
        DSP.ifft(DSP.fftshift(S))
    end

    function curveletTransform(S::Matrix{Complex{Float64}}, plan::CurveletPlan)
        tf = Array{Matrix{Complex{Float64}}}(plan.angles)
        if(Curvelet_verbose)
            print("Transforming")
        end
        for i=1:plan.angles
            if(Curvelet_verbose)
                print(".")
            end
            tf[i] = S.*plan.windows[i]
        end
        if(Curvelet_verbose)
            println()
        end

        ds = Array{Matrix{Complex{Float64}}}(plan.angles)
        if(Curvelet_verbose)
            print("Downsampling")
        end
        for i=1:plan.angles
            ds[i] = downsample2(tf[i])
            if(Curvelet_verbose)
                print(".")
            end
        end
        if(Curvelet_verbose)
            println()
        end

        coeff = Array{Matrix{Complex{Float64}}}(plan.angles)
        if(Curvelet_verbose)
            print("creating coeff")
        end
        for i=1:plan.angles
            coeff[i] = DSP.ifft(DSP.fftshift(ds[i]))
            if(Curvelet_verbose)
                print(".")
            end
        end
        if(Curvelet_verbose)
            println()
        end

        CurveletCoeff(coeff[2:end], curveletTransform(ds[1], plan.subplan))
    end

    inverseCurveletTransform(x) = inverseCurveletTransform(x, planCurveletTransform(size(x.HP[1],1)*2, C1))

    function inverseCurveletTransform(cc::CurveletCoeff, plan::CurveletPlan)
        Wsize = plan.dim

        us = Array{Matrix{Complex{Float64}}}(plan.angles)
        if(Curvelet_verbose)
            print("upsampling")
        end
        for i=1:plan.angles-1
            us[i+1] = upsample2(DSP.fftshift(DSP.fft(cc.HP[i])))
            if(Curvelet_verbose)
                print(".")
            end
        end
        if(Curvelet_verbose)
            println()
        end
        us[1] = upsample2(inverseCurveletTransform(cc.LP, plan.subplan, true))
        if(Curvelet_verbose)
            println()
        end

        prev = zeros(Wsize,Wsize)
        if(Curvelet_verbose)
            print("Inverting")
        end
        for i=1:plan.angles
            if(Curvelet_verbose)
                print(".")
            end
            prev += us[i].*plan.windows[i]
        end
        if(Curvelet_verbose)
            println()
        end

        if(Curvelet_verbose)
            println("IFFT")
            println("imaginary : ",norm(imag(DSP.ifft(DSP.ifftshift(prev)))))
        end
        real(DSP.ifft(DSP.ifftshift(prev)))
    end

    function inverseCurveletTransform(cc::CurveletCoeff, plan::CurveletPlan, asdf::Bool)
        Wsize = plan.dim

        us = Array{Matrix{Complex{Float64}}}(plan.angles)
        if(Curvelet_verbose)
            print("upsampling")
        end
        for i=1:plan.angles-1
            us[i+1] = upsample2(DSP.fftshift(DSP.fft(cc.HP[i])))
            if(Curvelet_verbose)
                print(".")
            end
        end
        if(Curvelet_verbose)
            println()
        end
        #println("types are = ", typeof(cc.LP), " ", typeof(plan.subplan))
        us[1] = upsample2(inverseCurveletTransform(cc.LP, plan.subplan, true))
        if(Curvelet_verbose)
            println()
        end

        prev = zeros(Wsize,Wsize)
        if(Curvelet_verbose)
            print("Inverting")
        end
        for i=1:plan.angles
            if(Curvelet_verbose)
                print(".")
            end
            prev += us[i].*plan.windows[i]
        end
        if(Curvelet_verbose)
            println()
        end

        prev
    end

    function inverseCurveletTransform(cc::Matrix{Complex{Float64}}, plan::Void, asdf::Bool)
        DSP.fftshift(DSP.fft(cc))
    end

    #### Optimization methods
    ### While the overall transform as laid out above is quick, the generation of
    ### necessary windowing functions is not very fast if it is done with only first
    ### principle definitions.
    ### The overall windowing functions have a Region Of Support (ROS) which is
    ### mostly convex and slowly changing.
    ### This feature can be exploited to only calculate the points that are non-zero
    ### within the [-π,π]^2 square
    ###
    ### Through symmetry we can limit this optimization to the windows on the
    ### eastern side.
    ### From N&l we can estimate the approximate angle and use that to find a
    ### point within the ROS.
    ### From there the upper and lower bounds can be adjusted as the window is
    ### calculated column-wise.

    ### Fit the window bounds starting from the previous bounds
    function refit(f, low::Int, high::Int)
        while(f(low-1) != 0)
            low -= 1
        end
        while(f(low) == 0 && low != high)
            low += 1
        end
        while(f(high+1) != 0)
            high += 1
        end
        while(f(high) == 0 && high != low)
            high -= 1
        end

        return (low, high)
    end

    ## The oddness of the below code is either the result of some strange bug
    ## that I introduced or some WEIRD compiler bug with reguards to
    ## lambda precision...
    function run_cols(f, x, b, A, dir, N, l)
        #V = x::Float64->-pi+2.0*pi/size(A,1)*(x-1)
        M = size(A,1)
        (low, high) = refit(y->f(grid(x,M),grid(y,M)), b[1], b[2])

        for y=low:high
            x_::Int = mod(x-1, size(A,1))+1
            y_::Int = mod(y-1, size(A,2))+1
            A[y_, x_] = _ul((grid(x_,M),grid(y_,M)), N, l)
        end

        if(low != high)
            run_cols(f, x+dir, (low, high), A, dir, N, l)
        end
    end

    function fast_window_definition(size::Int, N::Int, l::Int)
        if(l>N)
            result = flipdim(fast_window_definition(size,N,N-l),1)'
            result[:,[2:end, 1]] = result[:, [1:end]]
            return result
        end
        fn = (x,y)->_ul((x,y),N,l)

        result = zeros(size, size)

        V = x->-pi+2*pi/size*(x-1)
        #f = (x,y)->fn(idxToVal(x), idxToVal(y))

        #-pi/2..3*pi/2 FIXME see eqn 28. T(theta)
        interp = x->(-1/2+1*x/N)
        slope  = (atan(pi*interp(l-1))+atan(pi*interp(l)))/2

        x::Int = l != 0 ? floor(size*4/5) : size/2 #TODO enhance this
        y::Int = round((size/(2*pi)*(V(x)*slope+pi))) + 1

        if(fn(V(x),V(y)) == 0)
            println("Warning: bad y value(",y,"), correcting...")
            ## Bad guess, which should not happen if the above math is
            ## right
            ## oh well, do a quick search to find a valid point
            for ny = 1:size
                if(fn(V(x),V(ny)) != 0)
                    y = ny
                    break
                end
            end
        end

        ## Real failure happens here
        @assert fn(V(x),V(y)) != 0

        bounds = refit(y->fn(V(x),V(y)), y, y)
        run_cols(fn, x, (y,y), result, -1, N, l)
        run_cols(fn, x+1, (y,y), result, +1, N, l)

        return result
    end
end