text:
# Test infinite parameter function methods
@testset "Infinite Parameter Function Methods" begin
    # setup the needed info
    m = InfiniteModel()
    @infinite_parameter(m, t in [0, 1])
    @infinite_parameter(m, x[1:2] in [-1, 1])
    func = InfiniteParameterFunction(sin, IC.VectorTuple(t), [1], [1])
    object = ParameterFunctionData(func, "test")
    idx = ParameterFunctionIndex(1)
    fref = ParameterFunctionRef(m, idx)
    gvref = GeneralVariableRef(m, 1, ParameterFunctionIndex)
    bad_vref = ParameterFunctionRef(m, ParameterFunctionIndex(-1))
    # JuMP.owner_model
    @testset "JuMP.owner_model" begin
        @test owner_model(fref) === m
        @test owner_model(gvref) === m
    end
    # JuMP.index
    @testset "JuMP.index" begin
        @test index(fref) == idx
        @test index(gvref) == idx
    end
    # dispatch_variable_ref
    @testset "dispatch_variable_ref" begin
        @test dispatch_variable_ref(m, idx) == fref
        @test dispatch_variable_ref(gvref) == fref
    end
    # _add_data_object
    @testset "_add_data_object" begin
        @test InfiniteOpt._add_data_object(m, object) == idx
    end
    # _data_dictionary
    @testset "_data_dictionary" begin
        @test InfiniteOpt._data_dictionary(fref) === m.param_functions
        @test InfiniteOpt._data_dictionary(gvref) === m.param_functions
    end
    # JuMP.is_valid
    @testset "JuMP.is_valid" begin
        @test is_valid(m, fref)
        @test is_valid(m, gvref)
    end
    # _data_object
    @testset "_data_object" begin
        @test InfiniteOpt._data_object(fref) === object
        @test InfiniteOpt._data_object(gvref) === object
        @test_throws ErrorException InfiniteOpt._data_object(bad_vref)
    end
    # _core_variable_object
    @testset "_core_variable_object" begin
        @test InfiniteOpt._core_variable_object(fref) === func
        @test InfiniteOpt._core_variable_object(gvref) === func
    end
    # _object_numbers
    @testset "_object_numbers" begin
        @test InfiniteOpt._object_numbers(fref) == [1]
    end
    # _parameter_numbers
    @testset "_parameter_numbers" begin
        @test InfiniteOpt._parameter_numbers(fref) == [1]
    end
    # JuMP.name
    @testset "JuMP.name" begin
        @test name(fref) == "test"
        @test name(gvref) == "test"
        @test name(bad_vref) == ""
    end
    # JuMP.set_name
    @testset "JuMP.set_name" begin
        @test isa(set_name(fref, "new"), Nothing)
        @test name(fref) == "new"
        @test isa(set_name(gvref, "new2"), Nothing)
        @test name(fref) == "new2"
        @test_throws ErrorException set_name(bad_vref, "")
    end
    # _measure_dependencies
    @testset "_measure_dependencies" begin
        @test InfiniteOpt._measure_dependencies(fref) == MeasureIndex[]
        @test InfiniteOpt._measure_dependencies(gvref) == MeasureIndex[]
    end
    # _constraint_dependencies
    @testset "_constraint_dependencies" begin
        @test InfiniteOpt._constraint_dependencies(fref) == ConstraintIndex[]
        @test InfiniteOpt._constraint_dependencies(gvref) == ConstraintIndex[]
    end
    # _reduced_variable_dependencies
    @testset "_reduced_variable_dependencies" begin
        @test InfiniteOpt._reduced_variable_dependencies(fref) == ReducedVariableIndex[]
        @test InfiniteOpt._reduced_variable_dependencies(gvref) == ReducedVariableIndex[]
    end
    # _derivative_dependencies
    @testset "_derivative_dependencies" begin
        @test InfiniteOpt._derivative_dependencies(fref) == DerivativeIndex[]
        @test InfiniteOpt._derivative_dependencies(gvref) == DerivativeIndex[]
    end
    # raw_parameter_refs
    @testset "raw_parameter_refs" begin
        @test raw_parameter_refs(fref) == IC.VectorTuple(t)
        @test raw_parameter_refs(gvref) == IC.VectorTuple(t)
    end
    # parameter_refs
    @testset "parameter_refs" begin
        @test parameter_refs(fref) == (t, )
        @test parameter_refs(gvref) == (t, )
    end
    # parameter_list
    @testset "parameter_list" begin
        @test parameter_list(fref) == [t]
        @test parameter_list(gvref) == [t]
    end
    # raw_function
    @testset "raw_function" begin
        @test raw_function(fref) == sin
        @test raw_function(gvref) == sin
    end
    # test used_by_reduced_variable
    @testset "used_by_reduced_variable" begin
        @test !used_by_reduced_variable(fref)
        push!(InfiniteOpt._reduced_variable_dependencies(fref),
              ReducedVariableIndex(1))
        @test used_by_reduced_variable(fref)
        empty!(InfiniteOpt._reduced_variable_dependencies(fref))
    end
    # test used_by_derivative
    @testset "used_by_derivative" begin
        @test !used_by_derivative(fref)
        push!(InfiniteOpt._derivative_dependencies(fref), DerivativeIndex(1))
        @test used_by_derivative(fref)
        empty!(InfiniteOpt._derivative_dependencies(fref))
    end
    # test used_by_measure
    @testset "used_by_measure" begin
        @test !used_by_measure(fref)
        push!(InfiniteOpt._measure_dependencies(fref), MeasureIndex(1))
        @test used_by_measure(fref)
        empty!(InfiniteOpt._measure_dependencies(fref))
    end
    # test used_by_constraint
    @testset "used_by_constraint" begin
        @test !used_by_constraint(fref)
        push!(InfiniteOpt._constraint_dependencies(fref), ConstraintIndex(1))
        @test used_by_constraint(fref)
        empty!(InfiniteOpt._constraint_dependencies(fref))
    end
    # test is_used
    @testset "is_used" begin
        @test !is_used(fref)
        push!(InfiniteOpt._constraint_dependencies(fref), ConstraintIndex(1))
        @test is_used(fref)
        empty!(InfiniteOpt._constraint_dependencies(fref))
    end
    # test _update_param_var_mapping
    @testset "_update_param_var_mapping" begin
        @test InfiniteOpt._update_param_var_mapping(fref, IC.VectorTuple(t)) isa Nothing
        @test used_by_parameter_function(t)
        empty!(InfiniteOpt._parameter_function_dependencies(t))
    end
    # _delete_data_object
    @testset "_delete_data_object" begin
        @test InfiniteOpt._delete_data_object(fref) isa Nothing
        @test length(InfiniteOpt._data_dictionary(fref)) == 0
        @test !is_valid(m, fref)
    end
    # test build_parameter_function
    @testset "build_parameter_function" begin
        # test errors
        @test_throws ErrorException build_parameter_function(error, sin, x[1])
        @test_throws ErrorException build_parameter_function(error, sin, (t, x))
        # test normal
        @test build_parameter_function(error, (a, b) -> 2, (t, x)) isa InfiniteParameterFunction
    end
    # test add_parameter_function
    @testset "add_parameter_function" begin
        # test does not belong
        @infinite_parameter(InfiniteModel(), t2 in [0, 1])
        func = build_parameter_function(error, (a, b) -> 2, (t2, x))
        @test_throws VariableNotOwned{GeneralVariableRef} add_parameter_function(m, func)
        # test normal
        func = build_parameter_function(error, (a, b) -> 2, (t, x))
        fref = GeneralVariableRef(m, 2, ParameterFunctionIndex)
        @test add_parameter_function(m, func, "test") == fref
        @test name(fref) == "test"
        @test parameter_refs(fref) == (t, x)
        @test used_by_parameter_function(t)
        # test default name
        func = build_parameter_function(error, cos, t)
        fref = GeneralVariableRef(m, 3, ParameterFunctionIndex)
        @test add_parameter_function(m, func) == fref
        @test name(fref) == "cos"
    end
    # test parameter_function
    @testset "parameter_function" begin
        @test parameter_function(sin, t) isa GeneralVariableRef
        @test parameter_function(sin, t, name = "name") isa GeneralVariableRef
    end
    # test making other objects
    @testset "Other Objects" begin
        f = parameter_function((a,b) -> 2, (t, x))
        # test making reduced variable
        d = Dict{Int, Float64}(1 => 0)
        @test add_variable(m, build_variable(error, f, d)) isa GeneralVariableRef
        # test making derivative
        @test deriv(f, t) isa GeneralVariableRef
    end
end

# Test _all_function_variables
@testset "_all_function_variables" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    data = TestData(par, 0, 5)
    meas = Measure(hold, data, Int[], Int[], true)
    object = MeasureData(meas, "test")
    mindex = MeasureIndex(1)
    @test InfiniteOpt._add_data_object(m, object) == mindex
    meas = InfiniteOpt._make_variable_ref(m, mindex)
    dinf = @deriv(inf, par)
    # test for variable reference
    @testset "Variable" begin
        @test InfiniteOpt._all_function_variables(par) == [par]
        @test InfiniteOpt._all_function_variables(inf) == [inf]
        @test InfiniteOpt._all_function_variables(meas) == [meas]
        @test InfiniteOpt._all_function_variables(dinf) == [dinf]
    end
    # test for GenericAffExpr
    @testset "AffExpr" begin
        # make expressions
        aff1 = meas + 2par + hold - dinf
        aff2 = zero(GenericAffExpr{Float64, GeneralVariableRef})
        # test expressions
        @test isempty(setdiff(InfiniteOpt._all_function_variables(aff1),
                              [meas, par, hold, dinf]))
        @test InfiniteOpt._all_function_variables(aff2) == GeneralVariableRef[]
    end
    # test for GenericQuadExpr
    @testset "QuadExpr" begin
        # make expressions
        quad1 = pt^2 + inf * pt - meas + 2par + hold - dinf
        quad2 = pt^2 + inf * pt
        quad3 = zero(GenericQuadExpr{Float64, GeneralVariableRef})
        # test expressions
        @test isempty(setdiff(InfiniteOpt._all_function_variables(quad1),
                      [meas, par, hold, dinf, pt, inf]))
        @test isempty(setdiff(InfiniteOpt._all_function_variables(quad2),
                              [pt, inf]))
        @test InfiniteOpt._all_function_variables(quad3) == GeneralVariableRef[]
    end
    # test backup
    @testset "Fallback" begin
        @variable(Model(), x)
        @test_throws ErrorException InfiniteOpt._all_function_variables(x)
    end
end

# Test comparisons
@testset "Comparisons" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_parameter(m, 0 <= par2 <= 1)
    @infinite_variable(m, inf(par))
    @infinite_variable(m, inf2(par, par2))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    # test AffExpr comparison
    @testset "Base.:(==) AffExpr" begin
        @test par + par2 + inf - 2 == par + (par2 + inf) - 2
        @test 0.25par- inf == 0.25par - inf
        @test inf + 3 - inf != par + inf
        @test inf + 3 - inf != inf + 3 - hold
    end
    # test QuadExpr comparison
    @testset "Base.:(==) QuadExpr" begin
        @test par * inf + par2 + inf - 2 == par * inf + (par2 + inf) - 2
        @test inf * inf - inf == inf * inf - inf
        @test inf * inf + 3 - inf != par * inf2 + inf
        @test par * par2 + inf * inf2 == par * par2 + inf * inf2
        @test par * par2 + inf * inf2 != par * par2
        @test par * par2 + 2 * inf * inf2 !=  par * par2 + inf * inf2
        @test par * inf + par2 + inf - 2 != par * inf + (par2 + inf) - 3
    end
end

# Test _object_numbers
@testset "_object_numbers" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_parameter(m, 0 <= pars[1:2] <= 1)
    @infinite_variable(m, inf(par))
    @infinite_variable(m, inf2(par, pars))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    var = build_variable(error, inf2, Dict{Int, Float64}(1 => 0.5), check = false)
    red = add_variable(m, var)
    # test for finite variable reference
    @testset "FiniteVariable" begin
        @test InfiniteOpt._object_numbers(pt) == []
        @test InfiniteOpt._object_numbers(hold) == []
    end
    # test for infinite variable reference
    @testset "InfiniteVariable" begin
        @test InfiniteOpt._object_numbers(inf) == [1]
    end
    # test for parameter reference
    @testset "Parameter" begin
        @test InfiniteOpt._object_numbers(par) == [1]
        @test InfiniteOpt._object_numbers(pars[1]) == [2]
    end
    # test for reduced infinite variable reference
    @testset "ReducedInfinite" begin
        @test InfiniteOpt._object_numbers(red) == [2]
    end
    # test for GenericAffExpr
    @testset "AffExpr" begin
        # make expressions
        aff1 = inf + inf2 + pt - 3
        aff2 = pt + hold - 2
        # test expressions
        @test sort!(InfiniteOpt._object_numbers(aff1)) == [1, 2]
        @test InfiniteOpt._object_numbers(aff2) == []
    end
    # test for GenericQuadExpr
    @testset "QuadExpr" begin
        # make expressions
        quad1 = inf * inf2 + inf + inf2 + pt - 3 - par
        quad2 = pt * pt + pt + hold - 2
        # test expressions
        @test sort!(InfiniteOpt._object_numbers(quad1)) == [1, 2]
        @test InfiniteOpt._object_numbers(quad2) == []
    end
end

# Test _parameter_numbers
@testset "_parameter_numbers" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_parameter(m, 0 <= pars[1:2] <= 1)
    @infinite_variable(m, inf(par))
    @infinite_variable(m, inf2(par, pars))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    var = build_variable(error, inf2, Dict{Int, Float64}(1 => 0.5), check = false)
    red = add_variable(m, var)
    # test for finite variable reference
    @testset "FiniteVariable" begin
        @test InfiniteOpt._parameter_numbers(pt) == []
        @test InfiniteOpt._parameter_numbers(hold) == []
    end
    # test for infinite variable reference
    @testset "InfiniteVariable" begin
        @test InfiniteOpt._parameter_numbers(inf) == [1]
        @test InfiniteOpt._parameter_numbers(inf2) == [1, 2, 3]
    end
    # test for parameter reference
    @testset "Parameter" begin
        @test InfiniteOpt._parameter_numbers(par) == [1]
        @test InfiniteOpt._parameter_numbers(pars[2]) == [3]
    end
    # test for reduced infinite variable reference
    @testset "ReducedInfinite" begin
        @test InfiniteOpt._parameter_numbers(red) == [2, 3]
    end
    # test for GenericAffExpr
    @testset "AffExpr" begin
        # make expressions
        aff1 = inf + inf2 + pt - 3
        aff2 = pt + hold - 2
        # test expressions
        @test sort!(InfiniteOpt._parameter_numbers(aff1)) == [1, 2, 3]
        @test InfiniteOpt._parameter_numbers(aff2) == []
    end
    # test for GenericQuadExpr
    @testset "QuadExpr" begin
        # make expressions
        quad1 = inf * inf2 + inf + inf2 + pt - 3 - par
        quad2 = pt * pt + pt + hold - 2
        # test expressions
        @test sort!(InfiniteOpt._parameter_numbers(quad1)) == [1, 2, 3]
        @test InfiniteOpt._parameter_numbers(quad2) == []
    end
end

# Test _model_from_expr
@testset "_model_from_expr" begin
    # initialize model and references
    m = InfiniteModel()
    @hold_variable(m, hd)
    # test for variable reference
    @testset "Variable" begin
        @test InfiniteOpt._model_from_expr(hd) === m
    end
    # test for GenericAffExpr
    @testset "AffExpr" begin
        # make expressions
        aff1 = hd + 2
        aff2 = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
        # test expressions
        @test InfiniteOpt._model_from_expr(aff1) === m
        @test InfiniteOpt._model_from_expr(aff2) isa Nothing
    end
    # test for GenericQuadExpr
    @testset "QuadExpr" begin
        # make expressions
        quad1 = hd * hd + hd + 1
        quad2 = zero(JuMP.GenericQuadExpr{Float64, GeneralVariableRef})
        quad3 = hd * hd + 1
        # test expressions
        @test InfiniteOpt._model_from_expr(quad1) === m
        @test InfiniteOpt._model_from_expr(quad2) isa Nothing
        @test InfiniteOpt._model_from_expr(quad3) === m
    end
    # test for Vector{GeneralVariableRef}
    @testset "Vector{GeneralVariableRef}" begin
        vrefs1 = GeneralVariableRef[]
        vrefs2 = [hd]
        # test expressions
        @test InfiniteOpt._model_from_expr(vrefs1) isa Nothing
        @test InfiniteOpt._model_from_expr(vrefs2) === m
    end
    # test Fallback
    @testset "Fallback" begin
        @test_throws ErrorException InfiniteOpt._model_from_expr(42)
    end
end

# Test _remove_variable
@testset "_remove_variable" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    # test for GenericAffExpr
    @testset "AffExpr" begin
        # make expressions
        aff1 = pt + 2par + hold
        aff2 = zero(GenericAffExpr{Float64, GeneralVariableRef})
        # test expressions
        @test isa(InfiniteOpt._remove_variable(aff1, hold), Nothing)
        @test !haskey(aff1.terms, hold)
        @test isa(InfiniteOpt._remove_variable(aff2, inf), Nothing)
    end
    # test for GenericQuadExpr
    @testset "QuadExpr" begin
        # make expressions
        quad1 = pt^2 + inf * pt + 2par + hold
        quad2 = zero(GenericQuadExpr{Float64, GeneralVariableRef})
        # test expressions
        quad = copy(quad1)
        @test isa(InfiniteOpt._remove_variable(quad, inf), Nothing)
        @test !haskey(quad.aff.terms, inf)
        @test !haskey(quad.terms, UnorderedPair{GeneralVariableRef}(inf, pt))
        quad = copy(quad1)
        @test isa(InfiniteOpt._remove_variable(quad, pt), Nothing)
        @test !haskey(quad.terms, UnorderedPair{GeneralVariableRef}(inf, pt))
        @test !haskey(quad.terms, UnorderedPair{GeneralVariableRef}(pt, pt))
        @test isa(InfiniteOpt._remove_variable(quad2, inf), Nothing)
    end
end

# Test _set_variable_coefficient!
@testset "_set_variable_coefficient!" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_parameter(m, 0 <= par2 <= 1)
    @infinite_variable(m, x(par))
    @infinite_variable(m, y(par, par2))
    @hold_variable(m, z)
    # test with GeneralVariableRef
    @testset "GeneralVariableRef" begin
        @test InfiniteOpt._set_variable_coefficient!(x, x, 2) == 2 * x
        @test InfiniteOpt._set_variable_coefficient!(z, x, 2) == z + 2x
    end
    # test with GenericAffExpr
    @testset "AffExpr" begin
        @test InfiniteOpt._set_variable_coefficient!(x + z, x, 2) == 2x + z
        @test InfiniteOpt._set_variable_coefficient!(x + z, y, 2) == x + z + 2y
    end
    # test with GenericQuadExpr
    @testset "QuadExpr" begin
        @test InfiniteOpt._set_variable_coefficient!(y ^2 + x, x, 2) == y^2 + 2x
        @test InfiniteOpt._set_variable_coefficient!(y^2 + x, y, 2) == y^2 + x + 2y
    end
    # test fallabck
    @testset "Fallback" begin
        @test_throws ErrorException InfiniteOpt._set_variable_coefficient!(2, x, 2)
    end
end

# Test parameter reference methods
@testset "Parameter References" begin
    m = InfiniteModel()
    @independent_parameter(m, t in [0, 1])
    @independent_parameter(m, y in [0, 1])
    @dependent_parameters(m, x[1:3] in [0, 1])
    @hold_variable(m, z)
    @expression(m, c1, 2z)
    @expression(m, c2, z + t + x[1])
    # test _make_param_tuple_element (IndependentParameterIndex)
    @testset "_make_param_tuple_element (Independent)" begin
        @test InfiniteOpt._make_param_tuple_element(m, index(t)) == t
        @test InfiniteOpt._make_param_tuple_element(m, index(y)) == y
    end
    # test _make_param_tuple_element (DependentParametersIndex)
    @testset "_make_param_tuple_element (Dependent)" begin
        obj_idx = index(first(x)).object_index
        @test InfiniteOpt._make_param_tuple_element(m, obj_idx) == x
    end
    # test parameter_refs
    @testset "parameter_refs" begin
        @test parameter_refs(c1) == ()
        @test parameter_refs(c2) == (t, x)
        @test parameter_refs(zero(AffExpr)) == ()
    end
end


text:
include("common.jl")

function getThreads(dbg)
    vo = run_expr(dbg,"kproc->p_threads.arr.num")
    num = ValueObjectToJulia(vo)
    threads = Array(Uint32,num)
    for i = 0:num-1
        threads[i+1] = ValueObjectToJulia(run_expr(dbg,"(uint32_t)kproc->p_threads.arr.v[$i]"))
    end
    threads
end

immutable Thread
    name
end

all_threads = Thread[]
active_threads = Dict{Uint32,Any}()
blocked_threads = Dict{Uint32,Int}()
block_events = Array((Int,Int,Int),0)

function thread_get_name(thread)
    return "no name"
end

# Record process creations
function f1(env, ctx, id, loc)
    ctx = pcpp"lldb_private::StoppointCallbackContext"(ctx)
    exe_ctx = @cxx &ctx->exe_ctx_ref
    try
        proc = run_expr(dbg,"(uint32_t)proc",exe_ctx)
        thread = run_expr(dbg,"(uint32_t)t",exe_ctx)
        if haskey(active_threads, thread)
            # Got back to the same place
            # (can happen if an exception happens at inconvenient places)
            return false
        end
        println("Adding new thread $thread to $proc")
        push!(all_threads, Thread(thread_get_name(thread)))
        id = length(all_threads)
        active_threads[thread] = id
    catch e
        println(e)
        return true
    end
    return false
end
SetBreakpoint(f1,target,"proc_addthread")

# Record process destruction
function f2(env, ctx, id, loc)
    ctx = pcpp"lldb_private::StoppointCallbackContext"(ctx)
    try
        thread = run_expr(dbg,"(uint32_t)t",@cxx &ctx->exe_ctx_ref)
        println("Removing thread $thread")
        @assert haskey(active_threads, thread)
        delete!(active_threads, thread)
    catch e
        println(e)
        return true
    end
    return false
end
SetBreakpoint(f2,target,"proc_remthread")

events = Array((Int,Int),0)

const S_RUN     = 0
const S_READY   = 1
const S_SLEEP   = 2
const S_ZOMBIE  = 3

function f3(env, ctx, id, loc)
    ctx = pcpp"lldb_private::StoppointCallbackContext"(ctx)
    exe_ctx = @cxx &ctx->exe_ctx_ref
    try
        curp = bswap(run_expr(dbg,"\$s7",exe_ctx))
        next = run_expr(dbg,"(uint32_t)next",exe_ctx)
        state = run_expr(dbg,"(uint32_t)newstate",exe_ctx)
        if !haskey(active_threads,curp) || !haskey(active_threads, next)
            return false
        end
        cur = active_threads[curp]
        next = active_threads[next]
        push!(events,(cur,next))
        if state == S_SLEEP
            id = length(events)
            blocked_threads[curp] = id
        end
    catch e
        println(e)
        return true
    end
    if (length(events) % 100) == 0
        process_and_plot() |> display
        println()
    end
    return false
end
SetBreakpointAtLoc(f3,target,"thread.c",657)

function f4(env, ctx, id, loc)
    ctx = pcpp"lldb_private::StoppointCallbackContext"(ctx)
    exe_ctx = @cxx &ctx->exe_ctx_ref
    try
        thread = run_expr(dbg,"(uint32_t)target",exe_ctx)
        if !haskey(blocked_threads,thread)
            return false
        end
        start = blocked_threads[thread]
        delete!(blocked_threads, thread)
        last = length(events)
        push!(block_events,(active_threads[thread],start,last))
    catch e
        println(e)
        return true
    end
    return false
end
SetBreakpoint(f4,target,"thread_make_runnable")

using Compose

function plot_data(events, block_events, all_threads; block_color="red", thread_color="white")
    nthreads = length(all_threads)
    npoints = length(events)
    points = Array((Float64,Float64),2*npoints)
    twidth = 1/(nthreads+2)
    theight = 1/(npoints+2)
    for (i,e) in enumerate(events)
        (a,b) = e
        points[2*i-1] = (twidth + a*twidth, i*theight)
        points[2*i] = (twidth + b*twidth, i*theight)
    end
    x = []
    for (t,s,e) in block_events
        push!(x,Compose.line([(twidth + t*twidth, s*theight),
                     (twidth + t*twidth, e*theight)]))
    end
    Compose.compose(Compose.context(),
        Compose.compose(Compose.context(),x...,Compose.stroke(block_color)),
        Compose.compose(Compose.context(),Compose.line(points),Compose.stroke(thread_color)),
        )
end

function find_mapping(events2,bmap,id)
    while !haskey(bmap,id)
        id = id+1
        if id > length(events2)
            return length(events2)
        end
    end
    bmap[id]
end

function process_data()
    global events, block_events
    events3 = collect(enumerate(events))
    events4 = filter(events3) do t
        (a,(b,c)) = t
        b != c
    end

    bmap = Dict{Int,Int}()
    for (i,(a,_)) in enumerate(events4)
        bmap[a] = i
    end

    block_events3 = filter(block_events) do t
        (a,b,c) = t
        b != c
    end

    events2 = [x for (a,x) in events4]
    block_events2 = [(a,find_mapping(events2,bmap,b),find_mapping(events2,bmap,c)) for (a,b,c) in block_events3]

    events2, block_events2
end

function process_and_plot(;args...)
    events2, block_events2 = process_data()
    plot_data(events2, block_events2, all_threads; args...)
end

#=
compose(context(),line(plot_data(events2[i:(i+200)],all_threads)),stroke("black"))

using Reel
frames = Frames(MIME("image/png"), fps=10)

for i=1:length(events2)-200
    push!(frames, compose(context(),line(plot_data(events2[i:(i+200)],all_threads)),stroke("black")))
end

frames

film = roll(render, fps=30, duration=2)
=#

ext: jl


text:
# Delete folder
if !isdir("./downloads")
    # Remove any file already in location
    rm("./downloads", force=true, recursive=true)

    # Create downloads folder
    mkdir("./downloads")

    # Download required data
    @debug "Download GSHHG data products..."
    download("http://www.soest.hawaii.edu/pwessel/gshhg/gshhg-shp-2.3.7.zip", "./downloads/gshhg-shp-2.3.7.zip")
end

# Ensure output folder exists
if !isdir("../data")
    mkdir("../data")

    # Unzip data
    @debug "Unzipping GSHHG data..."
    run(`unzip ./downloads/gshhg-shp-2.3.7.zip -d ../data`)
end


text:
#####
##### `AbstractDifferential`
#####

"""
The subtypes of `AbstractDifferential` define a custom \"algebra\" for chain
rule evaluation that attempts to factor various features like complex derivative
support, broadcast fusion, zero-elision, etc. into nicely separated parts.

In general a differential type is the type of a derivative of a value.
The type of the value is for contrast called the primal type.
Differential types correspond to primal types, although the relation is not one-to-one.
Subtypes of  `AbstractDifferential` are not the only differential types.
In fact for the most common primal types, such as `Real` or `AbstractArray{Real}` the
the differential type is the same as the primal type.

In a circular definition: the most important property of a differential is that it should
be able to be added (by defining `+`) to another differential of the same primal type.
That allows for gradients to be accumulated.

It generally also should be able to be added to a primal to give back another primal, as
this facilitates gradient descent.

All subtypes of `AbstractDifferential` implement the following operations:

 - `+(a, b)`: linearly combine differential `a` and differential `b`
 - `*(a, b)`: multiply the differential `b` by the scaling factor `a`
 -  `Base.zero(x) = Zero()`: a zero.

Further, they often implement other linear operators, such as `conj`, `adjoint`, `dot`.
Pullbacks/pushforwards are linear operators, and their inputs are often
`AbstractDifferential` subtypes.
Pullbacks/pushforwards in-turn call other linear operators on those inputs.
Thus it is desirable to have all common linear operators work on `AbstractDifferential`s.
"""
abstract type AbstractDifferential end

Base.:+(x::AbstractDifferential) = x

"""
    extern(x)

Makes a best effort attempt to convert a differential into a primal value.
This is not always a well-defined operation.
For two reasons:
 - It may not be possible to determine the primal type for a given differential.
 For example, `Zero` is a valid differential for any primal.
 - The primal type might not be a vector space, thus might not be a valid differential type.
 For example, if the primal type is `DateTime`, it's not a valid differential type as two
 `DateTime` can not be added (fun fact: `Milisecond` is a differential for `DateTime`).

Where it is defined the operation of `extern` for a primal type `P` should be
`extern(x) = zero(P) + x`.

!!! note
    Because of its limitations, `extern` should only really be used for testing.
    It can be useful, if you know what you are getting out, as it recursively removes
    thunks, and otherwise makes outputs more consistent with finite differencing.

    The more useful action in general is to call `+`, or in the case of a [`Thunk`](@ref)
    to call [`unthunk`](@ref).

!!! warning
    `extern` may return an alias (not necessarily a copy) to data
    wrapped by `x`, such that mutating `extern(x)` might mutate `x` itself.
"""
@inline extern(x) = x

@inline Base.conj(x::AbstractDifferential) = x


text:
#
# A toy database backend.
#

abstract AbstractDatabase

#
# An opaque unit of data such as an object, a person or a thing.
#

immutable Entity{name}
    id::Int
end

isentity{name}(::Type{Entity{name}}) = true
isentity(::Type) = false
classname{name}(::Type{Entity{name}}) = name
convert{name}(::Type{Entity{name}}, id::Int) = Entity{name}(id)
show{name}(io::IO, e::Entity{name}) = print(io, typeof(e), "(", e.id, ")")
show{name}(io::IO, ::Type{Entity{name}}) = print(io, ucfirst(string(name)))

#
# Describes an attribute or a relation defined on a class of entities.
#

const SelectType = Union{Void, Symbol, Tuple{Vararg{Symbol}}}
const InverseType = Union{Void, Symbol}

immutable Arrow
    name::Symbol
    output::Output
    select::SelectType
    inverse::InverseType
end

Arrow(
    name::Symbol, domain::Type;
    lunique=true, ltotal=true, runique=false, rtotal=false,
    select=nothing, inverse=nothing) =
    Arrow(name, Output(domain, OutputMode(lunique, ltotal, runique, rtotal)), select, inverse)
Arrow(name::Symbol, targetname::Symbol; props...) =
    Arrow(name, Entity{targetname}; props...)
Arrow(name::Symbol; props...) =
    Arrow(name, Entity{name}; props...)

output(a::Arrow) = a.output

show(io::IO, a::Arrow) =
    print(io,
        a.name,
        a.inverse != nothing ? " (inverse of $(odomain(a)).$(a.inverse))" : "",
        " :: ", a.output)

#
# Describes a group of homogeneous entities.
#

immutable Class
    name::Symbol
    arrows::Vector{Arrow}
    select::SelectType
    name2arrow::Dict{Symbol,Arrow}
end

Class(name::Symbol, as::Arrow...; select=nothing) =
    Class(name, collect(Arrow, as), select, Dict{Symbol,Arrow}([a.name => a for a in as]))

show(io::IO, c::Class) =
    begin
        print(io, Entity{c.name}, ":")
        for a in c.arrows
            print(io, "\n  ", a)
        end
    end

#
# All the classes in the database.
#

immutable Schema
    classes::Vector{Class}
    name2class::Dict{Symbol,Class}
end

Schema(cs::Class...) =
    Schema(collect(Class, cs), Dict{Symbol,Class}([c.name => c for c in cs]))

show(io::IO, s::Schema) = print_joined(io, s.classes, "\n")

#
# Data for a specific database instance.  In the toy backend,
# we keep all data in memory.
#

immutable Instance
    sets::Dict{Symbol, Vector}
    maps::Dict{Tuple{Symbol, Symbol}, Vector}
end

immutable ToyDatabase <: AbstractDatabase
    schema::Schema
    instance::Instance
end

show(io::IO, db::ToyDatabase) = show(io, db.schema)

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
module Normalize

using Base: String
using CSV
using DataFrames
using ExcelFiles
using OdsIO
using SciPy
using StatFiles
using XLSX

export apply, normalize, record, record_all
export are_negative_skews, are_positive_skews, is_negative_skew, is_positive_skew
export are_normal, is_normal
export get_negative_skew_transformations, get_positive_skew_transformations
export get_skew_transformations, get_stretch_skew_transformations
export kurtosis, kurtosis_error, kurtosis_stat, kurtosis_variance
export print_findings, print_skewness_kurtosis
export skewness, skewness_error, skewness_stat, skewness_variance
export replace_blank!, replace_missing!
export normal_to_csv, tabular_to_dataframe

"""
    normalize(
        gdf; normal_ratio::Real=2, dependent::Bool=false, marker::AbstractString="__"
    ) -> AbstractDict

Normalize a (grouped) data frame `gdf` for its skewness and kurtosis ratios to be within
the range of ±`normal_ratio`.

If `dependent` is `true`, then the difference between the columns of `gdf` will be
normalized.

Functions will be separately applied to `gdf`, and the column names of the transformed
(grouped) data frame will be suffixed with a string `marker` and the applied functions.

A dictionary is returned with the following key-value pairs:
- `"normalized"` => a dictionary of a collection value of named tuples containing
    skewnesses and kurtoses for each normal transformation key
- `"normal gdf"` => a list of transformed (grouped) data frames whose data are normal
- `"nonnormal gdf"` => a list of transformed (grouped) data frames whose data are nonnormal

!!! compat "Julia 1.5"
    `normalize` uses `mergewith!`, which requires at least Julia 1.5.

!!! compat "Julia 1.5"
    `normalize` uses `contains`, which requires at least Julia 1.5.

# Examples
```jldoctest
julia> df = DataFrame(a=[-0.5, 2.47, 2.54, 2.91, 3.13]);

julia> results = normalize(df);

julia> keys(results["normalized"])
KeySet for a Dict{Any, Any} with 5 entries. Keys:
  "cube"
  "square"
  "antilog"
  "reflect and invert"
  "reflect and log base 10"

julia> length(results["normal gdf"])
5

julia> results["normal gdf"][1]
5×1 DataFrame
 Row │ a__square
     │ Float64
─────┼───────────
   1 │    0.25
   2 │    6.1009
   3 │    6.4516
   4 │    8.4681
   5 │    9.7969

julia> pairs(results["normalized"]["square"][1])
pairs(::NamedTuple) with 7 entries:
  :name           => "a"
  :skewness_stat  => -1.315
  :skewness_error => 0.913
  :skewness_ratio => -1.44
  :kurtosis_stat  => 2.149
  :kurtosis_error => 2.0
  :kurtosis_ratio => 1.074

julia> typeof(results["nonnormal gdf"][1])
DataFrame

julia> df2 = DataFrame(b=[6.24, 6.34, 6.39, 2.6, 8.35], c=[5.5, 5.8, 5, 5, 7.5]);

julia> results2 = normalize(df2, dependent=true);

julia> names(results2["normal gdf"][1])
1-element Vector{String}:
 "c_minus_b__add_n_logit"

julia> df3 = DataFrame(g=[1,2,1,2,1,2,1,2], kel=[5.5,5.8,5,5,7.5,1.08,0.94,0.56]);

julia> results3 = normalize(groupby(df3, :g));

julia> results3["nonnormal gdf"]
1-element Vector{Any}:
 GroupedDataFrame with 2 groups based on key: g
First Group (4 rows): g = 1
 Row │ g      kel__add_n_logit
     │ Int64  Float64
─────┼─────────────────────────
   1 │     1         0.0829742
   2 │     1         0.0917704
   3 │     1         0.0599979
   4 │     1         0.796793
⋮
Last Group (4 rows): g = 2
 Row │ g      kel__add_n_logit
     │ Int64  Float64
─────┼─────────────────────────
   1 │     2         0.078464
   2 │     2         0.0917704
   3 │     2         0.605338
   4 │     2         0.629731

julia> names(results3["normal gdf"][1])
2-element Vector{String}:
 "g"
 "kel__square_root"

julia> for nt in results3["normalized"]["square root"]
           println(nt[:name])
       end
kel (g = 1)
kel (g = 2)
```
"""
function normalize(
    gdf; normal_ratio::Real=2, dependent::Bool=false, marker::AbstractString="__"
)
    try
        if dependent
            if _has_only_one_group_and_two_vars(gdf)
                difference = _diff_pair(gdf)
                gdf = DataFrame([difference])
            else
                throw(
                    ArgumentError(
                        "One group with only two dependent variables can be normalized."
                    )
                )
            end
        end
        transformations = get_skew_transformations(gdf; normal_ratio)
        return record_all(gdf, transformations; normal_ratio, marker)
    catch e
        println(e)
        return Dict(
            "normalized" => Dict(),
            "normal gdf" => [],
            "nonnormal gdf" => [],
        )
    end
end

function _has_only_one_group_and_two_vars(gdf)
    return gdf isa AbstractDataFrame && ncol(gdf) === 2
end

function _diff_pair(df)
    var = names(df)[1]
    var2 = names(df)[2]
    return "$(var2)_minus_$(var)" => df[!, var2] .- df[!, var]
end

"""
    get_skew_transformations(gdf; normal_ratio::Real=2) -> AbstractDict

Return functions to apply to a (grouped) data frame `gdf` if its skewness and kurtosis
ratios are not within the range of ±`normal_ratio`.

The returned dictionary has the following key-value pairs:
- `"one arg"` => a collection of functions that require one argument
- `"two args"` => a collection of tuples containing functions that require two arguments,
    and their second arguments

!!! compat "Julia 1.5"
    `get_skew_transformations` uses `mergewith!`, which requires at least Julia 1.5.

Throw `ArgumentError` if there is a column that does not have more than 3 non-`NaN` values.

See also: [`get_negative_skew_transformations`](@ref),
[`get_positive_skew_transformations`](@ref), [`get_stretch_skew_transformations`](@ref)

# Examples
```jldoctest
julia> df = DataFrame(a=[1.074, 0.058, 2.321, 0.466, 0.226, 0.171, 0.439]);

julia> transformations = get_skew_transformations(df);

julia> transformations["one arg"]
3-element Vector{Function}:
 square_root (generic function with 1 method)
 add_n_logit (generic function with 1 method)
 logit (generic function with 1 method)

julia> transformations["two args"]
5-element Vector{Any}:
 (Normalize.square_root_n_add_n_invert, 0.058)
 (Normalize.add_n_invert, 0.058)
 (Normalize.add_n_log_base_10, 0.058)
 (Normalize.add_n_natural_log, 0.058)
 (Normalize.square_n_add_n_invert, 0.058)
```
"""
function get_skew_transformations(gdf; normal_ratio::Real=2)
    transformations = Dict(
        "one arg" => Function[],
        "two args" => Any[],
    )
    skews_and_kurts = _get_skewness_kurtosis(gdf)
    df_new = DataFrame(gdf)
    if _are_nonnormal(skews_and_kurts; normal_ratio)
        if are_positive_skews(skews_and_kurts)
            positive = get_positive_skew_transformations(df_new)
            mergewith!(append!, transformations, positive)
        elseif are_negative_skews(skews_and_kurts)
            negative = get_negative_skew_transformations(df_new)
            mergewith!(append!, transformations, negative)
        end
        stretch = get_stretch_skew_transformations(df_new)
        mergewith!(append!, transformations, stretch)
    end

    return transformations
end

function _get_skewness_kurtosis(df::AbstractDataFrame)
    skews = skewness.(eachcol(df))
    kurts = kurtosis.(eachcol(df))
    return merge.(
        skews,
        kurts,
    )
end

function _get_skewness_kurtosis(gd)
    gd_new = combine(gd, valuecols(gd) .=> x -> [skewness(x), kurtosis(x),]; ungroup=false)
    return [
        merge(gd_new[group][!, var]...)
        for group in keys(gd_new)
            for var in valuecols(gd_new)
    ]
end

"""
    skewness(col, round_to::Integer=3) -> NamedTuple

Compute the skewness of a collection `col`, rounded to `round_to` digits.

The returned named tuple has the following fields:
- `skewness_stat`: skewness statistic
- `skewness_error`: skewness standard error
- `skewness_ratio`: skewness ratio

Skewness statistic is corrected for bias with the adjusted Fisher-Pearson standardized
moment coefficient.

Throw `ArgumentError` if `col` does not have more than 2 non-`NaN` values.

Source: [`scipy.stats.skew`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.skew.html)

See also: [`skewness_stat`](@ref), [`skewness_error`](@ref)

# Examples
```jldoctest
julia> skewness([1.072, 0.486, 0.133, 0.223, 1.139])
(skewness_stat = 0.306, skewness_error = 0.913, skewness_ratio = 0.335)
```
"""
function skewness(col, round_to::Integer=3)
    stat = skewness_stat(col, round_to)
    error = skewness_error(col, round_to)
    ratio = round(stat * invert(error); digits=round_to)
    return (
        skewness_stat=stat,
        skewness_error=error,
        skewness_ratio=ratio,
    )
end

"""
    skewness_stat(col, round_to::Integer=3)

Compute the skewness statistic of a collection `col`, rounded to `round_to` digits.

Skewness statistic is corrected for bias with the adjusted Fisher-Pearson standardized
moment coefficient.

Throw `ArgumentError` if `col` is empty or contains only `NaN`s.

Source: [`scipy.stats.skew`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.skew.html)

See also: [`skewness`](@ref)

# Examples
```jldoctest
julia> skewness_stat([1.072, 0.486, 0.133, 0.223, 1.139])
0.306
```
"""
function skewness_stat(col, round_to::Integer=3)
    if iszero(_length_with_nan_excluded(col))
        throw(ArgumentError("col cannot be empty nor contain only NaNs."))
    elseif any(isnan, col)
        return round(
            SciPy.stats.skew(col, bias=false, nan_policy="omit")[1]; digits=round_to
        )
    else
        return round(SciPy.stats.skew(col, bias=false, nan_policy="omit"); digits=round_to)
    end
end

function _length_with_nan_excluded(col)
    omission = count(isnan.(col))
    return length(col) - omission
end

"""
    skewness_error(col, round_to::Integer=3)

Compute the skewness error of a collection `col`, rounded to `round_to` digits.

Throw `ArgumentError` if the number of items excluding `NaN` values is 2.

See also: [`skewness`](@ref)

# Examples
```jldoctest
julia> skewness_error([1.072, 0.486, 0.133, 0.223, 1.139])
0.913
```
"""
function skewness_error(col, round_to::Integer=3)
    variance = skewness_variance(col)
    return round(√variance; digits=round_to)
end

"""
    skewness_variance(col)

Compute the skewness variance of a collection `col`.

It is calculated from the formula:

```math
\\frac{6N(N-1)}{(N-2)(N+1)(N+3)}
```

Throw `ArgumentError` if the number of items excluding `NaN` values is 2.

# Examples
```jldoctest
julia> skewness_variance([1.072, 0.486, 0.133, 0.223, 1.139])
0.8333333333333333
```
"""
function skewness_variance(col)
    N = _length_with_nan_excluded(col)
    if N === 2
        throw(ArgumentError("The number of items excluding `NaN` values cannot be 2."))
    else
        return 6 * N * (N - 1) * invert((N - 2) * (N + 1) * (N + 3))
    end
end

"""
    invert(x::Real)

Compute ``\\frac{1}{x}``.

Throw `DivideError` if `x` is 0.
"""
function invert(x::Real)
    if iszero(x)
        throw(DivideError())
    else
        return inv(x)
    end
end

"""
    kurtosis(col, round_to::Integer=3) -> NamedTuple

Compute the kurtosis of a collection `col`, rounded to `round_to` digits.

The returned named tuple has the following fields:
- `kurtosis_stat`: kurtosis statistic
- `kurtosis_error`: kurtosis standard error
- `kurtosis_ratio`: kurtosis ratio

Kurtosis statistic is corrected for bias with k statistics, and Fisher's definition is
used.

Throw `ArgumentError` if `col` does not have more than 3 non-`NaN` values.

Source: [`scipy.stats.kurtosis`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kurtosis.html)

See also: [`kurtosis_stat`](@ref), [`kurtosis_error`](@ref)

# Examples
```jldoctest
julia> kurtosis([1.072, 0.486, 0.133, 0.223, 1.139])
(kurtosis_stat = -2.952, kurtosis_error = 2.0, kurtosis_ratio = -1.476)
```
"""
function kurtosis(col, round_to::Integer=3)
    stat = kurtosis_stat(col, round_to)
    error = kurtosis_error(col, round_to)
    ratio = round(stat * invert(error); digits=round_to)
    return (
        kurtosis_stat=stat,
        kurtosis_error=error,
        kurtosis_ratio=ratio,
    )
end

"""
    kurtosis_stat(col, round_to::Integer=3)

Compute the kurtosis statistic of a collection `col`, rounded to `round_to` digits.

Kurtosis statistic is corrected for bias with k statistics, and Fisher's definition is
used.

Throw `ArgumentError` if `col` is empty or contains only `NaN`s.

Source: [`scipy.stats.kurtosis`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.kurtosis.html)

See also: [`kurtosis`](@ref)

# Examples
```jldoctest
julia> kurtosis_stat([1.072, 0.486, 0.133, 0.223, 1.139])
-2.952
```
"""
function kurtosis_stat(col, round_to::Integer=3)
    if iszero(_length_with_nan_excluded(col))
        throw(ArgumentError("col cannot be empty nor contain only NaNs."))
    elseif any(isnan, col)
        return round(
            SciPy.stats.kurtosis(col, bias=false, nan_policy="omit")[1]; digits=round_to
        )
    else
        return round(
            SciPy.stats.kurtosis(col, bias=false, nan_policy="omit"); digits=round_to
        )
    end
end

"""
    kurtosis_error(col, round_to::Integer=3)

Compute the kurtosis error of a collection `col`, rounded to `round_to` digits.

Throw `ArgumentError` if the number of items excluding `NaN` values is 2 or 3.

See also: [`kurtosis`](@ref)

# Examples
```jldoctest
julia> kurtosis_error([1.072, 0.486, 0.133, 0.223, 1.139])
2.0
```
"""
function kurtosis_error(col, round_to::Integer=3)
    variance = kurtosis_variance(col)
    return round(√variance; digits=round_to)
end

"""
    kurtosis_variance(col)

Compute the kurtosis variance of a collection `col`.

It is calculated from the formula:

```math
\\frac{4(N^2-1)*skewnessVariance}{(N-3)(N+5)}
```

Throw `ArgumentError` if the number of items excluding `NaN` values is 2 or 3.

# Examples
```jldoctest
julia> kurtosis_variance([1.072, 0.486, 0.133, 0.223, 1.139])
4.0
```
"""
function kurtosis_variance(col)
    N = _length_with_nan_excluded(col)
    if N === 3
        throw(ArgumentError("The number of items excluding `NaN` values cannot be 3."))
    end
    skewness_var = skewness_variance(col)
    return 4 * (N^2 - 1) * skewness_var * invert((N - 3) * (N + 5))
end

function _are_nonnormal(skew_kurt_ratios; normal_ratio::Real=2)
    return !are_normal(skew_kurt_ratios; normal_ratio)
end

"""
    are_normal(skew_kurt_ratios; normal_ratio::Real=2) -> Bool

Return `true` if ratios are within the range of ±`normal_ratio` in collections
`skew_kurt_ratios` containing named tuples, and `false` otherwise.

The named tuples in `skew_kurt_ratios` have the following fields:
- `skewness_ratio`: skewness ratio
- `kurtosis_ratio`: kurtosis ratio

# Examples
```jldoctest
julia> sk = [
           (skewness_ratio=1.75, kurtosis_ratio=0.5),
           (skewness_ratio=0.3, kurtosis_ratio=0.6),
       ];

julia> are_normal(sk)
true

julia> sk2 = (
           (skewness_ratio=-3.14, kurtosis_ratio=1.11),
           (skewness_ratio=1.38, kurtosis_ratio=0.66),
       );

julia> are_normal(sk2)
false
```
"""
function are_normal(skew_kurt_ratios; normal_ratio::Real=2)
    return all(
        i -> is_normal(i[:skewness_ratio], i[:kurtosis_ratio]; normal_ratio),
        skew_kurt_ratios,
    )
end

"""
    is_normal(skewness_ratio::Real, kurtosis_ratio::Real; normal_ratio::Real=2) -> Bool

Return `true` if `skewness_ratio` and `kurtosis_ratio` are within the range of
±`normal_ratio`, and `false` otherwise.

# Examples
```jldoctest
julia> is_normal(1.38, 0.66)
true

julia> is_normal(-3.14, 1.11)
false
```
"""
function is_normal(skewness_ratio::Real, kurtosis_ratio::Real; normal_ratio::Real=2)
    return (-normal_ratio ≤ skewness_ratio ≤ normal_ratio) &&
        (-normal_ratio ≤ kurtosis_ratio ≤ normal_ratio)
end

"""
    are_positive_skews(skewnesses) -> Bool

Return `true` if there are any positive, non-zero ratios in a collection `skewnesses`
containing named tuples, and `false` otherwise.

`skewnesses` has named tuples with a field `skewness_ratio`.

# Examples
```jldoctest
julia> are_positive_skews([(skewness_ratio=0.513,), (skewness_ratio=-2.049,)])
true

julia> are_positive_skews([(skewness_ratio=-1.452,), (skewness_ratio=0,)])
false
```
"""
function are_positive_skews(skewnesses)
    return _contains(is_positive_skew, skewnesses)
end

function _contains(func, skewnesses)
    return any(i -> func(i[:skewness_ratio]), skewnesses)
end

function _contains(func, df::AbstractDataFrame)
    filtered_cols = _filter_cols(func, df)
    filtered_lengths = length.(filtered_cols)
    return sum(filtered_lengths) > 0
end

function _contains(func, df::AbstractDataFrame, func_other_arg)
    filtered_cols = _filter_cols(func, df, func_other_arg)
    filtered_lengths = length.(filtered_cols)
    return sum(filtered_lengths) > 0
end

function _filter_cols(func, df)
    return filter.(func, eachcol(df))
end

function _filter_cols(func, df, func_other_arg)
    return filter.(x -> func(x, func_other_arg), eachcol(df))
end

"""
    is_positive_skew(skewness_ratio::Real) -> Bool

Return `true` if `skewness_ratio` is positive and not zero, and `false` otherwise.

# Examples
```jldoctest
julia> is_positive_skew(1.25)
true

julia> is_positive_skew(-1.25)
false

julia> is_positive_skew(0)
false
```
"""
function is_positive_skew(skewness_ratio::Real)
    return skewness_ratio > 0
end

"""
    get_positive_skew_transformations(df) -> AbstractDict

Return functions to transform positive skew for a data frame `df` based on its minimum
value.

The returned dictionary has the following key-value pairs:
- `"one arg"` => a collection of functions that require one argument
- `"two args"` => a collection of tuples containing functions that require two arguments,
    and their second arguments

# Examples
```jldoctest
julia> df = DataFrame(a=[1.072, 0.486, 0.133, 0.223, 1.139]);

julia> positive = get_positive_skew_transformations(df);

julia> positive["one arg"]
1-element Vector{Function}:
 square_root (generic function with 1 method)

julia> positive["two args"]
5-element Vector{Tuple{Function, Float64}}:
 (Normalize.square_root_n_add_n_invert, 0.133)
 (Normalize.add_n_invert, 0.133)
 (Normalize.add_n_log_base_10, 0.133)
 (Normalize.add_n_natural_log, 0.133)
 (Normalize.square_n_add_n_invert, 0.133)
```
"""
function get_positive_skew_transformations(df)
    min = _extrema(df)[1]
    if min < 0
        positive = Dict(
            "one arg" => Function[],
            "two args" => Any[
                add_n_square_root,
                add_n_square_root_n_invert,
                add_n_invert,
                add_n_log_base_10,
                add_n_natural_log,
                square_n_add_n_invert,
            ],
        )
        if _contains(_cannot_square_n_add_n_invert, df, min)
            pop!(positive["two args"])
        end
    elseif 0 ≤ min < 1
        positive = Dict(
            "one arg" => Function[
                square_root,
            ],
            "two args" => Any[
                square_root_n_add_n_invert,
                add_n_invert,
                add_n_log_base_10,
                add_n_natural_log,
                square_n_add_n_invert,
            ],
        )
    else
        positive = Dict(
            "one arg" => Function[
                square_root,
                square_root_n_invert,
                invert,
                square_n_invert,
                log_base_10,
                natural_log,
            ],
            "two args" => Any[],
        )
    end
    positive["two args"] = map(func -> (func, min), positive["two args"])
    return positive
end

function _extrema(df)
    filtered_cols = _filter_cols(isfinite, df)
    mins_from_cols = minimum.(filtered_cols)
    min = minimum(mins_from_cols)
    maxes_from_cols = maximum.(filtered_cols)
    max = maximum(maxes_from_cols)
    return (min, max)
end

# min < 0
# (add_n_square_root, add_n_square_root_n_invert, add_n_invert,
# add_n_log_base_10, add_n_natural_log, square_n_add_n_invert)

"""
    add_n_square_root(x::Real, min::Real)

Compute ``\\sqrt{x + 1 - min}``.

Throw `ArgumentError` if `min > x`.
"""
function add_n_square_root(x::Real, min::Real)
    if min > x
        throw(ArgumentError("min must be smaller."))
    else
        return √(x + 1 - min)
    end
end

"""
    add_n_square_root_n_invert(x::Real, min::Real)

Compute ``\\frac{1}{\\sqrt{x + 1 - min}}``.

Throw `ArgumentError` if `min > x`.
"""
function add_n_square_root_n_invert(x::Real, min::Real)
    return invert(add_n_square_root(x, min))
end

"""
    add_n_invert(x::Real, min::Real)

Compute ``\\frac{1}{x + 1 - min}``

Throw `ArgumentError` if `min > x`.
"""
function add_n_invert(x::Real, min::Real)
    if min > x
        throw(ArgumentError("min must be smaller."))
    else
        return invert(x + 1 - min)
    end
end

"""
    add_n_log_base_10(x::Real, min::Real)

Compute ``\\log_{10}(x + 1 - min)``.

Throw `ArgumentError` if `min > x`.
"""
function add_n_log_base_10(x::Real, min::Real)
    if min > x
        throw(ArgumentError("min must be smaller."))
    else
        return log_base_10(x + 1 - min)
    end
end

"""
    log_base_10(x::Real)

Compute the logarithm of `x` to base 10, in other words ``\\log_{10}(x)``.

Throw `DomainError` if `x ≤ 0`.
"""
function log_base_10(x::Real)
    if iszero(x)
        throw(DomainError(x))
    else
        return log10(x)
    end
end

"""
    add_n_natural_log(x::Real, min::Real)

Compute ``\\ln(x + 1 - min)``.

Throw `ArgumentError` if `min > x`.
"""
function add_n_natural_log(x::Real, min::Real)
    if min > x
        throw(ArgumentError("min must be smaller."))
    else
        return natural_log(x + 1 - min)
    end
end

"""
    natural_log(x::Real)

Compute the natural logarithm of `x`, in other words ``\\ln(x)``.

Throw `DomainError` if `x ≤ 0`.
"""
function natural_log(x::Real)
    if iszero(x)
        throw(DomainError(x))
    else
        return log(x)
    end
end

"""
    square_n_add_n_invert(x::Real, min::Real)

Compute ``\\frac{1}{x^2 + 1 - min^2}``.

Throw `DivideError` if ``x^2 + 1 - min^2`` is 0, and `ArgumentError` if `min > x`.
"""
function square_n_add_n_invert(x::Real, min::Real)
    if min > x
        throw(ArgumentError("min must be smaller."))
    else
        return invert(x^2 + 1 - min^2)
    end
end

function _cannot_square_n_add_n_invert(x::Real, min::Real)
    return x^2 + 1 - min^2 === 0
end

# 0 ≤ min < 1
# (square_root, square_root_n_add_n_invert, add_n_invert, add_n_log_base_10,
# add_n_natural_log, square_n_add_n_invert)

"""
    square_root(x::Real)

Compute ``\\sqrt{x}``.

Throw `DomainError` if `x` is negative.
"""
function square_root(x::Real)
    return √x
end

"""
    square_root_n_add_n_invert(x::Real, min::Real)

Compute ``\\frac{1}{\\sqrt{x} + 1 - \\sqrt{min}}``.

Throw `DomainError` if `x` or `min` is negative, and `ArgumentError` if `min > x`.
"""
function square_root_n_add_n_invert(x::Real, min::Real)
    if min > x
        throw(ArgumentError("min must be smaller."))
    else
        return invert(√x + 1 - √min)
    end
end

# min ≥ 1
# (square_root, square_root_n_invert, invert, square_n_invert, log_base_10,
# natural_log, invert)

"""
    square_root_n_invert(x::Real)

Compute ``\\frac{1}{\\sqrt{x}}``.

Throw `DomainError` if `x` is negative, and `DivideError` if `x` is 0.
"""
function square_root_n_invert(x::Real)
    return invert(√x)
end

"""
    square_n_invert(x::Real)

Compute ``\\frac{1}{x^2}``.

Throw `DivideError` if `x` is 0.
"""
function square_n_invert(x::Real)
    return invert(x^2)
end

"""
    are_negative_skews(skewnesses) -> Bool

Return `true` if there are any negative ratios in a collection `skewnesses` containing
named tuples, and `false` otherwise.

`skewnesses` has named tuples with a field `skewness_ratio`.

# Examples
```jldoctest
julia> are_negative_skews([(skewness_ratio=0.513,), (skewness_ratio=-2.049,)])
true

julia> are_negative_skews([(skewness_ratio=0.0,), (skewness_ratio=0.978,)])
false
```
"""
function are_negative_skews(skewnesses)
    return _contains(is_negative_skew, skewnesses)
end

"""
    is_negative_skew(skewness_ratio::Real) -> Bool

Return `true` if `skewness_ratio` is negative, and `false` otherwise.

# Examples
```jldoctest
julia> is_negative_skew(1.25)
false

julia> is_negative_skew(-1.25)
true

julia> is_negative_skew(0)
false
```
"""
function is_negative_skew(skewness_ratio::Real)
    return skewness_ratio < 0
end

"""
    get_negative_skew_transformations(df) -> AbstractDict

Return functions to transform negative skew for a data frame `df` based on its maximum
value.

The returned dictionary has the following key-value pairs:
- `"one arg"` => a collection of functions that require one argument
- `"two args"` => a collection of tuples containing functions that require two arguments,
    and their second arguments

# Examples
```jldoctest
julia> df = DataFrame(a=[1.072, 0.486, 0.133, 0.223, 1.139]);

julia> negative = get_negative_skew_transformations(df);

julia> negative["one arg"]
3-element Vector{Function}:
 square (generic function with 1 method)
 cube (generic function with 1 method)
 antilog (generic function with 1 method)

julia> negative["two args"]
3-element Vector{Tuple{Function, Float64}}:
 (Normalize.reflect_n_square_root, 1.139)
 (Normalize.reflect_n_log_base_10, 1.139)
 (Normalize.reflect_n_invert, 1.139)
```
"""
function get_negative_skew_transformations(df)
    max = _extrema(df)[2]
    negative = Dict(
        "one arg" => Function[
            square,
            cube,
            antilog,
        ],
        "two args" => Any[
            reflect_n_square_root,
            reflect_n_log_base_10,
            reflect_n_invert,
        ],
    )
    negative["two args"] = map(func -> (func, max), negative["two args"])
    return negative
end

"""
    square(x::Real)

Compute ``x^2``.
"""
function square(x::Real)
    return x^2
end

"""
    cube(x::Real)

Compute ``x^3``.
"""
function cube(x::Real)
    return x^3
end

"""
    antilog(x::Real)

Compute the antilog of `x` to base 10, in other words ``10^x``.
"""
function antilog(x::Real)
    return exp10(x)
end

"""
    reflect_n_square_root(x::Real, max::Real)

Compute ``\\sqrt{max + 1 - x}``.

Throw `ArgumentError` if `max < x`.
"""
function reflect_n_square_root(x::Real, max::Real)
    if max < x
        throw(ArgumentError("max must be greater."))
    else
        return √(max + 1 - x)
    end
end

"""
    reflect_n_log_base_10(x::Real, max::Real)

Compute ``\\log_{10}(max + 1 - x)``.

Throw `ArgumentError` if `max < x`.
"""
function reflect_n_log_base_10(x::Real, max::Real)
    if max < x
        throw(ArgumentError("max must be greater."))
    else
        return log_base_10(max + 1 - x)
    end
end

"""
    reflect_n_invert(x::Real, max::Real)

Compute ``\\frac{1}{max + 1 - x}``.

Throw `ArgumentError` if `max < x`.
"""
function reflect_n_invert(x::Real, max::Real)
    if max < x
        throw(ArgumentError("max must be greater."))
    else
        return invert(max + 1 - x)
    end
end

"""
    get_stretch_skew_transformations(df) -> AbstractDict

Return functions to stretch skew for a data frame `df`.

The returned dictionary has the following key-value pairs:
- `"one arg"` => a collection of functions that require one argument
- `"two args"` => a collection of functions that require two arguments

# Examples
```jldoctest
julia> df = DataFrame(a=[1.072, 0.486, 0.133, 0.223, 1.139]);

julia> get_stretch_skew_transformations(df)
Dict{String, Vector{T} where T} with 2 entries:
  "one arg"  => Function[add_n_logit, logit]
  "two args" => Any[]
```
"""
function get_stretch_skew_transformations(df)
    stretch = Dict(
        "one arg" => Function[
            add_n_logit,
            logit,
        ],
        "two args" => Any[],
    )
    if _contains(_cannot_add_n_logit, df)
        popfirst!(stretch["one arg"])
    end

    if _contains(_cannot_logit, df)
        pop!(stretch["one arg"])
    end

    return stretch
end

"""
    add_n_logit(x::Real)

Compute the logit of `x + 0.25` in base 10, in other words
``\\log_{10}|\\frac{x + 0.25}{1 - (x + 0.25)}|``.

Throw `DomainError` if `x` is -0.25, and `DivideError` if `x` is 0.75.
"""
function add_n_logit(x::Real)
    return logit(x + 0.25)
end

"""
    logit(x::Real)

Compute the logit of `x` to base 10, in other words ``\\log_{10}|\\frac{x}{1 - x}|``.

Throw `DomainError` if `x` is 0, and `DivideError` if `x` is 1.
"""
function logit(x::Real)
    return log_base_10(abs(x * invert(1 - x)))
end

function _cannot_add_n_logit(x::Real)
    return x === -0.25 || x === 0.75
end

function _cannot_logit(x::Real)
    return iszero(x) || isone(x)
end

"""
    record_all(
        gdf, transform_series; normal_ratio::Real=2, marker::AbstractString="__"
    ) -> AbstractDict

Compute the skewness and kurtosis for each function in a dictionary `transform_series`
applied to a (grouped) data frame `gdf`, and return a dictionary of transformed (grouped)
data frames, and of any resulting skewness and kurtosis whose ratios are within the range
of ±`normal_ratio`.

`transform_series` has the following key-value pairs:
- `"one arg"` => a collection of functions that require one argument
- `"two args"` => a collection of tuples containing functions that require two arguments,
    and their second arguments

Functions will be separately applied to `gdf`, and the column names of the transformed
(grouped) data frame will be suffixed with a string `marker` and the applied functions.

A dictionary is returned with the following key-value pairs:
- `"normalized"` => a dictionary of a collection value of named tuples containing
    skewnesses and kurtoses for each normal transformation key
- `"normal gdf"` => a list of transformed (grouped) data frames whose data are normal
- `"nonnormal gdf"` => a list of transformed (grouped) data frames whose data are nonnormal

!!! compat "Julia 1.5"
    `record_all` uses `contains`, which requires at least Julia 1.5.

Throw `ArgumentError` if there is a column that does not have more than 3 non-`NaN` values,
or if `marker` is shorter than 2 characters, or does not only consist of punctuations; and
`OverflowError` if transforming `df` or `gd` results in a value too large to represent.

See also: [`record`](@ref)

# Examples
```jldoctest
julia> df = DataFrame(a=[-0.5, 2.47, 2.54, 2.91, 3.13]);

julia> series = Dict(
           "one arg" => Function[Normalize.antilog],
           "two args" => [(Normalize.reflect_n_square_root, 3.13)],
       );

julia> records = record_all(df, series);

julia> records["normal gdf"]
1-element Vector{Any}:
 5×1 DataFrame
 Row │ a__antilog
     │ Float64
─────┼─────────────
   1 │    0.316228
   2 │  295.121
   3 │  346.737
   4 │  812.831
   5 │ 1348.96

julia> pairs(records["normalized"]["antilog"][1])
pairs(::NamedTuple) with 7 entries:
  :name           => "a"
  :skewness_stat  => 0.852
  :skewness_error => 0.913
  :skewness_ratio => 0.933
  :kurtosis_stat  => -0.029
  :kurtosis_error => 2.0
  :kurtosis_ratio => -0.014

julia> records["nonnormal gdf"]
1-element Vector{Any}:
 5×1 DataFrame
 Row │ a__reflect_n_square_root
     │ Float64
─────┼──────────────────────────
   1 │                  2.15174
   2 │                  1.28841
   3 │                  1.26095
   4 │                  1.10454
   5 │                  1.0

julia> df2 = DataFrame(g=[1,2,1,2,1,2,1,2], kel=[5.5,5.8,5,5,7.5,1.08,0.94,0.56]);

julia> series2 = Dict(
           "two args" => [(Normalize.reflect_n_log_base_10, 7.5)],
           "one arg" => [],
       );

julia> records2 = record_all(groupby(df2, :g), series2);

julia> records2["normal gdf"][1]
GroupedDataFrame with 2 groups based on key: g
First Group (4 rows): g = 1
 Row │ g      kel__reflect_n_log_base_10
     │ Int64  Float64
─────┼───────────────────────────────────
   1 │     1                    0.477121
   2 │     1                    0.544068
   3 │     1                    0.0
   4 │     1                    0.878522
⋮
Last Group (4 rows): g = 2
 Row │ g      kel__reflect_n_log_base_10
     │ Int64  Float64
─────┼───────────────────────────────────
   1 │     2                    0.431364
   2 │     2                    0.544068
   3 │     2                    0.870404
   4 │     2                    0.899821

julia> pairs(records2["normalized"]["reflect and log base 10"][1])
pairs(::NamedTuple) with 7 entries:
  :name           => "kel (g = 1)"
  :skewness_stat  => -0.577
  :skewness_error => 1.014
  :skewness_ratio => -0.569
  :kurtosis_stat  => 1.523
  :kurtosis_error => 2.619
  :kurtosis_ratio => 0.582

julia> pairs(records2["normalized"]["reflect and log base 10"][2])
pairs(::NamedTuple) with 7 entries:
  :name           => "kel (g = 2)"
  :skewness_stat  => -0.183
  :skewness_error => 1.014
  :skewness_ratio => -0.18
  :kurtosis_stat  => -4.806
  :kurtosis_error => 2.619
  :kurtosis_ratio => -1.835

julia> records2["nonnormal gdf"]
Any[]
```
"""
function record_all(
    gdf, transform_series; normal_ratio::Real=2, marker::AbstractString="__"
)
    first_record = record(gdf, transform_series["one arg"]; normal_ratio, marker)
    other_record = record(gdf, transform_series["two args"]; normal_ratio, marker)
    return _merge_results!(first_record, other_record)
end

"""
    record(
        gdf, transformations::AbstractVector{Function};
        normal_ratio::Real=2, marker::AbstractString="__"
    ) -> AbstractDict
    record(
        gdf, transformations; normal_ratio::Real=2, marker::AbstractString="__"
    ) -> AbstractDict

Compute the skewness and kurtosis for each function in `transformations` applied to a
(grouped) data frame `gdf`, and return a dictionary of transformed (grouped) data frames,
and of any resulting skewness and kurtosis whose ratios are within the range of
±`normal_ratio`.

If `transformations` is not of type `AbstractVector{Function}`, then it is a collection of
function and second argument pairings.

Functions will be separately applied to `gdf`, and the column names of the transformed
(grouped) data frame will be suffixed with a string `marker` and the applied functions.

A dictionary is returned with the following key-value pairs:
- `"normalized"` => a dictionary of a collection value of named tuples containing
    skewnesses and kurtoses for each normal transformation key
- `"normal gdf"` => a list of transformed (grouped) data frames whose data are normal
- `"nonnormal gdf"` => a list of transformed (grouped) data frames whose data are nonnormal

!!! compat "Julia 1.5"
    `record` uses `contains`, which requires at least Julia 1.5.

Throw `ArgumentError` if there is a column that does not have more than 3 non-`NaN` values,
or if `marker` is shorter than 2 characters, or does not only consist of punctuations; and
`OverflowError` if transforming `df` or `gd` results in a value too large to represent.

See also: [`record_all`](@ref)

# Examples
```jldoctest
julia> df = DataFrame(g=[1,2,1,2,1,2,1,2], kel=[5.5,5.8,5,5,7.5,1.08,0.94,0.56]);

julia> catalog = record(groupby(df, :g), [(Normalize.add_n_natural_log, 0.56)]);

julia> catalog["normal gdf"]
1-element Vector{Any}:
 GroupedDataFrame with 2 groups based on key: g
First Group (4 rows): g = 1
 Row │ g      kel__add_n_natural_log
     │ Int64  Float64
─────┼───────────────────────────────
   1 │     1                1.78171
   2 │     1                1.69378
   3 │     1                2.07191
   4 │     1                0.322083
⋮
Last Group (4 rows): g = 2
 Row │ g      kel__add_n_natural_log
     │ Int64  Float64
─────┼───────────────────────────────
   1 │     2                 1.83098
   2 │     2                 1.69378
   3 │     2                 0.41871
   4 │     2                 0.0

julia> pairs(catalog["normalized"]["add and natural log"][1])
pairs(::NamedTuple) with 7 entries:
  :name           => "kel (g = 1)"
  :skewness_stat  => -1.737
  :skewness_error => 1.014
  :skewness_ratio => -1.713
  :kurtosis_stat  => 3.271
  :kurtosis_error => 2.619
  :kurtosis_ratio => 1.249

julia> pairs(catalog["normalized"]["add and natural log"][2])
pairs(::NamedTuple) with 7 entries:
  :name           => "kel (g = 2)"
  :skewness_stat  => -0.159
  :skewness_error => 1.014
  :skewness_ratio => -0.157
  :kurtosis_stat  => -4.877
  :kurtosis_error => 2.619
  :kurtosis_ratio => -1.862

julia> catalog["nonnormal gdf"]
Any[]

julia> df2 = DataFrame(a=[-0.5, 2.47, 2.54, 2.91, 3.13]);

julia> catalog2 = record(df2, Function[Normalize.logit]);

julia> catalog2["normalized"]
Dict{Any, Any}()

julia> catalog2["normal gdf"]
Any[]

julia> catalog2["nonnormal gdf"][1]
5×1 DataFrame
 Row │ a__logit
     │ Float64
─────┼───────────
   1 │ -0.477121
   2 │  0.22538
   3 │  0.217313
   4 │  0.18286
   5 │  0.167165
```
"""
function record(
    gdf, transformations::AbstractVector{Function};
    normal_ratio::Real=2, marker::AbstractString="__"
)
    main_record = Dict(
        "normalized" => Dict(),
        "normal gdf" => [],
        "nonnormal gdf" => [],
    )
    for transformation in transformations
        results = apply(transformation, gdf; marker)
        _update_normal!(main_record, results; normal_ratio, marker)
    end

    return main_record
end

function record(gdf, transformations; normal_ratio::Real=2, marker::AbstractString="__")
    main_record = Dict(
        "normalized" => Dict(),
        "normal gdf" => [],
        "nonnormal gdf" => [],
    )
    for (transformation, extremum) in transformations
        results = apply(transformation, gdf, extremum; marker)
        _update_normal!(main_record, results; normal_ratio, marker)
    end

    return main_record
end

"""
    apply(
        func, df::AbstractDataFrame[, func_other_arg]; marker::AbstractString="__"
    ) -> NamedTuple
    apply(func, gd[, func_other_arg]; marker::AbstractString="__") -> NamedTuple

Apply a function `func` to `df` or grouped data frame `gd`, and return the skewness,
kurtosis, and transformed (grouped) data frame.

`func` is passed up to two arguments. If `func` takes two arguments, then `func_other_arg`
is the second argument.

Functions will be separately applied to `df` or `gd`, and the column names of the
transformed (grouped) data frame will be suffixed with a string `marker` and the applied
functions.

A named tuple is returned with the following fields:
- `skewness_and_kurtosis`: a named tuple with the fields `skewness_stat`,
    `skewness_error`, `skewness_ratio`, `kurtosis_stat`, `kurtosis_error`, and
    `kurtosis_ratio`
- `transformed_gdf`: a transformed (grouped) data frame

!!! compat "Julia 1.5"
    `apply` uses `contains`, which requires at least Julia 1.5.

Throw `ArgumentError` if there is a column that does not have more than 3 non-`NaN` values,
or if `marker` is shorter than 2 characters, or does not only consist of punctuations; and
`OverflowError` if transforming `df` or `gd` results in a value too large to represent.

# Examples
```jldoctest
julia> df = DataFrame(a=1:4);

julia> results = apply(Normalize.square, df);

julia> pairs(results[:skewness_and_kurtosis][1])
pairs(::NamedTuple) with 6 entries:
  :skewness_stat  => 0.709
  :skewness_error => 1.014
  :skewness_ratio => 0.699
  :kurtosis_stat  => -0.592
  :kurtosis_error => 2.619
  :kurtosis_ratio => -0.226

julia> results[:transformed_gdf]
4×1 DataFrame
 Row │ a__square
     │ Int64
─────┼───────────
   1 │         1
   2 │         4
   3 │         9
   4 │        16

julia> df2 = DataFrame(group=[1, 1, 1, 1, 2, 2, 2, 2], a=1:8);

julia> results2 = apply(Normalize.add_n_invert, groupby(df2, :group), 1);

julia> pairs(results2[:skewness_and_kurtosis][1])
pairs(::NamedTuple) with 6 entries:
  :skewness_stat  => 1.469
  :skewness_error => 1.014
  :skewness_ratio => 1.449
  :kurtosis_stat  => 2.031
  :kurtosis_error => 2.619
  :kurtosis_ratio => 0.775

julia> pairs(results2[:skewness_and_kurtosis][2])
pairs(::NamedTuple) with 6 entries:
  :skewness_stat  => 0.574
  :skewness_error => 1.014
  :skewness_ratio => 0.566
  :kurtosis_stat  => -0.625
  :kurtosis_error => 2.619
  :kurtosis_ratio => -0.239

julia> results2[:transformed_gdf]
GroupedDataFrame with 2 groups based on key: group
First Group (4 rows): group = 1
 Row │ group  a__add_n_invert
     │ Int64  Float64
─────┼────────────────────────
   1 │     1         1.0
   2 │     1         0.5
   3 │     1         0.333333
   4 │     1         0.25
⋮
Last Group (4 rows): group = 2
 Row │ group  a__add_n_invert
     │ Int64  Float64
─────┼────────────────────────
   1 │     2         0.2
   2 │     2         0.166667
   3 │     2         0.142857
   4 │     2         0.125
```
"""
function apply(func, df::AbstractDataFrame; marker::AbstractString="__")
    if _is_not_long_punctuation(marker)
        throw(
            ArgumentError(
                "marker should contain at least 2 characters, and only punctuations."
            )
        )
    end
    df_altered = rename(colname -> _rename_with("-+", colname; marker=""), df)
    df_altered = _select_finite(func, df_altered)
    rename!(colname -> _rename_with(func, colname; marker), df_altered)
    return (
        skewness_and_kurtosis=_get_skewness_kurtosis(df_altered),
        transformed_gdf=df_altered,
    )
end

function apply(func, df::AbstractDataFrame, func_other_arg; marker::AbstractString="__")
    if _is_not_long_punctuation(marker)
        throw(
            ArgumentError(
                "marker should contain at least 2 characters, and only punctuations."
            )
        )
    end
    df_altered = rename(colname -> _rename_with("-+", colname; marker=""), df)
    df_altered = _select_finite(func, df_altered, func_other_arg)
    rename!(colname -> _rename_with(func, colname; marker), df_altered)
    return (
        skewness_and_kurtosis=_get_skewness_kurtosis(df_altered),
        transformed_gdf=df_altered,
    )
end

function apply(func, gd; marker::AbstractString="__")
    if _is_not_long_punctuation(marker)
        throw(
            ArgumentError(
                "marker should contain at least 2 characters, and only punctuations."
            )
        )
    end
    gd_new = _rename_valuecols(gd, "-+"; marker="")
    gd_new = _select_finite(func, gd_new)
    gd_new = _rename_valuecols(gd_new, func; marker)
    return (
        skewness_and_kurtosis=_get_skewness_kurtosis(gd_new),
        transformed_gdf=gd_new,
    )
end

function apply(func, gd, func_other_arg; marker::AbstractString="__")
    if _is_not_long_punctuation(marker)
        throw(
            ArgumentError(
                "marker should contain at least 2 characters, and only punctuations."
            )
        )
    end
    gd_new = _rename_valuecols(gd, "-+"; marker="")
    gd_new = _select_finite(func, gd_new, func_other_arg)
    gd_new = _rename_valuecols(gd_new, func; marker)
    return (
        skewness_and_kurtosis=_get_skewness_kurtosis(gd_new),
        transformed_gdf=gd_new,
    )
end

function _is_not_long_punctuation(separator)
    return length(separator) < 2 || !all(ispunct, separator)
end

function _rename_with(fragment, name; marker::AbstractString="__")
    old_name = _get_original(name; marker)
    return "$(old_name)$(marker)$(fragment)"
end

function _get_original(colname; marker::AbstractString="__")
    return SubString(colname, 1, _find_original_end(colname; marker))
end

function _find_original_end(colname; marker::AbstractString="__")
    if length(marker) < 2
        return length(colname)
    elseif Base.contains(colname, "-+_")
        starting = findlast("-+_", colname)[1]
        return starting - 1
    else
        starting = findfirst(marker, colname)[1]
        return starting - 1
    end
end

function _rename_valuecols(gd, fragment; marker::AbstractString="__")
    grouping = groupcols(gd)
    new_names = [
        var => _rename_with(fragment, string(var); marker)
        for var in valuecols(gd)
    ]
    df_altered = rename(DataFrame(gd), new_names)
    return groupby(df_altered, grouping)
end

function _select_finite(func, df::AbstractDataFrame)
    df = select(df, names(df) .=> ByRow(func))
    if _contains(isinf, df)
        throw(OverflowError("Transforming df results in a value too large to represent."))
    else
        return df
    end
end

function _select_finite(func, df::AbstractDataFrame, func_other_arg)
    df = select(df, names(df) .=> ByRow(x -> func(x, func_other_arg)))
    if _contains(isinf, df)
        throw(OverflowError("Transforming df results in a value too large to represent."))
    else
        return df
    end
end

function _select_finite(func, gd)
    grouping = groupcols(gd)
    gd = select(gd, valuecols(gd) .=> ByRow(func); ungroup=false)
    df_example = DataFrame(gd)
    if _contains(isinf, df_example[!, Not(grouping)])
        throw(OverflowError("Transforming gd results in a value too large to represent."))
    else
        return gd
    end
end

function _select_finite(func, gd, func_other_arg)
    grouping = groupcols(gd)
    gd = select(
        gd, valuecols(gd) .=> ByRow(x -> func(x, func_other_arg)); ungroup=false
    )
    df_example = DataFrame(gd)
    if _contains(isinf, df_example[!, Not(grouping)])
        throw(OverflowError("Transforming gd results in a value too large to represent."))
    else
        return gd
    end
end

function _update_normal!(
    entries, applied; normal_ratio::Real=2, marker::AbstractString="__"
)
    if are_normal(applied[:skewness_and_kurtosis]; normal_ratio)
        _label_findings!(entries, applied; marker)
        _store_transformed!(entries, "normal gdf", applied[:transformed_gdf])
    else
        _store_transformed!(entries, "nonnormal gdf", applied[:transformed_gdf])
    end

    return nothing
end

function _label_findings!(entries, applied; marker::AbstractString="__")
    gdf_altered = applied[:transformed_gdf]
    transformations = _get_applied_function_names(gdf_altered; marker)
    original_colnames = _get_original_colnames(gdf_altered; marker)
    entries["normalized"][transformations] = _label.(
        applied[:skewness_and_kurtosis], original_colnames
    )
    return entries["normalized"][transformations]
end

function _get_applied_function_names(df::AbstractDataFrame; marker::AbstractString="__")
    new_colname = names(df)[1]
    before_marker = _find_original_end(new_colname; marker)
    after_marker = before_marker + length(marker) + 1
    function_names = SubString(new_colname, after_marker)
    return _make_legible(function_names, " n "; marker)
end

function _get_applied_function_names(gd; marker::AbstractString="__")
    new_colname = string(valuecols(gd)[1])
    before_marker = _find_original_end(new_colname; marker)
    after_marker = before_marker + length(marker) + 1
    function_names = SubString(new_colname, after_marker)
    return _make_legible(function_names, " n "; marker)
end

function _make_legible(snake_case, listing_split_on=" "; marker::AbstractString="__")
    words = _make_phrase(snake_case, "_", " ")
    listable = _make_listing(words, listing_split_on)
    extra_chars = _underscores_to_spaces(marker)
    return _make_phrase(listable, extra_chars, "; ")
end

function _make_phrase(fragment, split_on, join_with)
    words = split(fragment, split_on)
    return join(words, join_with)
end

function _make_listing(phrase, split_on=" ")
    words = split(phrase, split_on)
    if length(words) ≤ 2
        return join(words, " and ")
    else
        return join(words, ", ", ", and ")
    end
end

function _underscores_to_spaces(separator)
    chars = length(separator)
    if separator === repeat("_", chars)
        return repeat(" ", chars)
    else
        return separator
    end
end

function _get_original_colnames(gdf::AbstractDataFrame; marker::AbstractString="__")
    return _get_original.(names(gdf); marker)
end

function _get_original_colnames(gdf; marker::AbstractString="__")
    originals = _get_original.(string(var) for var in valuecols(gdf); marker)
    return _concat_groupnames(gdf, originals)
end

function _concat_groupnames(gd, colnames)
    new_names = []
    for varname in colnames
        for group in keys(gd)
            grouping = _get_groupnames(group)
            push!(new_names, "$varname ($grouping)")
        end
    end

    return new_names
end

function _get_groupnames(group)
    if isone(length(group))
        punctuations = 2
    else
        punctuations = 1
    end

    return chop(string(group), head=length("GroupKey: ("), tail=punctuations)
end

function _label(prelim, tag)
    title = (
        name=tag,
    )
    return merge(title, prelim)
end

function _store_transformed!(entries, key, value::AbstractVector)
    return append!(entries[key], value)
end

function _store_transformed!(entries, key, value)
    return push!(entries[key], value)
end

function _merge_results!(entries, other)
    merge!(entries["normalized"], other["normalized"])
    _store_transformed!(entries, "normal gdf", other["normal gdf"])
    _store_transformed!(entries, "nonnormal gdf", other["nonnormal gdf"])
    return entries
end

"""
    tabular_to_dataframe(path, sheet::AbstractString="") -> AbstractDataFrame

Load tabular data from a string `path`, and converts it to a data frame.

Acceptable file extensions are .csv, .dta, .ods, .sav, .xls, and .xlsx. If `path` ends in
.ods, .xls, or .xlsx, the worksheet that is named `sheet` is converted.
"""
function tabular_to_dataframe(path, sheet::AbstractString="")
    if endswith(path, ".csv")
        return CSV.read(path, DataFrame)
    elseif endswith(path, ".xlsx")
        return DataFrame(XLSX.readtable(path, sheet)...)
    elseif endswith(path, ".sav") || endswith(path, ".dta")
        return DataFrame(load(path))
    elseif endswith(path, ".xls")
        return DataFrame(load(path, sheet))
    elseif endswith(path, ".ods")
        return ods_read(path; sheetName=sheet, retType="DataFrame")
    else
        println(
            "Files ending with .csv, .dta, .ods, .sav, .xls, or .xlsx can only be used."
        )
        return DataFrame()
    end
end

"""
    replace_blank!(df, colname; new_value::Real) -> AbstractVector

Replace all occurrences of blank values (i.e., " ") with `new_value` in a column named
`colname` in a data frame `df`. The column will be converted to element type `Float64`.

If a value cannot be converted to `Float64`, an error is raised.

# Examples
```jldoctest
julia> df = DataFrame(a=[" ", "7"], b=[6.39, " "], c=[" ", " "])
2×3 DataFrame
 Row │ a       b     c
     │ String  Any   String
─────┼──────────────────────
   1 │         6.39
   2 │ 7

julia> replace_blank!(df, :a, new_value=1.25)
2-element Vector{Float64}:
 1.25
 7.0

julia> replace_blank!(df, "b", new_value=8)
2-element Vector{Float64}:
 6.39
 8.0

julia> replace_blank!(df, 3, new_value=5)
2-element Vector{Float64}:
 5.0
 5.0

julia> df
2×3 DataFrame
 Row │ a        b        c
     │ Float64  Float64  Float64
─────┼───────────────────────────
   1 │    1.25     6.39      5.0
   2 │    7.0      8.0       5.0
```
"""
function replace_blank!(df, colname; new_value::Real)
    col = df[!, colname]
    if _is_from_csv_excel_or_opendoc(col)
        if _is_from_csv(col)
            replace!(col, " " => "$(new_value)")
            df[!, colname] = parse.(Float64, col)
        else
            replace!(col, " " => new_value)
            df[!, colname] = convert(Vector{Float64}, col)
        end
    end

    return df[!, colname]
end

function _is_from_csv_excel_or_opendoc(col)
    return eltype(col) === String || eltype(col) === Any
end

function _is_from_csv(col)
    nonblanks = col[col .!== " "]
    return eltype(nonblanks) === String
end

"""
    replace_missing!(df, colname; new_value::Real) -> AbstractVector

Replace all occurrences of `missing` with `new_value` in a column named `colname` in a data
frame `df`. The column will be converted to element type `Float64`.

If a value cannot be converted to `Float64`, an error is raised.

!!! compat "Julia 1.3"
    `replace_missing!` uses `nonmissingtype`, which is exported as of Julia 1.3.

# Examples
```jldoctest
julia> df = DataFrame(a=[missing, 10], b=[12.12, missing], c=[missing, missing])
2×3 DataFrame
 Row │ a        b           c
     │ Int64?   Float64?    Missing
─────┼──────────────────────────────
   1 │ missing       12.12  missing
   2 │      10  missing     missing

julia> replace_missing!(df, :a, new_value=2.5)
2-element Vector{Float64}:
  2.5
 10.0

julia> replace_missing!(df, "b", new_value=6)
2-element Vector{Float64}:
 12.12
  6.0

julia> replace_missing!(df, 3, new_value=9)
2-element Vector{Float64}:
 9.0
 9.0

julia> df
2×3 DataFrame
 Row │ a        b        c
     │ Float64  Float64  Float64
─────┼───────────────────────────
   1 │     2.5    12.12      9.0
   2 │    10.0     6.0       9.0
```
"""
function replace_missing!(df, colname; new_value::Real)
    col = df[!, colname]
    nonmissing = nonmissingtype(eltype(col))
    if nonmissing === String
        replace!(col, missing => "$(new_value)")
        df[!, colname] = parse.(Float64, col)
    else
        col = map(float, col)
        df[!, colname] = coalesce.(col, float(new_value))
    end

    return df[!, colname]
end

"""
    print_skewness_kurtosis(df::AbstractDataFrame; dependent::Bool=false)
    print_skewness_kurtosis(gd)

Print the skewness and kurtosis (statistic, standard error, ratio) of each column in `df`
or grouped data frame `gd`.

If `dependent` is `true`, then the skewness and kurtosis of the difference between the
columns of `df` will be printed.

# Examples
```jldoctest
julia> df = DataFrame(a=[1.072, 0.486, 0.133, 0.223, 1.139]);

julia> print_skewness_kurtosis(df)
--------
Variable: a
Skewness Statistic: 0.306
Skewness Std. Error: 0.913
Skewness Ratio: 0.335

Kurtosis Statistic: -2.952
Kurtosis Std. Error: 2.0
Kurtosis Ratio: -1.476

julia> df2 = DataFrame(b=[6.24, 6.34, 6.39, 2.6, 8.35], c=[5.5, 5.8, 5, 5, 7.5]);

julia> print_skewness_kurtosis(df2, dependent=true)
--------
Variable: c_minus_b
Skewness Statistic: 1.983
Skewness Std. Error: 0.913
Skewness Ratio: 2.172

Kurtosis Statistic: 4.212
Kurtosis Std. Error: 2.0
Kurtosis Ratio: 2.106

julia> df3 = DataFrame(
           group=[1,1,1,1,2,2,2,2],
           dd=[0.042, 0.849, 0.342, 0.771, 1.879, 1.289, 0.981, 1.073],
       );

julia> print_skewness_kurtosis(groupby(df3, :group))
--------
Variable: dd (group = 1)
Skewness Statistic: -0.478
Skewness Std. Error: 1.014
Skewness Ratio: -0.471

Kurtosis Statistic: -2.951
Kurtosis Std. Error: 2.619
Kurtosis Ratio: -1.127

--------
Variable: dd (group = 2)
Skewness Statistic: 1.439
Skewness Std. Error: 1.014
Skewness Ratio: 1.419

Kurtosis Statistic: 1.859
Kurtosis Std. Error: 2.619
Kurtosis Ratio: 0.71

```
"""
function print_skewness_kurtosis(df::AbstractDataFrame; dependent::Bool=false)
    try
        if dependent
            if ncol(df) === 2
                _print_dependent(df)
            else
                println("""
                    There can only be two dependent variables to measure
                    skewness and kurtosis.
                    """
                )
            end
        else
            _print_independent(df)
        end
    catch e
        println(e)
    end
end

function print_skewness_kurtosis(gd)
    try
        gd_new = combine(
            gd, valuecols(gd) .=> x -> [skewness(x), kurtosis(x),]; ungroup=false
        )
        for group in keys(gd_new)
            for var in valuecols(gd_new)
                tagged = _describe_var(gd_new, group, var)
                _print_summary(tagged)
            end
        end
    catch e
        println(e)
    end
end

function _print_dependent(df)
    difference = _diff_pair(df)
    tagged = _describe_var(difference.first, difference.second)
    _print_summary(tagged)
end

function _describe_var(colname, col)
    tag = (
        name=colname,
    )
    skew = skewness(col)
    kurt = kurtosis(col)
    return merge(tag, skew, kurt)
end

function _describe_var(gd, group_column, value_column)
    varname = chop(string(value_column), tail=length("_function"))
    grouping = _get_groupnames(group_column)
    tag = (
        name="$varname ($grouping)",
    )
    skews_and_kurts = gd[group_column][!, value_column]
    return merge(tag, skews_and_kurts...)
end

function _print_summary(factor::NamedTuple)
    println("""
        --------
        Variable: $(factor.name)
        Skewness Statistic: $(factor.skewness_stat)
        Skewness Std. Error: $(factor.skewness_error)
        Skewness Ratio: $(factor.skewness_ratio)

        Kurtosis Statistic: $(factor.kurtosis_stat)
        Kurtosis Std. Error: $(factor.kurtosis_error)
        Kurtosis Ratio: $(factor.kurtosis_ratio)
        """
    )
end

function _print_independent(df)
    for colname in names(df)
        tagged = _describe_var(colname, df[!, colname])
        _print_summary(tagged)
    end
end

"""
    print_findings(findings)

Print the normal skewness and kurtosis in a dictionary `findings`.

`findings` has the following key-value pairs:
- `"normalized"` => a dictionary of a collection value of named tuples containing
    skewnesses and kurtoses for each normal transformation key
- `"normal gdf"` => a list of transformed (grouped) data frames whose data are normal
- `"nonnormal gdf"` => a list of transformed (grouped) data frames whose data are nonnormal

# Examples
```jldoctest
julia> df = DataFrame(a=[-0.5, 2.47, 2.54, 2.91, 3.13]);

julia> print_findings(normalize(df))
APPLIED: cube
--------
Variable: a
Skewness Statistic: -0.678
Skewness Std. Error: 0.913
Skewness Ratio: -0.743

Kurtosis Statistic: 0.669
Kurtosis Std. Error: 2.0
Kurtosis Ratio: 0.334

APPLIED: square
--------
Variable: a
Skewness Statistic: -1.315
Skewness Std. Error: 0.913
Skewness Ratio: -1.44

Kurtosis Statistic: 2.149
Kurtosis Std. Error: 2.0
Kurtosis Ratio: 1.074

APPLIED: antilog
--------
Variable: a
Skewness Statistic: 0.852
Skewness Std. Error: 0.913
Skewness Ratio: 0.933

Kurtosis Statistic: -0.029
Kurtosis Std. Error: 2.0
Kurtosis Ratio: -0.014

APPLIED: reflect and invert
--------
Variable: a
Skewness Statistic: -0.626
Skewness Std. Error: 0.913
Skewness Ratio: -0.686

Kurtosis Statistic: 0.843
Kurtosis Std. Error: 2.0
Kurtosis Ratio: 0.422

APPLIED: reflect and log base 10
--------
Variable: a
Skewness Statistic: 1.567
Skewness Std. Error: 0.913
Skewness Ratio: 1.716

Kurtosis Statistic: 2.896
Kurtosis Std. Error: 2.0
Kurtosis Ratio: 1.448

```
"""
function print_findings(findings)
    normal = findings["normalized"]
    for (transformation, altered) in normal
        println("APPLIED: $transformation")
        _print_summary.(altered)
        println()
    end

    if _has_no_changes(findings)
        println("No transformations applied, so unable to normalize.")
    elseif isempty(normal)
        println("Transformations applied, but no normal results.")
    end
end

function _has_no_changes(findings)
    return isempty(findings["nonnormal gdf"]) && isempty(findings["normal gdf"])
end

"""
    normal_to_csv(path, findings; dependent::Bool=false)

Export normal data from a dictionary `findings` to a CSV file in a string `path`.

`findings` has a key `"normal gdf"` with a collection value of (grouped) data frames.

If `path` already exists (i.e., a file has the same name), then the file will be
overwritten. Otherwise, a new file will be created.

If `dependent` is `true`, then a column of zeros will also be created in the CSV file for
dependent testing of data with only one group and two dependent variables.

If there are any `NaN`, they will be replaced with `missing` in the CSV file.
"""
function normal_to_csv(path, findings; dependent::Bool=false)
    if !endswith(path, ".csv")
        println("Only a filename ending in .csv can be used.")
        return nothing
    end
    normal_data = findings["normal gdf"]
    if isempty(normal_data)
        println("There is no normal data to export.")
        return nothing
    end

    if normal_data[1] isa AbstractDataFrame
        df_normal = _flatten_dataframes(normal_data)
        if dependent
            df_normal = _store_zeros(df_normal)
        end
    else
        df_normal = _flatten_grouped_dataframes(normal_data)
    end
    CSV.write(path, df_normal, transform=(col, val) -> _nan_to(missing, val))
    return nothing
end

function _flatten_dataframes(gdfs)
    df_new = DataFrame()
    for df in gdfs
        df_new = hcat(df_new, df)
    end

    return df_new
end

function _store_zeros(df)
    df_zero = DataFrame(["for_dependent_test" => zeros(nrow(df))])
    return hcat(df, df_zero)
end

function _flatten_grouped_dataframes(gdfs)
    df_new = DataFrame()
    df_new = _store_groupcols(df_new, gdfs)
    return _store_valuecols(df_new, gdfs)
end

function _store_groupcols(df, gdfs)
    gd_example = gdfs[1]
    df_example = DataFrame(gd_example)
    for group in groupcols(gd_example)
        df = hcat(df, df_example[!, [group]])
    end

    return df
end

function _store_valuecols(df, gdfs)
    for grouped in gdfs
        df_convert = DataFrame(grouped)
        for var in valuecols(grouped)
            df = hcat(df, df_convert[!, [var]])
        end
    end

    return df
end

function _nan_to(new_value, original)
    if original === NaN
        return new_value
    else
        return original
    end
end

end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

polynomials over real domains
=#

import DynamicPolynomials
const DP = DynamicPolynomials
import Hypatia.PolyUtils: interpolate, BoxDomain, BallDomain, EllipsoidDomain

# get interpolation for a predefined real poly
function get_interp_data(
    ::Type{T},
    poly_name::Symbol,
    halfdeg::Int,
    ) where {T <: Real}
    (x, fn, dom, true_min) = real_poly_data(poly_name, T)
    (U, pts, Ps) = interpolate(dom, halfdeg)
    interp_vals = T[fn(pts[j, :]...) for j in 1:U]
    return (interp_vals, Ps, true_min)
end

# get interpolation for a random real poly in n variables of half degree
# halfdeg and use a box domain
function random_interp_data(
    ::Type{T},
    n::Int,
    halfdeg::Int,
    dom = BoxDomain{T}(-ones(T, n), ones(T, n)),
    ) where {T <: Real}
    (U, pts, Ps) = interpolate(dom, halfdeg)
    interp_vals = randn(T, U)
    true_min = T(NaN)
    return (interp_vals, Ps, true_min)
end

# predefined polynomials with known bounds
# see https://people.sc.fsu.edu/~jburkardt/py_src/polynomials/polynomials.html
function real_poly_data(polyname::Symbol, T::Type{<:Real} = Float64)
    if polyname == :butcher
        DP.@polyvar x[1:6]
        f = x[6]*x[2]^2+x[5]*x[3]^2-x[1]*x[4]^2+x[4]^3+x[4]^2-1/3*x[1]+4/3*x[4]
        dom = BoxDomain{T}(T[-1,-0.1,-0.1,-1,-0.1,-0.1],
            T[0,0.9,0.5,-0.1,-0.05,-0.03])
        true_obj = -1.4393333333
    elseif polyname == :caprasse
        DP.@polyvar x[1:4]
        f = -x[1]*x[3]^3+4x[2]*x[3]^2*x[4]+4x[1]*x[3]*x[4]^2+2x[2]*x[4]^3+
            4x[1]*x[3]+4x[3]^2-10x[2]*x[4]-10x[4]^2+2
        dom = BoxDomain{T}(fill(-T(0.5), 4), fill(T(0.5), 4))
        true_obj = -3.1800966258
    elseif polyname == :goldsteinprice
        DP.@polyvar x[1:2]
        f = (1+(x[1]+x[2]+1)^2*(19-14x[1]+3x[1]^2-14x[2]+6x[1]*x[2]+3x[2]^2))*
            (30+(2x[1]-3x[2])^2*(18-32x[1]+12x[1]^2+48x[2]-36x[1]*x[2]+27x[2]^2))
        dom = BoxDomain{T}(fill(-2, 2), fill(2, 2))
        true_obj = 3
    elseif polyname == :goldsteinprice_ball
        DP.@polyvar x[1:2]
        f = (1+(x[1]+x[2]+1)^2*(19-14x[1]+3x[1]^2-14x[2]+6x[1]*x[2]+3x[2]^2))*
            (30+(2x[1]-3x[2])^2*(18-32x[1]+12x[1]^2+48x[2]-36x[1]*x[2]+27x[2]^2))
        dom = BallDomain{T}(zeros(T, 2), 2*sqrt(T(2)))
        true_obj = 3
    elseif polyname == :goldsteinprice_ellipsoid
        DP.@polyvar x[1:2]
        f = (1+(x[1]+x[2]+1)^2*(19-14x[1]+3x[1]^2-14x[2]+6x[1]*x[2]+3x[2]^2))*
            (30+(2x[1]-3x[2])^2*(18-32x[1]+12x[1]^2+48x[2]-36x[1]*x[2]+27x[2]^2))
        centers = zeros(T, 2)
        Q = Diagonal(T(0.25) * ones(T, 2))
        dom = EllipsoidDomain{T}(centers, Q)
        true_obj = 3
    elseif polyname == :heart
        DP.@polyvar x[1:8]
        f = x[1]*x[6]^3-3x[1]*x[6]*x[7]^2+x[3]*x[7]^3-3x[3]*x[7]*x[6]^2+x[2]*
            x[5]^3-3*x[2]*x[5]*x[8]^2+x[4]*x[8]^3-3x[4]*x[8]*x[5]^2+0.9563453
        dom = BoxDomain{T}(T[-0.1,0.4,-0.7,-0.7,0.1,-0.1,-0.3,-1.1],
            T[0.4,1,-0.4,0.4,0.2,0.2,1.1,-0.3])
        true_obj = -1.36775
    elseif polyname == :lotkavolterra
        DP.@polyvar x[1:4]
        f = x[1]*(x[2]^2+x[3]^2+x[4]^2-1.1)+1
        dom = BoxDomain{T}(fill(-2, 4), fill(2, 4))
        true_obj = -20.8
    elseif polyname == :magnetism7
        DP.@polyvar x[1:7]
        f = x[1]^2+2x[2]^2+2x[3]^2+2x[4]^2+2x[5]^2+2x[6]^2+2x[7]^2-x[1]
        dom = BoxDomain{T}(-ones(T, 7), ones(T, 7))
        true_obj = -0.25
    elseif polyname == :magnetism7_ball
        DP.@polyvar x[1:7]
        f = x[1]^2+2x[2]^2+2x[3]^2+2x[4]^2+2x[5]^2+2x[6]^2+2x[7]^2-x[1]
        dom = BallDomain{T}(zeros(T, 7), sqrt(T(7)))
        true_obj = -0.25
    elseif polyname == :motzkin
        DP.@polyvar x[1:2]
        f = 1-48x[1]^2*x[2]^2+64x[1]^2*x[2]^4+64x[1]^4*x[2]^2
        dom = BoxDomain{T}(-ones(T, 2), ones(T, 2))
        true_obj = 0
    elseif polyname == :motzkin_ball
        DP.@polyvar x[1:2]
        f = 1-48x[1]^2*x[2]^2+64x[1]^2*x[2]^4+64x[1]^4*x[2]^2
        dom = BallDomain{T}(zeros(T, 2), sqrt(T(2)))
        true_obj = 0
    elseif polyname == :motzkin_ellipsoid
        # ellipsoid contains two local minima in opposite orthants
        DP.@polyvar x[1:2]
        f = 1-48x[1]^2*x[2]^2+64x[1]^2*x[2]^4+64x[1]^4*x[2]^2
        Q = T[1 1; 1 -1]
        D = T[1 0; 0 0.1]
        dom = EllipsoidDomain{T}(zeros(T, 2), Symmetric(Q * D * Q))
        true_obj = 0
    elseif polyname == :reactiondiffusion
        DP.@polyvar x[1:3]
        f = -x[1]+2x[2]-x[3]-0.835634534x[2]*(1+x[2])
        dom = BoxDomain{T}(fill(-5, 3), fill(5, 3))
        true_obj = -36.71269068
    elseif polyname == :robinson
        DP.@polyvar x[1:2]
        f = 1+x[1]^6+x[2]^6-x[1]^4*x[2]^2+x[1]^4-x[1]^2*x[2]^4+x[2]^4-x[1]^2+
            x[2]^2+3x[1]^2*x[2]^2
        dom = BoxDomain{T}(-ones(T, 2), ones(T, 2))
        true_obj = 0.814814
    elseif polyname == :robinson_ball
        DP.@polyvar x[1:2]
        f = 1+x[1]^6+x[2]^6-x[1]^4*x[2]^2+x[1]^4-x[1]^2*x[2]^4+x[2]^4-x[1]^2+
            x[2]^2+3x[1]^2*x[2]^2
        dom = BallDomain{T}(zeros(T, 2), sqrt(T(2)))
        true_obj = 0.814814
    elseif polyname == :rosenbrock
        DP.@polyvar x[1:2]
        f = (1-x[1])^2+100*(x[1]^2-x[2])^2
        dom = BoxDomain{T}(fill(-5, 2), fill(10, 2))
        true_obj = 0
    elseif polyname == :rosenbrock_ball
        DP.@polyvar x[1:2]
        f = (1-x[1])^2+100*(x[1]^2-x[2])^2
        dom = BallDomain{T}(T(2.5) * ones(T, 2), T(7.5) * sqrt(T(2)))
        true_obj = 0
    elseif polyname == :schwefel
        DP.@polyvar x[1:3]
        f = (x[1]-x[2]^2)^2+(x[2]-1)^2+(x[1]-x[3]^2)^2+(x[3]-1)^2
        dom = BoxDomain{T}(fill(-10, 3), fill(10, 3))
        true_obj = 0
    elseif polyname == :schwefel_ball
        DP.@polyvar x[1:3]
        f = (x[1]-x[2]^2)^2+(x[2]-1)^2+(x[1]-x[3]^2)^2+(x[3]-1)^2
        dom = BallDomain{T}(zeros(T, 3), 10 * sqrt(T(3)))
        true_obj = 0
    else
        error("poly $polyname not recognized")
    end

    return (x, f, dom, true_obj)
end


text:
pruneable(m) = Flux.trainable(m)
pruneable(m::Dense) = (weight = m.weight,)
pruneable(m::Conv) = (weight = m.weight,)

function mappruneable1(f, x)
    children, re = functor(x)
    tchildren = pruneable(x)
    re(map(x -> x ∈ tchildren ? f(x) : x, children))
end

  # See https://github.com/FluxML/Functors.jl/issues/2 for a discussion regarding the need for
  # cache.
function mappruneable(f, x; exclude = isleaf, cache = IdDict())
    haskey(cache, x) && return cache[x]
    y = exclude(x) ? f(x) :
                     mappruneable1(x -> mappruneable(f, x, cache = cache, exclude = exclude), x)
    cache[x] = y

    return y
end