text:
module TestCallbacks

using CPLEX
using Random
using Test

const MOI = CPLEX.MOI

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function callback_simple_model()
    model = CPLEX.Optimizer()
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.NumberOfThreads(), 1)
    MOI.set(model, MOI.RawOptimizerAttribute("CPX_PARAM_PREIND"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("CPX_PARAM_HEURFREQ"), -1)

    MOI.Utilities.loadfromstring!(
        model,
        """
    variables: x, y
    maxobjective: y
    c1: x in Integer()
    c2: y in Integer()
    c3: x in Interval(0.0, 2.5)
    c4: y in Interval(0.0, 2.5)
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    y = MOI.get(model, MOI.VariableIndex, "y")
    return model, x, y
end

function callback_knapsack_model()
    model = CPLEX.Optimizer()
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.NumberOfThreads(), 1)
    MOI.set(model, MOI.RawOptimizerAttribute("CPX_PARAM_PREIND"), 0)
    MOI.set(model, MOI.RawOptimizerAttribute("CPX_PARAM_HEURFREQ"), -1)

    N = 30
    x = MOI.add_variables(model, N)
    MOI.add_constraints(model, x, MOI.ZeroOne())
    MOI.set.(model, MOI.VariablePrimalStart(), x, 0.0)
    Random.seed!(1)
    item_weights, item_values = rand(N), rand(N)
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_weights, x), 0.0),
        MOI.LessThan(10.0),
    )
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(item_values, x), 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    return model, x, item_weights
end

function test_LazyConstraint()
    model, x, y = callback_simple_model()
    lazy_called = false
    MOI.set(
        model,
        MOI.LazyConstraintCallback(),
        cb_data -> begin
            lazy_called = true
            x_val = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), x)
            y_val = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), y)
            status = MOI.get(
                model,
                MOI.CallbackNodeStatus(cb_data),
            )::MOI.CallbackNodeStatusCode
            if round.(Int, [x_val, y_val]) ≈ [x_val, y_val]
                atol = 1e-6
                @test status == MOI.CALLBACK_NODE_STATUS_INTEGER
            else
                @test status == MOI.CALLBACK_NODE_STATUS_FRACTIONAL
            end
            @test MOI.supports(model, MOI.LazyConstraint(cb_data))
            if y_val - x_val > 1 + 1e-6
                MOI.submit(
                    model,
                    MOI.LazyConstraint(cb_data),
                    MOI.ScalarAffineFunction{Float64}(
                        MOI.ScalarAffineTerm.([-1.0, 1.0], [x, y]),
                        0.0,
                    ),
                    MOI.LessThan{Float64}(1.0),
                )
            elseif y_val + x_val > 3 + 1e-6
                MOI.submit(
                    model,
                    MOI.LazyConstraint(cb_data),
                    MOI.ScalarAffineFunction{Float64}(
                        MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
                        0.0,
                    ),
                    MOI.LessThan{Float64}(3.0),
                )
            end
        end,
    )
    @test MOI.supports(model, MOI.LazyConstraintCallback())
    MOI.optimize!(model)
    @test lazy_called
    @test MOI.get(model, MOI.VariablePrimal(), x) == 1
    @test MOI.get(model, MOI.VariablePrimal(), y) == 2
end

function test_OptimizeInProgress()
    model, x, y = callback_simple_model()
    MOI.set(
        model,
        MOI.LazyConstraintCallback(),
        cb_data -> begin
            @test_throws(
                MOI.OptimizeInProgress(MOI.VariablePrimal()),
                MOI.get(model, MOI.VariablePrimal(), x)
            )
            @test_throws(
                MOI.OptimizeInProgress(MOI.ObjectiveValue()),
                MOI.get(model, MOI.ObjectiveValue())
            )
            @test_throws(
                MOI.OptimizeInProgress(MOI.ObjectiveBound()),
                MOI.get(model, MOI.ObjectiveBound())
            )
        end,
    )
    return MOI.optimize!(model)
end

function test_LazyConstraint_UserCut()
    model, x, y = callback_simple_model()
    cb = nothing
    MOI.set(
        model,
        MOI.LazyConstraintCallback(),
        cb_data -> begin
            cb = cb_data
            MOI.submit(
                model,
                MOI.UserCut(cb_data),
                MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
                MOI.LessThan(2.0),
            )
        end,
    )
    @test_throws(
        MOI.InvalidCallbackUsage(MOI.LazyConstraintCallback(), MOI.UserCut(cb)),
        MOI.optimize!(model)
    )
end

function test_LazyConstraint_HeuristicSolution()
    model, x, y = callback_simple_model()
    cb = nothing
    MOI.set(
        model,
        MOI.LazyConstraintCallback(),
        cb_data -> begin
            cb = cb_data
            MOI.submit(model, MOI.HeuristicSolution(cb_data), [x], [2.0])
        end,
    )
    @test_throws(
        MOI.InvalidCallbackUsage(
            MOI.LazyConstraintCallback(),
            MOI.HeuristicSolution(cb),
        ),
        MOI.optimize!(model)
    )
end

function test_UserCut()
    model, x, item_weights = callback_knapsack_model()
    user_cut_submitted = false
    MOI.set(
        model,
        MOI.UserCutCallback(),
        cb_data -> begin
            terms = MOI.ScalarAffineTerm{Float64}[]
            accumulated = 0.0
            for (i, xi) in enumerate(x)
                if MOI.get(model, MOI.CallbackVariablePrimal(cb_data), xi) > 0.0
                    push!(terms, MOI.ScalarAffineTerm(1.0, xi))
                    accumulated += item_weights[i]
                end
            end
            @test MOI.supports(model, MOI.UserCut(cb_data))
            if accumulated > 10.0
                MOI.submit(
                    model,
                    MOI.UserCut(cb_data),
                    MOI.ScalarAffineFunction{Float64}(terms, 0.0),
                    MOI.LessThan{Float64}(length(terms) - 1),
                )
                user_cut_submitted = true
            end
        end,
    )
    @test MOI.supports(model, MOI.UserCutCallback())
    MOI.optimize!(model)
    @test user_cut_submitted
end

function test_UserCut_LazyConstraint()
    model, x, item_weights = callback_knapsack_model()
    cb = nothing
    MOI.set(
        model,
        MOI.UserCutCallback(),
        cb_data -> begin
            cb = cb_data
            MOI.submit(
                model,
                MOI.LazyConstraint(cb_data),
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
                MOI.LessThan(5.0),
            )
        end,
    )
    @test_throws(
        MOI.InvalidCallbackUsage(MOI.UserCutCallback(), MOI.LazyConstraint(cb)),
        MOI.optimize!(model)
    )
end

function test_UserCut_HeuristicSolution()
    model, x, item_weights = callback_knapsack_model()
    cb = nothing
    MOI.set(
        model,
        MOI.UserCutCallback(),
        cb_data -> begin
            cb = cb_data
            MOI.submit(model, MOI.HeuristicSolution(cb_data), [x[1]], [0.0])
        end,
    )
    @test_throws(
        MOI.InvalidCallbackUsage(
            MOI.UserCutCallback(),
            MOI.HeuristicSolution(cb),
        ),
        MOI.optimize!(model)
    )
end

function test_Heuristic()
    model, x, item_weights = callback_knapsack_model()
    callback_called = false
    MOI.set(
        model,
        MOI.HeuristicCallback(),
        cb_data -> begin
            x_vals = MOI.get.(model, MOI.CallbackVariablePrimal(cb_data), x)
            status = MOI.get(
                model,
                MOI.CallbackNodeStatus(cb_data),
            )::MOI.CallbackNodeStatusCode
            if round.(Int, x_vals) ≈ x_vals
                atol = 1e-6
                @test status == MOI.CALLBACK_NODE_STATUS_INTEGER
            else
                @test status == MOI.CALLBACK_NODE_STATUS_FRACTIONAL
            end
            @test MOI.supports(model, MOI.HeuristicSolution(cb_data))
            @test MOI.submit(
                model,
                MOI.HeuristicSolution(cb_data),
                x,
                floor.(x_vals),
            ) == MOI.HEURISTIC_SOLUTION_UNKNOWN
            callback_called = true
        end,
    )
    @test MOI.supports(model, MOI.HeuristicCallback())
    MOI.optimize!(model)
    @test callback_called
end

function test_Heuristic_LazyConstraint()
    model, x, item_weights = callback_knapsack_model()
    cb = nothing
    MOI.set(
        model,
        MOI.HeuristicCallback(),
        cb_data -> begin
            cb = cb_data
            MOI.submit(
                model,
                MOI.LazyConstraint(cb_data),
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
                MOI.LessThan(5.0),
            )
        end,
    )
    @test_throws(
        MOI.InvalidCallbackUsage(
            MOI.HeuristicCallback(),
            MOI.LazyConstraint(cb),
        ),
        MOI.optimize!(model)
    )
end

function test_Heuristic_UserCut()
    model, x, item_weights = callback_knapsack_model()
    cb = nothing
    MOI.set(
        model,
        MOI.HeuristicCallback(),
        cb_data -> begin
            cb = cb_data
            MOI.submit(
                model,
                MOI.UserCut(cb_data),
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
                MOI.LessThan(5.0),
            )
        end,
    )
    @test_throws(
        MOI.InvalidCallbackUsage(MOI.HeuristicCallback(), MOI.UserCut(cb)),
        MOI.optimize!(model)
    )
end

function test_CallbackFunction_OptimizeInProgress()
    model, x, y = callback_simple_model()
    MOI.set(
        model,
        CPLEX.CallbackFunction(),
        (cb_data, cb_context) -> begin
            @test_throws(
                MOI.OptimizeInProgress(MOI.VariablePrimal()),
                MOI.get(model, MOI.VariablePrimal(), x)
            )
            @test_throws(
                MOI.OptimizeInProgress(MOI.ObjectiveValue()),
                MOI.get(model, MOI.ObjectiveValue())
            )
            @test_throws(
                MOI.OptimizeInProgress(MOI.ObjectiveBound()),
                MOI.get(model, MOI.ObjectiveBound())
            )
        end,
    )
    @test MOI.supports(model, CPLEX.CallbackFunction())
    return MOI.optimize!(model)
end

function test_CallbackFunction_LazyConstraint()
    model, x, y = callback_simple_model()
    cb_calls = Clong[]
    function callback_function(
        cb_data::CPLEX.CallbackContext,
        cb_context::Clong,
    )
        push!(cb_calls, cb_context)
        if cb_context != CPX_CALLBACKCONTEXT_CANDIDATE
            return
        end
        CPLEX.load_callback_variable_primal(cb_data, cb_context)
        x_val = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), x)
        y_val = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), y)
        if y_val - x_val > 1 + 1e-6
            MOI.submit(
                model,
                MOI.LazyConstraint(cb_data),
                MOI.ScalarAffineFunction{Float64}(
                    MOI.ScalarAffineTerm.([-1.0, 1.0], [x, y]),
                    0.0,
                ),
                MOI.LessThan{Float64}(1.0),
            )
        elseif y_val + x_val > 3 + 1e-6
            MOI.submit(
                model,
                MOI.LazyConstraint(cb_data),
                MOI.ScalarAffineFunction{Float64}(
                    MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
                    0.0,
                ),
                MOI.LessThan{Float64}(3.0),
            )
        end
    end
    MOI.set(model, CPLEX.CallbackFunction(), callback_function)
    MOI.optimize!(model)
    @test MOI.get(model, MOI.VariablePrimal(), x) == 1
    @test MOI.get(model, MOI.VariablePrimal(), y) == 2
    @test length(cb_calls) > 0
end

function test_CallbackFunction_UserCut()
    model, x, item_weights = callback_knapsack_model()
    user_cut_submitted = false
    cb_calls = Clong[]
    MOI.set(
        model,
        CPLEX.CallbackFunction(),
        (cb_data, cb_context) -> begin
            push!(cb_calls, cb_context)
            if cb_context != CPX_CALLBACKCONTEXT_RELAXATION
                return
            end
            CPLEX.load_callback_variable_primal(cb_data, cb_context)
            terms = MOI.ScalarAffineTerm{Float64}[]
            accumulated = 0.0
            for (i, xi) in enumerate(x)
                if MOI.get(model, MOI.CallbackVariablePrimal(cb_data), xi) > 0.0
                    push!(terms, MOI.ScalarAffineTerm(1.0, xi))
                    accumulated += item_weights[i]
                end
            end
            if accumulated > 10.0
                MOI.submit(
                    model,
                    MOI.UserCut(cb_data),
                    MOI.ScalarAffineFunction{Float64}(terms, 0.0),
                    MOI.LessThan{Float64}(length(terms) - 1),
                )
                user_cut_submitted = true
            end
        end,
    )
    MOI.optimize!(model)
    @test user_cut_submitted
end

function test_CallbackFunction_HeuristicSolution()
    model, x, item_weights = callback_knapsack_model()
    callback_called = false
    cb_calls = Clong[]
    MOI.set(
        model,
        CPLEX.CallbackFunction(),
        (cb_data, cb_context) -> begin
            push!(cb_calls, cb_context)
            if cb_context != CPX_CALLBACKCONTEXT_RELAXATION
                return
            end
            CPLEX.load_callback_variable_primal(cb_data, cb_context)
            x_vals =
                MOI.get.(model, MOI.CallbackVariablePrimal(cb_data), x)
            @test MOI.submit(
                model,
                MOI.HeuristicSolution(cb_data),
                x,
                floor.(x_vals),
            ) == MOI.HEURISTIC_SOLUTION_UNKNOWN
            callback_called = true
        end,
    )
    MOI.optimize!(model)
    @test callback_called
end

function test_CallbackFunction_CallbackNodeStatus()
    model, x, item_weights = callback_knapsack_model()
    unknown_reached = false
    MOI.set(
        model,
        CPLEX.CallbackFunction(),
        (cb_data, cb_context) -> begin
            if MOI.get(model, MOI.CallbackNodeStatus(cb_data)) ==
               MOI.CALLBACK_NODE_STATUS_UNKNOWN
                unknown_reached = true
            end
        end,
    )
    MOI.optimize!(model)
    @test unknown_reached
end

function test_CPXcallbackabort()
    model = CPLEX.Optimizer()
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.Integer())
    MOI.set(model, MOI.NumberOfThreads(), 1)
    MOI.set(
        model,
        CPLEX.CallbackFunction(),
        (cb_data, context_id) -> begin
            @show context_id
            CPXcallbackabort(cb_data)
        end,
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INTERRUPTED
end

"""
    test_InterruptException()

This test simulates an InterruptException being thrown. It is a little
complicated due to the delayed handling of `terminate`, which _schedules_ a
request for termination, rather than terminating immediately. This means CPLEX
may continue to call the callback after the interruption.

First, we must ensure that InterruptException() is only thrown once. Double
interrupting would interrupt our handling of the first interrupt!

Second, if the model is too simplisitic, CPLEX may be able to prove optimality
after we have interrupted, but before it has decided to actually exit the solve.
"""
function test_InterruptException()
    model = CPLEX.Optimizer()
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.Integer())
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    MOI.set(model, MOI.NumberOfThreads(), 1)
    i = 0.0
    interrupt_thrown = false
    MOI.set(
        model,
        CPLEX.CallbackFunction(),
        (cb_data, cb_where) -> begin
            if cb_where != CPLEX.CPX_CALLBACKCONTEXT_CANDIDATE
                return
            end
            MOI.submit(
                model,
                MOI.LazyConstraint(cb_data),
                MOI.ScalarAffineFunction{Float64}(
                    [MOI.ScalarAffineTerm(1.0, x)],
                    0.0,
                ),
                MOI.GreaterThan{Float64}(i),
            )
            i += 1
            if !interrupt_thrown
                interrupt_thrown = true
                throw(InterruptException())
            end
        end,
    )
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.INTERRUPTED
end

function test_CallbackFunction_broadcast()
    model, x, _ = callback_knapsack_model()
    f(cb_data, x) = MOI.get(model, MOI.CallbackVariablePrimal(cb_data), x)
    solutions = Vector{Float64}[]
    MOI.set(
        model,
        CPLEX.CallbackFunction(),
        (cb_data, cb_where) -> begin
            if cb_where == CPLEX.CPX_CALLBACKCONTEXT_CANDIDATE
                CPLEX.load_callback_variable_primal(cb_data, cb_where)
                push!(solutions, f.(cb_data, x))
            end
        end,
    )
    MOI.optimize!(model)
    @test length(solutions) > 0
    @test length(solutions[1]) == length(x)
end

end  # module TestCallbacks

TestCallbacks.runtests()


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

# check that calling `plot` in the following format
#   (::Array{Any,1}, ::Gadfly.#plot, ::Gadfly.Geom.#density,
#          ::Gadfly.Guide.XLabel, ::Gadfly.Guide.YLabel)
# does not raise an ambiguity MethodError on Julia v0.5+
plot(x=1:10, y=1:10, color=1:10, Geom.density,
                   Guide.xlabel("foo"), Guide.ylabel("bar"))

ext: jl


text:
rng = StableRNG(51803)

const Vec = AbstractVector

@testset "misclassification_rate" begin
    y    = categorical(collect("asdfasdfaaassdd"))
    yhat = categorical(collect("asdfaadfaasssdf"))
    w = 1:15
    ym = vcat(y, [missing,])
    yhatm = vcat(yhat, [missing,])
    wm = 1:16
    @test misclassification_rate(yhat, y) ≈ 0.2
    @test misclassification_rate(yhatm, ym) ≈ 0.2
    @test misclassification_rate(yhat, y, w) ≈ (6*1 + 11*1 + 15*1) / 15
    @test misclassification_rate(yhatm, ym, wm) ≈ (6*1 + 11*1 + 15*1) / 15
end

@testset "mcr, acc, bacc, mcc" begin
    y = categorical(['m', 'f', 'n', 'f', 'm', 'n', 'n', 'm', 'f'])
    ŷ = categorical(['f', 'f', 'm', 'f', 'n', 'm', 'n', 'm', 'f'])
    @test accuracy(ŷ, y) == 1-mcr(ŷ,y) ==
        accuracy(MLJBase._confmat(ŷ, y, warn=false))  ==
        1-mcr(MLJBase._confmat(ŷ, y, warn=false))
    w = randn(rng,length(y))
    @test accuracy(ŷ, y, w) == 1-mcr(ŷ,y,w)

    ## balanced accuracy
    y = categorical([
        3, 4, 1, 1, 1, 4, 1, 3, 3, 1, 2, 3, 1, 3, 3, 3, 2, 4, 3, 2, 1, 3,
        3, 1, 1, 1, 2, 4, 1, 4, 4, 4, 1, 1, 4, 4, 3, 1, 2, 2, 3, 4, 2, 1,
        2, 2, 3, 2, 2, 3, 1, 2, 3, 4, 1, 2, 4, 2, 1, 4, 3, 2, 3, 3, 3, 1,
        3, 1, 4, 3, 1, 2, 3, 1, 2, 2, 4, 4, 1, 3, 2, 1, 4, 3, 3, 1, 3, 1,
        2, 2, 2, 2, 2, 3, 2, 1, 1, 4, 2, 2])
    ŷ = categorical([
        2, 3, 2, 1, 2, 2, 3, 3, 2, 4, 2, 3, 2, 4, 3, 4, 4, 2, 1, 3, 3, 3,
        3, 3, 2, 4, 4, 3, 4, 4, 1, 2, 3, 2, 4, 1, 2, 3, 1, 4, 2, 2, 1, 2,
        3, 2, 2, 4, 3, 2, 2, 2, 1, 2, 2, 1, 3, 1, 4, 1, 2, 1, 2, 4, 3, 2,
        4, 3, 2, 4, 4, 2, 4, 3, 2, 3, 1, 2, 1, 2, 1, 2, 3, 1, 1, 3, 4, 2,
        4, 4, 2, 1, 3, 2, 2, 4, 1, 1, 4, 1])
    w = [
        0.5, 1.4, 0.6, 1. , 0.1, 0.5, 1.2, 0.2, 1.8, 0.3, 0.6, 2.2, 0.1,
        1.4, 0.2, 0.4, 0.6, 2.1, 0.7, 0.2, 0.9, 0.4, 0.7, 0.3, 0.1, 1.7,
        0.2, 0.7, 1.2, 1. , 0.9, 0.4, 0.5, 0.5, 0.5, 1. , 0.3, 0.1, 0.2,
        0. , 2.2, 0.8, 0.9, 0.8, 1.3, 0.2, 0.4, 0.7, 1. , 0.7, 1.7, 0.7,
        1.1, 1.8, 0.1, 1.2, 1.8, 1. , 0.1, 0.5, 0.6, 0.7, 0.6, 1.2, 0.6,
        1.2, 0.5, 0.5, 0.8, 0.2, 0.6, 1. , 0.3, 1. , 0.2, 1.1, 1.1, 1.1,
        0.6, 1.4, 1.2, 0.3, 1.1, 0.2, 0.5, 1.6, 0.3, 1. , 0.3, 0.9, 0.9,
        0. , 0.6, 0.6, 0.4, 0.5, 0.4, 0.2, 0.9, 0.4]
    sk_bacc = 0.17493386243386244 # note: sk-learn reverses ŷ and y
    @test bacc(ŷ, y) ≈ sk_bacc
    sk_adjusted_bacc =  -0.10008818342151675
    @test BalancedAccuracy(adjusted=true)(ŷ, y) ≈ sk_adjusted_bacc
    sk_bacc_w = 0.1581913163016446
    @test bacc(ŷ, y, w) ≈ sk_bacc_w
    sk_adjusted_bacc_w = -0.1224115782644738
    @test BalancedAccuracy(adjusted=true)(ŷ, y, w) ≈ sk_adjusted_bacc_w

    ## matthews correlation
    sk_mcc = -0.09759509982785947
    @test mcc(ŷ, y) == matthews_correlation(ŷ, y) ≈ sk_mcc
    # invariance with respect to permutation ?
    cm = MLJBase._confmat(ŷ, y, perm=[3, 1, 2, 4])
    @test mcc(cm) ≈ sk_mcc

    # Issue #381
    cm = MLJBase.ConfusionMatrixObject([29488 13017; 12790 29753], ["0.0", "1.0"])
    @test mcc(cm) ≈ 0.39312321239417797
end

@testset "kappa" begin
    # Binary case
    y_b = categorical([2, 2, 2, 1, 1, 1, 1, 1, 2, 2, 1, 2, 2, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 2, 2, 1, 1, 1, 2, 2])
    ŷ_b = categorical([1, 1, 2, 2, 2, 2, 1, 1, 2, 1, 1, 2, 1, 2, 2, 2, 1, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 2, 2])
    cm_b = MLJBase._confmat(y_b, ŷ_b, warn=false)
    p0_b = (4+10)/30
    pe_b = (13*11 + 17*19)/(30*30)

    # Multiclass case
    y_m = categorical([5, 5, 3, 5, 4, 4, 2, 2, 3, 2, 5, 2, 4, 3, 2, 1, 1, 5, 1, 4, 2, 5, 4, 5, 2, 3, 3, 4, 2, 4])
    ŷ_m = categorical([1, 1, 1, 5, 4, 2, 1, 3, 4, 4, 2, 5, 4, 4, 1, 5, 5, 2, 3, 3, 1, 3, 2, 5, 5, 2, 3, 2, 5, 3])
    cm_m = MLJBase._confmat(ŷ_m, y_m, warn=false)
    p0_m = 5/30
    pe_m = (3*6 + 8*6 + 5*6 + 7*5 + 7*7)/(30*30)

    # Tests
    @test kappa(y_m, ŷ_m) ≈ (p0_m - pe_m)/(1 - pe_m)
    @test kappa(y_b, ŷ_b) ≈ (p0_b - pe_b)/(1 - pe_b)
    @test kappa(cm_m)     == kappa(y_m, ŷ_m)
    @test kappa(cm_b)     == kappa(y_b, ŷ_b)
    @test kappa(ŷ_m, y_m) == kappa(y_m, ŷ_m)
    @test kappa(ŷ_b, y_b) == kappa(y_b, ŷ_b)
    @test kappa(y_m, y_m) == 1.0
    @test kappa(y_b, y_b) == 1.0
end

@testset "confusion matrix {2}" begin
    # first class is 1 is assumed negative, second positive
    y = categorical([1, 2, 1, 2, 1, 1, 2])
    ŷ = categorical([1, 2, 2, 2, 2, 1, 2])
    cm = MLJBase._confmat(ŷ, y, warn=false)
    TN = sum(ŷ .== y .== 1) # pred and true = - (1)
    TP = sum(ŷ .== y .== 2) # pred and true = + (2)
    FP = sum(ŷ .!= y .== 1) # pred + (2) and true - (1)
    FN = sum(ŷ .!= y .== 2) # pred - (1) and true + (2)
    @test cm[1,1] == TN
    @test cm[2,2] == TP
    @test cm[1,2] == FN
    @test cm[2,1] == FP

    ym = categorical([1, missing, 2, 1, 2, 1, 1, 1, 2])
    ŷm = categorical([1, 2,       2, 2, 2, missing, 2, 1, 2])
    cm = MLJBase._confmat(ŷ, y, warn=false)
    TN = sum(skipmissing(ŷ .== y .== 1)) # pred and true = - (1)
    TP = sum(skipmissing(ŷ .== y .== 2)) # pred and true = + (2)
    FP = sum(skipmissing(ŷ .!= y .== 1)) # pred + (2) and true - (1)
    FN = sum(skipmissing(ŷ .!= y .== 2)) # pred - (1) and true + (2)
    @test cm[1,1] == TN
    @test cm[2,2] == TP
    @test cm[1,2] == FN
    @test cm[2,1] == FP

    cm2 = MLJBase._confmat(ŷ, y; rev=true)
    @test cm2[1,1] == cm[2,2]
    @test cm2[1,2] == cm[2,1]
    @test cm2[2,2] == cm[1,1]
    @test cm2[2,1] == cm[1,2]

    @test accuracy(ŷ, y) == accuracy(cm) == sum(y .== ŷ) / length(y)

    @test @test_logs((:warn, r"The classes are un-ordered"),
                     recall(ŷ, y) == TP / (TP + FN))

    ŷ = coerce(ŷ, Union{Missing,OrderedFactor})
    y = coerce(y, Union{Missing,OrderedFactor})

    @test precision(ŷ, y)   == TP / (TP + FP)
    @test specificity(ŷ, y) == TN / (TN + FP)
    @test f1score(ŷ, y) ≈
        2.0 / (1.0 / recall(ŷ, y) + 1.0 / precision(ŷ, y))

    recall_rev = Recall(rev=true)
    @test recall_rev(ŷ, y) ==
        TN / (TN + FP) # no warning because rev is specified
    precision_rev = Precision(rev=true)
    @test precision_rev(ŷ, y) == TN / (TN + FN)
    specificity_rev = Specificity(rev=true)
    @test specificity_rev(ŷ, y) == TP / (TP + FN)
    f1score_rev = FScore(rev=true)
    @test f1score_rev(ŷ, y) ≈
        2.0 / (1.0 / recall_rev(ŷ, y) + 1.0 / precision_rev(ŷ, y))
end

@testset "confusion matrix {n}" begin
    y = coerce([1, 2, 0, 2, 1, 0, 0, 1, 2, 2, 2, 1, 2,
                            2, 1, 0, 1, 1, 1, 2, 1, 2, 2, 1, 2, 1,
                            2, 2, 2], Multiclass)
    ŷ = coerce([2, 0, 2, 2, 2, 0, 1, 2, 1, 2, 0, 1, 2,
                            1, 1, 1, 2, 0, 1, 2, 1, 2, 2, 2, 1, 2,
                            1, 2, 2], Multiclass)
    class_w = Dict(0=>0,2=>2,1=>1)
    cm = MLJBase._confmat(ŷ, y, warn=false)

    #               ┌─────────────────────────────────────────┐
    #               │              Ground Truth               │
    # ┌─────────────┼─────────────┬─────────────┬─────────────┤
    # │  Predicted  │      0      │      1      │      2      │
    # ├─────────────┼─────────────┼─────────────┼─────────────┤
    # │      0      │      1      │      1      │      2      │
    # ├─────────────┼─────────────┼─────────────┼─────────────┤
    # │      1      │      2      │      4      │      4      │
    # ├─────────────┼─────────────┼─────────────┼─────────────┤
    # │      2      │      1      │      6      │      8      │
    # └─────────────┴─────────────┴─────────────┴─────────────┘

    cm_tp   = [1; 4; 8]
    cm_tn   = [22; 12; 8]
    cm_fp   = [1+2; 2+4; 1+6]
    cm_fn   = [2+1; 1+6; 2+4]
    cm_prec = cm_tp ./ ( cm_tp + cm_fp  )
    cm_rec  = cm_tp ./ ( cm_tp + cm_fn  )

    # Check if is positive
    m = MulticlassTruePositive(;return_type=Vector)
    @test  [0; 0; 0] <= m(ŷ, y) == cm_tp
    m = MulticlassTrueNegative(;return_type=Vector)
    @test  [0; 0; 0] <= m(ŷ, y) == cm_tn
    m = MulticlassFalsePositive(;return_type=Vector)
    @test  [0; 0; 0] <= m(ŷ, y) == cm_fp
    m = MulticlassFalseNegative(;return_type=Vector)
    @test  [0; 0; 0] <= m(ŷ, y) == cm_fn

    # Check if is in [0,1]
    m = MulticlassTruePositiveRate(average=no_avg;return_type=Vector)
    @test  [0; 0; 0] <= m(ŷ, y) == cm_tp ./ (cm_fn.+cm_tp) <= [1; 1; 1]
    m = MulticlassTrueNegativeRate(average=no_avg;return_type=Vector)
    @test  [0; 0; 0] <= m(ŷ, y) == cm_tn ./ (cm_tn.+cm_fp) <= [1; 1; 1]
    m = MulticlassFalsePositiveRate(average=no_avg;return_type=Vector)
    @test  [0; 0; 0] <= m(ŷ, y) == 1 .- cm_tn ./ (cm_tn.+cm_fp) <= [1; 1; 1]
    m = MulticlassFalseNegativeRate(average=no_avg;return_type=Vector)
    @test  [0; 0; 0] <= m(ŷ, y) == 1 .- cm_tp ./ (cm_fn.+cm_tp) <= [1; 1; 1]

    #`no_avg` and `LittleDict`
    @test collect(values(MulticlassPrecision(average=no_avg)(cm))) ≈
        collect(values(MulticlassPrecision(average=no_avg)(ŷ, y))) ≈
        cm_prec
    @test MulticlassPrecision(average=macro_avg)(cm) ≈
        MulticlassPrecision(average=macro_avg)(ŷ, y) ≈ mean(cm_prec)
    @test collect(keys(MulticlassPrecision(average=no_avg)(cm)))  ==
        collect(keys(MulticlassPrecision(average=no_avg)(ŷ, y))) ==
        ["0"; "1"; "2"]
    @test collect(values(MulticlassRecall(average=no_avg)(cm))) ≈
        collect(values(MulticlassRecall(average=no_avg)(ŷ, y))) ≈
        cm_rec
    @test collect(values(MulticlassFScore(average=no_avg)(cm))) ≈
        collect(values(MulticlassFScore(average=no_avg)(ŷ, y))) ≈
        2 ./ ( 1 ./ cm_prec + 1 ./ cm_rec )

    #`no_avg` and `LittleDict` with class weights
    @test collect(values(MulticlassPrecision(average=no_avg)(cm, class_w))) ≈
        collect(values(MulticlassPrecision(average=no_avg)(ŷ, y, class_w))) ≈
        cm_prec .* [0; 1; 2]
    @test collect(values(MulticlassRecall(average=no_avg)(cm, class_w))) ≈
        collect(values(MulticlassRecall(average=no_avg)(ŷ, y, class_w))) ≈
        cm_rec .* [0; 1; 2]
    @test collect(values(MulticlassFScore(average=no_avg)(cm, class_w))) ≈
        collect(values(MulticlassFScore(average=no_avg)(ŷ, y, class_w))) ≈
        2 ./ ( 1 ./ cm_prec + 1 ./ cm_rec ) .* [0; 1; 2]

    #`macro_avg` and `LittleDict`
    macro_prec = MulticlassPrecision(average=macro_avg)
    macro_rec  = MulticlassRecall(average=macro_avg)

    @test macro_prec(cm)    ≈ macro_prec(ŷ, y)    ≈ mean(cm_prec)
    @test macro_rec(cm)     ≈ macro_rec(ŷ, y)     ≈ mean(cm_rec)
    @test macro_f1score(cm) ≈ macro_f1score(ŷ, y) ≈ mean(2 ./ ( 1 ./ cm_prec + 1 ./ cm_rec ))

    #`micro_avg` and `LittleDict`
    micro_prec = MulticlassPrecision(average=micro_avg)
    micro_rec  = MulticlassRecall(average=micro_avg)

    @test micro_prec(cm)    == micro_prec(ŷ, y)    == sum(cm_tp) ./ sum(cm_fp.+cm_tp)
    @test micro_rec(cm)     == micro_rec(ŷ, y)     == sum(cm_tp) ./ sum(cm_fn.+cm_tp)
    @test micro_f1score(cm) == micro_f1score(ŷ, y) ==
    2 ./ ( 1 ./ ( sum(cm_tp) ./ sum(cm_fp.+cm_tp) ) + 1 ./ ( sum(cm_tp) ./ sum(cm_fn.+cm_tp) ) )

    #`no_avg` and `Vector` with class weights
    vec_precision = MulticlassPrecision(return_type=Vector)
    vec_recall    = MulticlassRecall(return_type=Vector)
    vec_f1score   = MulticlassFScore(return_type=Vector)

    @test vec_precision(cm, class_w) ≈ vec_precision(ŷ, y, class_w) ≈
        mean(cm_prec .* [0; 1; 2])
    @test vec_recall(cm, class_w)    ≈ vec_recall(ŷ, y, class_w)    ≈
        mean(cm_rec .* [0; 1; 2])
    @test vec_f1score(cm, class_w)   ≈ vec_f1score(ŷ, y, class_w)   ≈
        mean(2 ./ ( 1 ./ cm_prec + 1 ./ cm_rec ) .* [0; 1; 2])

    #`macro_avg` and `Vector`
    v_ma_prec = MulticlassPrecision(average=macro_avg,
                                    return_type=Vector)
    v_ma_rec  = MulticlassRecall(average=macro_avg, return_type=Vector)
    v_ma_f1   = MulticlassFScore(average=macro_avg, return_type=Vector)

    @test v_ma_prec(cm) ≈ v_ma_prec(ŷ, y) ≈ mean(cm_prec)
    @test v_ma_rec(cm)  ≈ v_ma_rec(ŷ, y)  ≈ mean(cm_rec)
    @test v_ma_f1(cm)   ≈ v_ma_f1(ŷ, y)   ≈ mean(2 ./ ( 1 ./ cm_prec + 1 ./ cm_rec ))

    #`macro_avg` and `Vector` with class weights
    @test v_ma_prec(cm, class_w) ≈ v_ma_prec(ŷ, y, class_w) ≈
        mean(cm_prec .* [0, 1, 2])
    @test v_ma_rec(cm, class_w)  ≈ v_ma_rec(ŷ, y, class_w)  ≈
        mean(cm_rec .* [0, 1, 2])
    @test v_ma_f1(cm, class_w)   ≈ v_ma_f1(ŷ, y, class_w)   ≈
        mean(2 ./ ( 1 ./ cm_prec + 1 ./ cm_rec ) .* [0, 1, 2])

    #`micro_avg` and `Vector`
    v_mi_prec = MulticlassPrecision(average=micro_avg, return_type=Vector)
    v_mi_rec  = MulticlassRecall(average=micro_avg, return_type=Vector)
    v_mi_f1   = MulticlassFScore(average=micro_avg, return_type=Vector)

    @test v_mi_prec(cm) == v_mi_prec(ŷ, y) == sum(cm_tp) ./ sum(cm_fp.+cm_tp)
    @test v_mi_rec(cm)  == v_mi_rec(ŷ, y)  == sum(cm_tp) ./ sum(cm_fn.+cm_tp)
    @test v_mi_f1(cm)   == v_mi_f1(ŷ, y)   ==
    2 ./ ( 1 ./ ( sum(cm_tp) ./ sum(cm_fp.+cm_tp) ) + 1 ./ ( sum(cm_tp) ./ sum(cm_fn.+cm_tp) ) )
end

@testset "issue #630" begin
    # multiclass fscore corner case of absent class

    y = coerce([1, 2, 2, 2, 3], OrderedFactor)[1:4]
    # [1, 2, 2, 2] # but 3 is in the pool
    yhat = reverse(y)
    # [2, 2, 2, 1]

    # In this case, assigning "3" as "positive" gives all true negative,
    # and so NaN for that class's contribution to the average F1Score,
    # which should accordingly be skipped.

    # postive class | TP | FP | FN | score for that class
    # --------------|----|----|----|---------------------
    #  1            | 0  | 1  | 2  | 0
    #  2            | 2  | 1  | 1  | 2/3
    #  3            | 0  | 0  | 0  | NaN

    # mean score with skippin NaN is 1/3
    @test MulticlassFScore()(yhat, y) ≈ 1/3
end

@testset "Metadata binary" begin
    for m in (accuracy, recall, Precision(), f1score, specificity)
        e = info(m)
        m == accuracy    && (@test e.name == "Accuracy")
        m == recall      && (@test e.name == "TruePositiveRate")
        m isa Precision  && (@test e.name == "Precision")
        m == f1score     && (@test e.name == "FScore")
        m == specificity && (@test e.name == "TrueNegativeRate")
        @test e.target_scitype <: AbstractArray{<:Union{Missing,Finite}}
        @test e.prediction_type == :deterministic
        @test e.orientation == :score
        @test e.reports_each_observation == false
        @test e.is_feature_dependent == false
        if m == accuracy
            @test e.supports_weights
        else
            @test !e.supports_weights
        end
    end
    e = info(auc)
    @test e.name == "AreaUnderCurve"
    @test e.target_scitype ==
        Union{AbstractArray{<:Union{Missing,Multiclass{2}}},
              AbstractArray{<:Union{Missing,OrderedFactor{2}}}}
    @test e.prediction_type == :probabilistic
    @test e.reports_each_observation == false
    @test e.is_feature_dependent == false
    @test e.supports_weights == false
end

@testset "Metadata multiclass" begin
    for m in (MulticlassRecall(), MulticlassPrecision(),
              MulticlassFScore(), MulticlassTrueNegativeRate())
        e = info(m)
        m isa MulticlassRecall &&
            (@test e.name == "MulticlassTruePositiveRate")
        m isa MulticlassPrecision   &&
            (@test e.name == "MulticlassPrecision")
        m isa MulticlassFScore &&
            (@test e.name == "MulticlassFScore")
        m isa MulticlassTrueNegativeRate &&
            (@test e.name == "MulticlassTrueNegativeRate")
        @test e.target_scitype <: AbstractArray{<:Union{Missing,Finite}}
        @test e.prediction_type == :deterministic
        @test e.orientation == :score
        @test e.reports_each_observation == false
        @test e.is_feature_dependent == false
        @test e.supports_weights == false
        @test e.supports_class_weights == true
    end
end

@testset "More binary metrics" begin
    y = coerce([missing, 1, 2, 1, 2, 1, 1, 2, 1, 2, 2, 2, 1, 2,
                2, 1, 2, 1, 1, 1, 2, 1, 2, 2, 1, 2, 1,
                2, 2, 2, 1], Union{Missing,OrderedFactor})
    ŷ = coerce([1, 1, 2, 2, 2, 2, 1, 2, 2, 1, 2, 2, 1, 2,
                1, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 1, 2,
                1, 2, 2, missing], Union{Missing,OrderedFactor})

    # check all constructors
    m = TruePositive()
    @test m(ŷ, y) == truepositive(ŷ, y)
    m = TruePositive(rev=true)
    @test m(ŷ, y) == truenegative(ŷ, y)
    m = TrueNegative()
    @test m(ŷ, y) == truenegative(ŷ, y)
    m = FalsePositive()
    @test m(ŷ, y) == falsepositive(ŷ, y)
    m = FalseNegative()
    @test m(ŷ, y) == falsenegative(ŷ, y)
    m = TruePositiveRate()
    @test m(ŷ, y) == tpr(ŷ, y) == truepositive_rate(ŷ, y)
    m = TrueNegativeRate()
    @test m(ŷ, y) == tnr(ŷ, y) == truenegative_rate(ŷ, y)
    m = FalsePositiveRate()
    @test m(ŷ, y) == fpr(ŷ, y) == falsepositive_rate(ŷ, y)
    m = FalseNegativeRate()
    @test m(ŷ, y) == fnr(ŷ, y) == falsenegative_rate(ŷ, y)
    m = FalseDiscoveryRate()
    @test m(ŷ, y) == fdr(ŷ, y) == falsediscovery_rate(ŷ, y)
    m = Precision()
    @test m(ŷ, y) == precision(ŷ, y)
    m = NPV()
    @test m(ŷ, y) == npv(ŷ, y)
    m = FScore()
    @test m(ŷ, y) == f1score(ŷ, y)
    # check synonyms
    m = TPR()
    @test m(ŷ, y) == tpr(ŷ, y)
    m = TNR()
    @test m(ŷ, y) == tnr(ŷ, y)
    m = FPR()
    @test m(ŷ, y) == fpr(ŷ, y) == fallout(ŷ, y)
    m = FNR()
    @test m(ŷ, y) == fnr(ŷ, y) == miss_rate(ŷ, y)
    m = FDR()
    @test m(ŷ, y) == fdr(ŷ, y)
    m = PPV()
    @test m(ŷ, y) == precision(ŷ, y) == ppv(ŷ, y)
    m = Recall()
    @test m(ŷ, y) == tpr(ŷ, y) == recall(ŷ, y) ==
        sensitivity(ŷ, y) == hit_rate(ŷ, y)
    m = Specificity()
    @test m(ŷ, y) == tnr(ŷ, y) == specificity(ŷ, y) == selectivity(ŷ, y)
    # 'higher order'
    m = BACC()
    @test m(ŷ, y) == bacc(ŷ, y) == (tpr(ŷ, y) + tnr(ŷ, y))/2

    ### External comparisons
    sk_prec = 0.6111111111111112 # m.precision_score(y, yhat, pos_label=2)
    @test precision(ŷ, y) ≈ sk_prec
    sk_rec = 0.6875
    @test recall(ŷ, y) == sk_rec # m.recall_score(y, yhat, pos_label=2)
    sk_f05 = 0.625
    f05 = FScore(β=0.5)
    @test f05(ŷ, y) ≈ sk_f05 # m.fbeta_score(y, yhat, 0.5, pos_label=2)

    # reversion mechanism
    sk_prec_rev = 0.5454545454545454
    prec_rev = Precision(rev=true)
    @test prec_rev(ŷ, y) ≈ sk_prec_rev
    sk_rec_rev = 0.46153846153846156
    rec_rev = Recall(rev=true)
    @test rec_rev(ŷ, y) ≈ sk_rec_rev
end

@testset "More multiclass metrics" begin
    y = coerce(categorical([missing, 1, 2, 0, 2, 1, 0, 0, 1, 2, 2, 2, 1, 2,
                            2, 1, 0, 1, 1, 1, 2, 1, 2, 2, 1, 2, 1,
                            2, 2, 2, 0]), Union{Missing,Multiclass})
    ŷ = coerce(categorical([0, 2, 0, 2, 2, 2, 0, 1, 2, 1, 2, 0, 1, 2,
                            1, 1, 1, 2, 0, 1, 2, 1, 2, 2, 2, 1, 2,
                            1, 2, 2, missing]), Union{Missing,Multiclass})
    w = Dict(0=>1, 1=>2, 2=>3) #class_w
    # check all constructors
    m = MulticlassTruePositive()
    @test m(ŷ, y) == multiclass_truepositive(ŷ, y)
    m = MulticlassTrueNegative()
    @test m(ŷ, y) == multiclass_truenegative(ŷ, y)
    m = MulticlassFalsePositive()
    @test m(ŷ, y) == multiclass_falsepositive(ŷ, y)
    m = MulticlassFalseNegative()
    @test m(ŷ, y) == multiclass_falsenegative(ŷ, y)
    m = MulticlassTruePositiveRate()
    @test m(ŷ, y) == multiclass_tpr(ŷ, y) ==
        multiclass_truepositive_rate(ŷ, y)
    @test m(ŷ, y, w) == multiclass_tpr(ŷ, y, w) ==
        multiclass_truepositive_rate(ŷ, y, w)
    m = MulticlassTrueNegativeRate()
    @test m(ŷ, y) == multiclass_tnr(ŷ, y) ==
        multiclass_truenegative_rate(ŷ, y)
    @test m(ŷ, y, w) == multiclass_tnr(ŷ, y, w) ==
        multiclass_truenegative_rate(ŷ, y, w)
    m = MulticlassFalsePositiveRate()
    @test m(ŷ, y) == multiclass_fpr(ŷ, y) ==
        multiclass_falsepositive_rate(ŷ, y)
    @test m(ŷ, y, w) == multiclass_fpr(ŷ, y, w) ==
        multiclass_falsepositive_rate(ŷ, y, w)
    m = MulticlassFalseNegativeRate()
    @test m(ŷ, y) == multiclass_fnr(ŷ, y) ==
        multiclass_falsenegative_rate(ŷ, y)
    @test m(ŷ, y, w) == multiclass_fnr(ŷ, y, w) ==
        multiclass_falsenegative_rate(ŷ, y, w)
    m = MulticlassFalseDiscoveryRate()
    @test m(ŷ, y) == multiclass_fdr(ŷ, y) ==
        multiclass_falsediscovery_rate(ŷ, y)
    @test m(ŷ, y, w) == multiclass_fdr(ŷ, y, w) ==
        multiclass_falsediscovery_rate(ŷ, y, w)
    m = MulticlassPrecision()
    @test m(ŷ, y) == multiclass_precision(ŷ, y)
    @test m(ŷ, y, w) == multiclass_precision(ŷ, y, w)
    m = MulticlassNegativePredictiveValue()
    @test m(ŷ, y) == multiclass_npv(ŷ, y)
    @test m(ŷ, y, w) == multiclass_npv(ŷ, y, w)
    m = MulticlassFScore()
    @test m(ŷ, y) == macro_f1score(ŷ, y)
    @test m(ŷ, y, w) == macro_f1score(ŷ, y, w)
    # check synonyms
    m = MTPR(return_type=Vector)
    @test m(ŷ, y) == multiclass_tpr(ŷ, y)
    @test m(ŷ, y, w) == multiclass_tpr(ŷ, y, w)
    m = MTNR(return_type=Vector)
    @test m(ŷ, y) == multiclass_tnr(ŷ, y)
    @test m(ŷ, y, w) == multiclass_tnr(ŷ, y, w)
    m = MFPR()
    @test m(ŷ, y) == multiclass_fpr(ŷ, y) == multiclass_fallout(ŷ, y)
    @test m(ŷ, y, w) == multiclass_fpr(ŷ, y, w) ==
        multiclass_fallout(ŷ, y, w)
    m = MFNR()
    @test m(ŷ, y) == multiclass_fnr(ŷ, y) ==
        multiclass_miss_rate(ŷ, y)
    @test m(ŷ, y, w) == multiclass_fnr(ŷ, y, w) ==
        multiclass_miss_rate(ŷ, y, w)
    m = MFDR()
    @test m(ŷ, y) == multiclass_fdr(ŷ, y)
    @test m(ŷ, y, w) == multiclass_fdr(ŷ, y, w)
    m = MPPV()
    @test m(ŷ, y) == MulticlassPrecision()(ŷ, y) ==
        multiclass_ppv(ŷ, y)
    @test m(ŷ, y, w) == MulticlassPrecision()(ŷ, y, w) ==
        multiclass_ppv(ŷ, y, w)
    m = MulticlassRecall()
    @test m(ŷ, y) == multiclass_tpr(ŷ, y)
    @test m(ŷ, y, w) == multiclass_tpr(ŷ, y, w)
    @test m(ŷ, y) == multiclass_sensitivity(ŷ, y) ==
        multiclass_hit_rate(ŷ, y)
    @test m(ŷ, y, w) == multiclass_sensitivity(ŷ, y, w) ==
        multiclass_hit_rate(ŷ, y, w)
    m = MulticlassSpecificity()
    @test m(ŷ, y) == multiclass_tnr(ŷ, y) == multiclass_specificity(ŷ, y) ==
        multiclass_selectivity(ŷ, y)
    @test m(ŷ, y, w) == multiclass_tnr(ŷ, y, w) ==
        multiclass_specificity(ŷ, y, w) == multiclass_selectivity(ŷ, y, w)
end

@testset "Additional multiclass tests" begin
    table = reshape(collect("aabbbccccddbabccbacccd"), 11, 2)
    table = coerce(table, Multiclass);
    yhat = table[:,1] # ['a', 'a', 'b', 'b', 'b', 'c', 'c', 'c', 'c', 'd', 'd']
    y    = table[:,2] # ['b', 'a', 'b', 'c', 'c', 'b', 'a', 'c', 'c', 'c', 'd']
    class_w = Dict('a'=>7, 'b'=>5, 'c'=>2, 'd'=> 0)

    # class | TP | FP | TP + FP | precision | FN | TP + FN | recall
    # ------|----|----|------------------------------------|-------
    # a     | 1  | 1  | 2       | 1/2       | 1  | 2       | 1/2
    # b     | 1  | 2  | 3       | 1/3       | 2  | 3       | 1/3
    # c     | 2  | 2  | 4       | 1/2       | 3  | 5       | 2/5
    # d     | 1  | 1  | 2       | 1/2       | 0  | 1       | 1

    # helper:
    inverse(x) = 1/x
    harmonic_mean(x, y; beta=1.0) =
        (1 + inverse(beta^2))*inverse(mean(inverse(beta^2*x)+ inverse(y)))

    # precision:
    p_macro = mean([1/2, 1/3, 1/2, 1/2])
    @test MulticlassPrecision()(yhat, y) ≈ p_macro
    p_macro_w = mean([7/2, 5/3, 2/2, 0/2])
    @test MulticlassPrecision()(yhat, y, class_w) ≈ p_macro_w
    @test p_macro_w ≈
        @test_logs((:warn, r"Using macro"),
                     MulticlassPrecision(average=micro_avg)(yhat, y, class_w))
    p_micro = (1 + 1 + 2 + 1)/(2 + 3 + 4 + 2)
    @test MulticlassPrecision(average=micro_avg)(yhat, y) ≈ p_micro

    # recall:
    r_macro = mean([1/2, 1/3, 2/5, 1])
    @test MulticlassRecall(average=macro_avg)(yhat, y) ≈ r_macro
    r_macro_w = mean([7/2, 5/3, 4/5, 0/1])
    @test MulticlassRecall(average=macro_avg)(yhat, y, class_w) ≈ r_macro_w
    @test r_macro_w ≈
        @test_logs((:warn, r"Using macro"),
                     MulticlassRecall(average=micro_avg)(yhat, y, class_w))
    r_micro = (1 + 1 + 2 + 1)/(2 + 3 + 5 + 1)
    @test MulticlassPrecision(average=micro_avg)(yhat, y) ≈ r_micro

    # fscore:
    harm_means = [harmonic_mean(1/2, 1/2),
                     harmonic_mean(1/3, 1/3),
                     harmonic_mean(1/2, 2/5),
                     harmonic_mean(1/2, 1)]
    f1_macro = mean(harm_means)
    @test MulticlassFScore(average=macro_avg)(yhat, y) ≈ f1_macro
    @test MulticlassFScore(average=no_avg,
                           return_type=Vector)(yhat, y, class_w) ≈
        [7, 5, 2, 0] .* harm_means
    f1_macro_w = mean([7, 5, 2, 0] .* harm_means)
    @test MulticlassFScore(average=macro_avg)(yhat, y, class_w) ≈ f1_macro_w
    @test f1_macro_w ≈
        @test_logs((:warn, r"Using macro"),
                     MulticlassFScore(average=micro_avg)(yhat, y, class_w))
    f1_micro = harmonic_mean(p_micro, r_micro)
    @test MulticlassFScore(average=micro_avg)(yhat, y) ≈ f1_micro

    # fscore, β=1/3:
    harm_means = [harmonic_mean(1/2, 1/2, beta=1/3),
                     harmonic_mean(1/3, 1/3, beta=1/3),
                     harmonic_mean(1/2, 2/5, beta=1/3),
                     harmonic_mean(1/2, 1, beta=1/3)]
    f1_macro = mean(harm_means)
    @test MulticlassFScore(β=1/3, average=macro_avg)(yhat, y) ≈ f1_macro
    @test MulticlassFScore(β=1/3,
                           average=no_avg,
                           return_type=Vector)(yhat, y, class_w) ≈
        [7, 5, 2, 0] .* harm_means
    f1_macro_w = mean([7, 5, 2, 0] .* harm_means)
    @test MulticlassFScore(β=1/3,
                           average=macro_avg)(yhat, y, class_w) ≈ f1_macro_w
    @test f1_macro_w ≈
        @test_logs((:warn, r"Using macro"),
                   MulticlassFScore(β=1/3,
                                    average=micro_avg)(yhat, y, class_w))
    f1_micro = harmonic_mean(p_micro, r_micro, beta=1/3)
    @test MulticlassFScore(β=1/3, average=micro_avg)(yhat, y) ≈ f1_micro
end

@testset "docstrings coverage" begin
    @test startswith(info(BrierScore()).docstring, "`BrierScore`")
end


text:
#///////////////////////////////////////
#// File Name: utils_test.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/02/25
#// Description: Test code for src/utils.jl
#///////////////////////////////////////

import StatsBase: fit
using Distributions
using LinearAlgebra
using PyCall
using Random
using RobotOS

# RSSAC Controller Setup Tests
@testset "Controller Setup: Synthetic + Gaussian" begin
    # Scene loader & Predictor parameters
    prediction_device = "cpu";
    prediction_steps = 12;
    ado_pos_init_dict = Dict("Pedestrian/1" => [0.0, -5.0]);
    ado_vel_dict = Dict("Pedestrian/1" => MvNormal([0.0, 1.0], Diagonal([0.01, 0.3])));
    dto = 0.4;
    deterministic = false;
    prediction_rng_seed = 1;
    num_samples = 50;
    # Cost Parameters
    Cep = Matrix(1.0I, 2, 2);
    Cu = Matrix(1.0I, 2, 2);
    β_pos = 0.6;
    β_col = 0.4;
    α_col = 40.0;
    λ_col = 2.0;
    σ_risk = 0.0;
    # Control Parameters
    dtc = 0.01;
    dtr = 0.1;
    dtexec = [0.02, 0.06];
    tcalc = 0.1;
    u_norm_max = 5.0;
    u_nominal_base = [0.0, 0.0];
    u_nominal_cand = append!([u_nominal_base],
                             [round.([a*cos(deg2rad(θ)), a*sin(deg2rad(θ))], digits=5)
                              for a = [2., 4.] for θ = 0.:45.:(360. - 45.)])            # nominal control candidate value [ax, ay] [m/s^2]
    nominal_search_depth = 1;
    # Ego initial state
    ego_pos_init_vec = [-5., 0.];
    ego_pos_goal_vec = [5., 0.];
    # Other parameters
    target_speed = 1.0;
    sim_horizon = 10.0;

    rng = MersenneTwister(prediction_rng_seed);
    scene_param = SyntheticSceneParameter(rng);
    predictor_param = GaussianPredictorParameter(prediction_steps,
                                                 num_samples, deterministic, rng);
    cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
    cnt_param = ControlParameter(u_norm_max, tcalc, dtexec, dtr,
                                 u_nominal_base, u_nominal_cand, nominal_search_depth);

    scene_loader, controller, w_init, measurement_schedule, target_trajectory, ~ =
        controller_setup(scene_param,
                         predictor_param,
                         prediction_device=prediction_device,
                         cost_param=cost_param,
                         cnt_param=cnt_param,
                         dtc=dtc,
                         dto=dto,
                         ado_pos_init_dict=ado_pos_init_dict,
                         ado_vel_dict=ado_vel_dict,
                         ego_pos_init_vec=ego_pos_init_vec,
                         ego_vel_init_vec=nothing,
                         ego_pos_goal_vec=ego_pos_goal_vec,
                         target_speed=target_speed,
                         sim_horizon=sim_horizon,
                         verbose=false)

    @test scene_loader.param.rng === controller.predictor.param.rng
    @test scene_loader.ado_pos_dict == ado_pos_init_dict
    @test controller.sim_param.dtc == dtc;
    @test controller.sim_param.dto == dto;
    @test controller.sim_param.prediction_steps == prediction_steps;
    @test controller.sim_param.num_samples == num_samples;
    @test controller.sim_param.cost_param.Cep == Cep;
    @test controller.sim_param.cost_param.Cu == Cu;
    @test controller.sim_param.cost_param.β_pos == β_pos;
    @test controller.sim_param.cost_param.α_col == α_col;
    @test controller.sim_param.cost_param.β_col == β_col;
    @test controller.sim_param.cost_param.λ_col == λ_col;
    @test controller.sim_param.cost_param.σ_risk == σ_risk;
    @test controller.cnt_param.eamax == u_norm_max;
    @test controller.cnt_param.tcalc == tcalc;
    @test controller.cnt_param.dtexec == dtexec;
    @test controller.cnt_param.dtr == dtr;
    @test controller.predictor.dto == dto;
    @test controller.predictor.param.deterministic == deterministic;
    @test controller.predictor.param.prediction_steps == prediction_steps;
    @test controller.predictor.param.num_samples == num_samples;
    @test controller.predictor.ado_vel_dict == ado_vel_dict;
    @test measurement_schedule[1] == Time(controller.sim_param.dto);
    @test measurement_schedule[end] == Time(sim_horizon);
    @test minimum(target_trajectory)[1] == Time(0.0);
    @test minimum(target_trajectory)[2] == ego_pos_init_vec;
    @test maximum(target_trajectory)[1] == Time(14.8);
    @test maximum(target_trajectory)[2] == ego_pos_goal_vec;
    @test get_position(w_init.e_state) == ego_pos_init_vec;
    @test get_velocity(w_init.e_state) ==
        (ego_pos_goal_vec - ego_pos_init_vec)./
        norm(ego_pos_goal_vec - ego_pos_init_vec)*target_speed;
end

@testset "Controller Setup: Data + Gaussian" begin
    # Scene loader & Predictor parameters
    conf_file_name = "config.json";
    test_data_name = "eth_test.pkl";
    test_scene_id = 0;
    start_time_idx = 3;
    incl_robot_node = false;
    prediction_device = "cpu";
    prediction_steps = 12;
    ado_vel_dict = Dict("Any" => MvNormal([0.0, 0.0], Diagonal([0.8, 0.8])));
    deterministic = false;
    prediction_rng_seed = 1;
    num_samples = 50;
    # Cost Parameters
    Cep = Matrix(1.0I, 2, 2);
    Cu = Matrix(1.0I, 2, 2);
    β_pos = 0.6;
    β_col = 0.4;
    α_col = 40.0;
    λ_col = 2.0;
    σ_risk = 0.0;
    # Control Parameters
    dtc = 0.01;
    dtr = 0.1;
    dtexec = [0.02, 0.06];;
    tcalc = 0.1;
    u_norm_max = 5.0;
    u_nominal_base = [0.0, 0.0];
    u_nominal_cand = append!([u_nominal_base],
                             [round.([a*cos(deg2rad(θ)), a*sin(deg2rad(θ))], digits=5)
                              for a = [2., 4.] for θ = 0.:45.:(360. - 45.)])            # nominal control candidate value [ax, ay] [m/s^2]
    nominal_search_depth = 1;
    # Ego initial state
    ego_pos_init_vec = [-5., 0.];
    ego_pos_goal_vec = [5., 0.];
    # Other parameters
    target_speed = 1.0;
    sim_horizon = 10.0;

    scene_param = TrajectronSceneParameter(conf_file_name, test_data_name,
                                           test_scene_id, start_time_idx,
                                           incl_robot_node);
    rng = MersenneTwister(prediction_rng_seed);
    predictor_param = GaussianPredictorParameter(prediction_steps,
                                                 num_samples, deterministic, rng);
    cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
    cnt_param = ControlParameter(u_norm_max, tcalc, dtexec, dtr,
                                 u_nominal_base, u_nominal_cand, nominal_search_depth);

    scene_loader, controller, w_init, measurement_schedule, target_trajectory, ~ =
    controller_setup(scene_param,
                     predictor_param,
                     prediction_device=prediction_device,
                     cost_param=cost_param,
                     cnt_param=cnt_param,
                     dtc=dtc,
                     ado_vel_dict=ado_vel_dict,
                     ego_pos_init_vec=ego_pos_init_vec,
                     ego_vel_init_vec=nothing,
                     ego_pos_goal_vec=ego_pos_goal_vec,
                     target_speed=target_speed,
                     sim_horizon=sim_horizon,
                     verbose=false)

    @test scene_loader.param.conf_file_name == conf_file_name;
    @test scene_loader.param.test_data_name == test_data_name;
    @test scene_loader.param.test_scene_id == test_scene_id;
    @test scene_loader.param.start_time_idx == start_time_idx;
    @test scene_loader.param.incl_robot_node == false;
    @test scene_loader.curr_time_idx == start_time_idx + 1;
    @test controller.sim_param.dtc == dtc;
    @test controller.sim_param.dto == scene_loader.dto;
    @test controller.sim_param.prediction_steps == prediction_steps;
    @test controller.sim_param.num_samples == num_samples;
    @test controller.sim_param.cost_param.Cep == Cep;
    @test controller.sim_param.cost_param.Cu == Cu;
    @test controller.sim_param.cost_param.β_pos == β_pos;
    @test controller.sim_param.cost_param.α_col == α_col;
    @test controller.sim_param.cost_param.β_col == β_col;
    @test controller.sim_param.cost_param.λ_col == λ_col;
    @test controller.sim_param.cost_param.σ_risk == σ_risk;
    @test controller.cnt_param.eamax == u_norm_max;
    @test controller.cnt_param.tcalc == tcalc;
    @test controller.cnt_param.dtexec == dtexec;
    @test controller.cnt_param.dtr == dtr;
    @test controller.predictor.dto == scene_loader.dto;
    @test controller.predictor.param.deterministic == deterministic;
    @test controller.predictor.param.prediction_steps == prediction_steps;
    @test controller.predictor.param.num_samples == num_samples;
    @test controller.predictor.ado_vel_dict == ado_vel_dict;
    @test measurement_schedule[1] == Time(controller.sim_param.dto);
    @test measurement_schedule[end] == Time(sim_horizon);
    @test minimum(target_trajectory)[1] == Time(0.0);
    @test minimum(target_trajectory)[2] == ego_pos_init_vec;
    @test maximum(target_trajectory)[1] == Time(14.8);
    @test maximum(target_trajectory)[2] == ego_pos_goal_vec;
    @test get_position(w_init.e_state) == ego_pos_init_vec;
    @test get_velocity(w_init.e_state) ==
        (ego_pos_goal_vec - ego_pos_init_vec)./
        norm(ego_pos_goal_vec - ego_pos_init_vec)*target_speed;
end

@testset "Controller Setup: Data + Trajectron" begin
    # Scene loader & Predictor parameters
    conf_file_name = "config.json";
    test_data_name = "eth_test.pkl";
    test_scene_id = 0;
    start_time_idx = 2;
    incl_robot_node = false;
    prediction_device = "cpu";
    deterministic = false;
    prediction_rng_seed = 1;
    prediction_steps = 12;
    num_samples = 50;
    use_robot_future = false;
    # Cost Parameters
    Cep = Matrix(1.0I, 2, 2);
    Cu = Matrix(1.0I, 2, 2);
    β_pos = 0.6;
    β_col = 0.4;
    α_col = 40.0;
    λ_col = 2.0;
    σ_risk = 0.0;
    # Control Parameters
    dtc = 0.01;
    dtr = 0.1;
    dtexec = [0.02, 0.06];;
    tcalc = 0.1;
    u_norm_max = 5.0;
    u_nominal_base = [0.0, 0.0];
    u_nominal_cand = append!([u_nominal_base],
                             [round.([a*cos(deg2rad(θ)), a*sin(deg2rad(θ))], digits=5)
                              for a = [2., 4.] for θ = 0.:45.:(360. - 45.)])            # nominal control candidate value [ax, ay] [m/s^2]
    nominal_search_depth = 1;
    # Ego initial state
    ego_pos_init_vec = [5.0, 0.5];
    ego_pos_goal_vec = [5.0, 7.5];
    # Other parameters
    target_speed = 0.7;
    sim_horizon = 10.0;

    scene_param = TrajectronSceneParameter(conf_file_name, test_data_name,
                                           test_scene_id, start_time_idx,
                                           incl_robot_node);
    predictor_param = TrajectronPredictorParameter(prediction_steps,
                                                   num_samples,
                                                   use_robot_future,
                                                   deterministic,
                                                   prediction_rng_seed);
    cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
    cnt_param = ControlParameter(u_norm_max, tcalc, dtexec, dtr,
                                 u_nominal_base, u_nominal_cand, nominal_search_depth);
    scene_loader, controller, w_init, measurement_schedule, target_trajectory, ~ =
    controller_setup(scene_param,
                     predictor_param,
                     prediction_device=prediction_device,
                     cost_param=cost_param,
                     cnt_param=cnt_param,
                     dtc=dtc,
                     ego_pos_init_vec=ego_pos_init_vec,
                     ego_vel_init_vec=nothing,
                     ego_pos_goal_vec=ego_pos_goal_vec,
                     target_speed=target_speed,
                     sim_horizon=sim_horizon,
                     verbose=false)

    @test scene_loader.param.conf_file_name == conf_file_name;
    @test scene_loader.param.test_data_name == test_data_name;
    @test scene_loader.param.test_scene_id == test_scene_id;
    @test scene_loader.param.start_time_idx == start_time_idx;
    @test scene_loader.param.incl_robot_node == false;
    @test scene_loader.curr_time_idx == start_time_idx + 1;
    @test controller.sim_param.dtc == dtc;
    @test controller.sim_param.dto == scene_loader.dto;
    @test controller.sim_param.prediction_steps == prediction_steps;
    @test controller.sim_param.num_samples == num_samples;
    @test controller.sim_param.cost_param.Cep == Cep;
    @test controller.sim_param.cost_param.Cu == Cu;
    @test controller.sim_param.cost_param.β_pos == β_pos;
    @test controller.sim_param.cost_param.α_col == α_col;
    @test controller.sim_param.cost_param.β_col == β_col;
    @test controller.sim_param.cost_param.λ_col == λ_col;
    @test controller.sim_param.cost_param.σ_risk == σ_risk;
    @test controller.cnt_param.eamax == u_norm_max;
    @test controller.cnt_param.tcalc == tcalc;
    @test controller.cnt_param.dtexec == dtexec;
    @test controller.cnt_param.dtr == dtr;
    @test controller.predictor.param.deterministic == deterministic;
    @test controller.predictor.param.prediction_steps == prediction_steps;
    @test controller.predictor.param.num_samples == num_samples;
    @test controller.predictor.param.use_robot_future == false;
    @test controller.predictor.param.deterministic == false;
    @test controller.predictor.param.rng_seed_py == prediction_rng_seed;
    @test istaskdone(controller.prediction_task);
    @test !isnothing(controller.prediction_dict_tmp);
    @test w_init.t == Time(0.0);
    @test keys(w_init.ap_dict) == keys(controller.prediction_dict_tmp);
    @test w_init.t_last_m == Time(0.0);
    @test measurement_schedule[1] == Time(controller.sim_param.dto);
    @test measurement_schedule[end] == Time(sim_horizon);
    @test minimum(target_trajectory)[1] == Time(0.0);
    @test minimum(target_trajectory)[2] == ego_pos_init_vec;
    @test maximum(target_trajectory)[1] == Time(14.8);
    @test maximum(target_trajectory)[2] == ego_pos_goal_vec;
    @test get_position(w_init.e_state) == ego_pos_init_vec;
    @test get_velocity(w_init.e_state) ==
        (ego_pos_goal_vec - ego_pos_init_vec)./
        norm(ego_pos_goal_vec - ego_pos_init_vec)*target_speed;
end

@testset "Controller Setup: Data (Ado Replaced) + Trajectron" begin
    # Scene loader & Predictor parameters
    conf_file_name = "config.json";
    test_data_name = "eth_test.pkl";
    test_scene_id = 0;
    start_time_idx = 2;
    incl_robot_node = false;
    prediction_device = "cpu";
    deterministic = false;
    prediction_rng_seed = 1;
    prediction_steps = 12;
    num_samples = 50;
    use_robot_future = false;
    # Cost Parameters
    Cep = Matrix(1.0I, 2, 2);
    Cu = Matrix(1.0I, 2, 2);
    β_pos = 0.6;
    β_col = 0.4;
    α_col = 40.0;
    λ_col = 2.0;
    σ_risk = 0.0;
    # Control Parameters
    dtc = 0.01;
    dtr = 0.1;
    dtexec = [0.02, 0.06];;
    tcalc = 0.1;
    u_norm_max = 5.0;
    u_nominal_base = [0.0, 0.0];
    u_nominal_cand = append!([u_nominal_base],
                             [round.([a*cos(deg2rad(θ)), a*sin(deg2rad(θ))], digits=5)
                              for a = [2., 4.] for θ = 0.:45.:(360. - 45.)])            # nominal control candidate value [ax, ay] [m/s^2]
    nominal_search_depth = 1;
    # Other parameters
    sim_horizon = 10.0;

    scene_param = TrajectronSceneParameter(conf_file_name, test_data_name,
                                           test_scene_id, start_time_idx,
                                           incl_robot_node);
    predictor_param = TrajectronPredictorParameter(prediction_steps,
                                                   num_samples,
                                                   use_robot_future,
                                                   deterministic,
                                                   prediction_rng_seed);
    cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
    cnt_param = ControlParameter(u_norm_max, tcalc, dtexec, dtr,
                                 u_nominal_base, u_nominal_cand, nominal_search_depth);
    scene_loader, controller, w_init, measurement_schedule, target_trajectory, target_speed =
    controller_setup(scene_param,
                     predictor_param,
                     prediction_device=prediction_device,
                     cost_param=cost_param,
                     cnt_param=cnt_param,
                     dtc=dtc,
                     sim_horizon=sim_horizon,
                     ado_id_to_replace="PEDESTRIAN/2",
                     verbose=false)

    @test scene_loader.param.conf_file_name == conf_file_name;
    @test scene_loader.param.test_data_name == test_data_name;
    @test scene_loader.param.test_scene_id == test_scene_id;
    @test scene_loader.param.start_time_idx == start_time_idx;
    @test scene_loader.param.incl_robot_node == false;
    @test scene_loader.curr_time_idx == start_time_idx + 1;
    @test controller.sim_param.dtc == dtc;
    @test controller.sim_param.dto == scene_loader.dto;
    @test controller.sim_param.prediction_steps == prediction_steps;
    @test controller.sim_param.num_samples == num_samples;
    @test controller.sim_param.cost_param.Cep == Cep;
    @test controller.sim_param.cost_param.Cu == Cu;
    @test controller.sim_param.cost_param.β_pos == β_pos;
    @test controller.sim_param.cost_param.α_col == α_col;
    @test controller.sim_param.cost_param.β_col == β_col;
    @test controller.sim_param.cost_param.λ_col == λ_col;
    @test controller.sim_param.cost_param.σ_risk == σ_risk;
    @test controller.cnt_param.eamax == u_norm_max;
    @test controller.cnt_param.tcalc == tcalc;
    @test controller.cnt_param.dtexec == dtexec;
    @test controller.cnt_param.dtr == dtr;
    @test controller.predictor.param.deterministic == deterministic;
    @test controller.predictor.param.prediction_steps == prediction_steps;
    @test controller.predictor.param.num_samples == num_samples;
    @test controller.predictor.param.use_robot_future == false;
    @test controller.predictor.param.deterministic == false;
    @test controller.predictor.param.rng_seed_py == prediction_rng_seed;
    @test istaskdone(controller.prediction_task);
    @test !isnothing(controller.prediction_dict_tmp);
    @test w_init.t == Time(0.0);
    @test !in("PEDESTRIAN/2", keys(w_init.ap_dict))
    @test keys(w_init.ap_dict) == keys(controller.prediction_dict_tmp);
    @test w_init.t_last_m == Time(0.0);
    @test measurement_schedule[1] == Time(controller.sim_param.dto);
    @test measurement_schedule[end] == Time(sim_horizon);
    @test minimum(target_trajectory)[1] == Time(0.0);
    @test maximum(target_trajectory)[1] == Time(14.8);
    @test minimum(target_trajectory)[2] == get_position(w_init.e_state);
    @test all(get_velocity(w_init.e_state) .≈
              (collect(values(target_trajectory))[2] - minimum(target_trajectory)[2])./scene_loader.dto)
    @test target_speed ≈ norm((maximum(target_trajectory)[2] - minimum(target_trajectory)[2])./
                               to_sec(collect(keys(target_trajectory))[end-1] - minimum(target_trajectory)[1]))
end

@testset "Controller Setup: Data + Oracle" begin
    # Scene loader & Predictor parameters
    conf_file_name = "config.json";
    test_data_name = "eth_test.pkl";
    test_scene_id = 0;
    start_time_idx = 3;
    incl_robot_node = false;
    prediction_device = "cpu";
    deterministic = true;
    prediction_steps = 12;
    num_samples = 1;
    # Cost Parameters
    Cep = Matrix(1.0I, 2, 2);
    Cu = Matrix(1.0I, 2, 2);
    β_pos = 0.6;
    β_col = 0.4;
    α_col = 40.0;
    λ_col = 2.0;
    σ_risk = 0.0;
    # Control Parameters
    dtc = 0.01;
    dtr = 0.1;
    dtexec = [0.02, 0.06];;
    tcalc = 0.1;
    u_norm_max = 5.0;
    u_nominal_base = [0.0, 0.0];
    u_nominal_cand = append!([u_nominal_base],
                             [round.([a*cos(deg2rad(θ)), a*sin(deg2rad(θ))], digits=5)
                              for a = [2., 4.] for θ = 0.:45.:(360. - 45.)])            # nominal control candidate value [ax, ay] [m/s^2]
    nominal_search_depth = 1;
    # Ego initial state
    ego_pos_init_vec = [5.0, 0.5];
    ego_pos_goal_vec = [5.0, 7.5];
    # Other parameters
    target_speed = 0.7;
    sim_horizon = 10.0;

    scene_param = TrajectronSceneParameter(conf_file_name, test_data_name,
                                           test_scene_id, start_time_idx,
                                           incl_robot_node);
    predictor_param = OraclePredictorParameter(prediction_steps);
    cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
    cnt_param = ControlParameter(u_norm_max, tcalc, dtexec, dtr,
                                 u_nominal_base, u_nominal_cand, nominal_search_depth);

    scene_loader, controller, w_init, measurement_schedule, target_trajectory, ~ =
    controller_setup(scene_param,
                     predictor_param,
                     prediction_device=prediction_device,
                     cost_param=cost_param,
                     cnt_param=cnt_param,
                     dtc=dtc,
                     ego_pos_init_vec=ego_pos_init_vec,
                     ego_vel_init_vec=nothing,
                     ego_pos_goal_vec=ego_pos_goal_vec,
                     target_speed=target_speed,
                     sim_horizon=sim_horizon,
                     verbose=false)

    @test scene_loader.param.conf_file_name == conf_file_name;
    @test scene_loader.param.test_data_name == test_data_name;
    @test scene_loader.param.test_scene_id == test_scene_id;
    @test scene_loader.param.start_time_idx == start_time_idx;
    @test scene_loader.param.incl_robot_node == false;
    @test scene_loader.curr_time_idx == start_time_idx + 1;
    @test controller.sim_param.dtc == dtc;
    @test controller.sim_param.dto == scene_loader.dto;
    @test controller.sim_param.prediction_steps == prediction_steps;
    @test controller.sim_param.num_samples == num_samples;
    @test controller.sim_param.cost_param.Cep == Cep;
    @test controller.sim_param.cost_param.Cu == Cu;
    @test controller.sim_param.cost_param.β_pos == β_pos;
    @test controller.sim_param.cost_param.α_col == α_col;
    @test controller.sim_param.cost_param.β_col == β_col;
    @test controller.sim_param.cost_param.λ_col == λ_col;
    @test controller.sim_param.cost_param.σ_risk == σ_risk;
    @test controller.cnt_param.eamax == u_norm_max;
    @test controller.cnt_param.tcalc == tcalc;
    @test controller.cnt_param.dtexec == dtexec;
    @test controller.cnt_param.dtr == dtr;
    @test controller.predictor.param.prediction_steps == prediction_steps;
    @test controller.predictor.eval_scene == scene_loader.eval_scene;
    @test controller.predictor.state_def == scene_loader.state_def;
    @test controller.predictor.curr_time_idx == scene_loader.curr_time_idx
    @test controller.predictor.max_time_idx == scene_loader.max_time_idx
    @test istaskdone(controller.prediction_task);
    @test !isnothing(controller.prediction_dict_tmp);
    @test w_init.t == Time(0.0);
    @test keys(w_init.ap_dict) == keys(controller.prediction_dict_tmp);
    @test w_init.t_last_m == Time(0.0);
    @test measurement_schedule[1] == Time(controller.sim_param.dto);
    @test measurement_schedule[end] == Time(sim_horizon);
    @test minimum(target_trajectory)[1] == Time(0.0);
    @test minimum(target_trajectory)[2] == ego_pos_init_vec;
    @test maximum(target_trajectory)[1] == Time(14.8);
    @test maximum(target_trajectory)[2] == ego_pos_goal_vec;
    @test get_position(w_init.e_state) == ego_pos_init_vec;
    @test get_velocity(w_init.e_state) ==
        (ego_pos_goal_vec - ego_pos_init_vec)./
        norm(ego_pos_goal_vec - ego_pos_init_vec)*target_speed;
end

# BICController Setup Tests
@testset "Controller Setup: Synthetic + BIC" begin
    # Scene loader & Predictor parameters
    prediction_steps = 12;
    ado_pos_init_dict = Dict("Pedestrian/1" => [0.0, -5.0]);
    ado_vel_dict = Dict("Pedestrian/1" => MvNormal([0.0, 1.0], Diagonal([0.01, 0.3])));
    dto = 0.4;
    sim_rng_seed = 1;
    num_samples = 50;
    # Cost Parameters
    Cep = Matrix(1.0I, 2, 2);
    Cu = Matrix(1.0I, 2, 2);
    β_pos = 0.6;
    β_col = 0.4;
    α_col = 40.0;
    λ_col = 2.0;
    σ_risk = 0.0;
    # Control Parameters
    dtc = 0.01;
    dtr = 0.1;
    tcalc = 0.1;
    u_norm_max = 5.0;
    min_dist = 0.5;
    # Ego initial state
    ego_pos_init_vec = [-5., 0.];
    ego_pos_goal_vec = [5., 0.];
    # Other parameters
    target_speed = 1.0;
    sim_horizon = 10.0;

    rng = MersenneTwister(sim_rng_seed);
    scene_param = SyntheticSceneParameter(rng);
    cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
    cnt_param = BICControlParameter(u_norm_max, tcalc, dtr, min_dist);

    scene_loader, controller, w_init, measurement_schedule, target_trajectory, ~ =
    controller_setup(scene_param,
                     cost_param=cost_param,
                     cnt_param=cnt_param,
                     dtc=dtc,
                     dto=dto,
                     prediction_steps=prediction_steps,
                     num_samples=num_samples,
                     sim_rng=rng,
                     ado_pos_init_dict=ado_pos_init_dict,
                     ado_vel_dict=ado_vel_dict,
                     ego_pos_init_vec=ego_pos_init_vec,
                     ego_vel_init_vec=nothing,
                     ego_pos_goal_vec=ego_pos_goal_vec,
                     target_speed=target_speed,
                     sim_horizon=sim_horizon,
                     verbose=false)

    @test scene_loader.ado_pos_dict == ado_pos_init_dict
    @test controller.sim_param.dtc == dtc;
    @test controller.sim_param.dto == dto;
    @test controller.sim_param.prediction_steps == prediction_steps;
    @test controller.sim_param.num_samples == num_samples;
    @test controller.sim_param.cost_param.Cep == Cep;
    @test controller.sim_param.cost_param.Cu == Cu;
    @test controller.sim_param.cost_param.β_pos == β_pos;
    @test controller.sim_param.cost_param.α_col == α_col;
    @test controller.sim_param.cost_param.β_col == β_col;
    @test controller.sim_param.cost_param.λ_col == λ_col;
    @test controller.sim_param.cost_param.σ_risk == σ_risk;
    @test controller.cnt_param.eamax == u_norm_max;
    @test controller.cnt_param.tcalc == tcalc;
    @test controller.cnt_param.dtr == dtr;
    @test controller.cnt_param.min_dist == min_dist;
    @test measurement_schedule[1] == Time(controller.sim_param.dto);
    @test measurement_schedule[end] == Time(sim_horizon);
    @test minimum(target_trajectory)[1] == Time(0.0);
    @test minimum(target_trajectory)[2] == ego_pos_init_vec;
    @test maximum(target_trajectory)[1] == Time(14.8);
    @test maximum(target_trajectory)[2] == ego_pos_goal_vec;
    @test get_position(w_init.e_state) == ego_pos_init_vec;
    @test get_velocity(w_init.e_state) ==
        (ego_pos_goal_vec - ego_pos_init_vec)./
        norm(ego_pos_goal_vec - ego_pos_init_vec)*target_speed;
end

@testset "Controller Setup: Data + BIC" begin
    # Scene loader & Predictor parameters
    conf_file_name = "config.json";
    test_data_name = "eth_test.pkl";
    test_scene_id = 0;
    start_time_idx = 3;
    incl_robot_node = false;
    prediction_steps = 12;
    # Cost Parameters
    Cep = Matrix(1.0I, 2, 2);
    Cu = Matrix(1.0I, 2, 2);
    β_pos = 0.6;
    β_col = 0.4;
    α_col = 40.0;
    λ_col = 2.0;
    σ_risk = 0.0;
    # Control Parameters
    dtc = 0.01;
    dtr = 0.1;
    tcalc = 0.1;
    u_norm_max = 5.0;
    min_dist = 0.5;
    # Ego initial state
    ego_pos_init_vec = [-5., 0.];
    ego_pos_goal_vec = [5., 0.];
    # Other parameters
    target_speed = 1.0;
    sim_horizon = 10.0;

    scene_param = TrajectronSceneParameter(conf_file_name, test_data_name,
                                           test_scene_id, start_time_idx,
                                           incl_robot_node);
    cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
    cnt_param = BICControlParameter(u_norm_max, tcalc, dtr, min_dist);

    scene_loader, controller, w_init, measurement_schedule, target_trajectory, ~ =
    controller_setup(scene_param,
                     cost_param=cost_param,
                     cnt_param=cnt_param,
                     dtc=dtc,
                     prediction_steps=prediction_steps,
                     ego_pos_init_vec=ego_pos_init_vec,
                     ego_pos_goal_vec=ego_pos_goal_vec,
                     target_speed=target_speed,
                     sim_horizon=sim_horizon,
                     verbose=false);

    @test scene_loader.param.conf_file_name == conf_file_name;
    @test scene_loader.param.test_data_name == test_data_name;
    @test scene_loader.param.test_scene_id == test_scene_id;
    @test scene_loader.param.start_time_idx == start_time_idx;
    @test scene_loader.param.incl_robot_node == false;
    @test scene_loader.curr_time_idx == start_time_idx + 1;
    @test controller.sim_param.dtc == dtc;
    @test controller.sim_param.dto == scene_loader.dto;
    @test controller.sim_param.prediction_steps == prediction_steps;
    @test controller.sim_param.num_samples == 0;
    @test controller.sim_param.cost_param.Cep == Cep;
    @test controller.sim_param.cost_param.Cu == Cu;
    @test controller.sim_param.cost_param.β_pos == β_pos;
    @test controller.sim_param.cost_param.α_col == α_col;
    @test controller.sim_param.cost_param.β_col == β_col;
    @test controller.sim_param.cost_param.λ_col == λ_col;
    @test controller.sim_param.cost_param.σ_risk == σ_risk;
    @test controller.cnt_param.eamax == u_norm_max;
    @test controller.cnt_param.tcalc == tcalc;
    @test controller.cnt_param.dtr == dtr;
    @test controller.cnt_param.min_dist == min_dist;
    @test measurement_schedule[1] == Time(controller.sim_param.dto);
    @test measurement_schedule[end] == Time(sim_horizon);
    @test minimum(target_trajectory)[1] == Time(0.0);
    @test minimum(target_trajectory)[2] == ego_pos_init_vec;
    @test maximum(target_trajectory)[1] == Time(14.8);
    @test maximum(target_trajectory)[2] == ego_pos_goal_vec;
    @test get_position(w_init.e_state) == ego_pos_init_vec;
    @test get_velocity(w_init.e_state) ==
        (ego_pos_goal_vec - ego_pos_init_vec)./
        norm(ego_pos_goal_vec - ego_pos_init_vec)*target_speed;
end

@testset "Controller Setup: Data + CrowdNav" begin
    # Scene loader & Predictor parameters
    conf_file_name = "config.json";
    test_data_name = "eth_test.pkl";
    test_scene_id = 0;
    start_time_idx = 3;
    incl_robot_node = false;
    prediction_steps = 12;
    # Cost Parameters
    Cep = Matrix(1.0I, 2, 2);
    Cu = Matrix(1.0I, 2, 2);
    β_pos = 0.6;
    β_col = 0.4;
    α_col = 40.0;
    λ_col = 2.0;
    σ_risk = 0.0;
    # Control Parameters
    model_dir = normpath(joinpath(@__DIR__,
                                  "../CrowdNav/crowd_nav/data/output_om_sarl_radius_0.4"));
    env_config = "env.config";
    policy_config = "policy.config";
    policy_name = "sarl";
    tcalc = 0.01;
    dtr = 0.4;
    # Ego initial state
    ego_pos_init_vec = [-5., 0.];
    ego_pos_goal_vec = [5., 0.];
    # Other parameters
    dtc = 0.01;
    target_speed = 1.0;
    sim_horizon = 10.0;

    scene_param = TrajectronSceneParameter(conf_file_name, test_data_name,
                                           test_scene_id, start_time_idx,
                                           incl_robot_node);
    cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
    cnt_param = CrowdNavControlParameter(model_dir, env_config, policy_config,
                                         policy_name, ego_pos_goal_vec,
                                         tcalc, dtr, target_speed);

    scene_loader, controller, w_init, ado_inputs, measurement_schedule,
    target_trajectory, target_speed =
    controller_setup(scene_param, cnt_param,
                     cost_param=cost_param,
                     dtc=dtc,
                     prediction_steps=prediction_steps,
                     ego_pos_init_vec=ego_pos_init_vec,
                     ego_pos_goal_vec=ego_pos_goal_vec,
                     target_speed=target_speed,
                     sim_horizon=sim_horizon,
                     verbose=false);
    @test scene_loader.param.conf_file_name == conf_file_name;
    @test scene_loader.param.test_data_name == test_data_name;
    @test scene_loader.param.test_scene_id == test_scene_id;
    @test scene_loader.param.start_time_idx == start_time_idx;
    @test scene_loader.param.incl_robot_node == false;
    @test scene_loader.curr_time_idx == start_time_idx + 1;
    @test controller.sim_param.dtc == dtc;
    @test controller.sim_param.dto == scene_loader.dto;
    @test controller.sim_param.prediction_steps == prediction_steps;
    @test controller.sim_param.num_samples == 0;
    @test controller.sim_param.cost_param.Cep == Cep;
    @test controller.sim_param.cost_param.Cu == Cu;
    @test controller.sim_param.cost_param.β_pos == β_pos;
    @test controller.sim_param.cost_param.α_col == α_col;
    @test controller.sim_param.cost_param.β_col == β_col;
    @test controller.sim_param.cost_param.λ_col == λ_col;
    @test controller.sim_param.cost_param.σ_risk == σ_risk;
    @test controller.cnt_param.model_dir == model_dir;
    @test controller.cnt_param.env_config == env_config;
    @test controller.cnt_param.policy_config == policy_config;
    @test controller.cnt_param.policy_name == policy_name;
    @test controller.cnt_param.goal_pos == ego_pos_goal_vec;
    @test controller.cnt_param.tcalc == tcalc;
    @test controller.cnt_param.dtr == dtr;
    @test measurement_schedule[1] == Time(controller.sim_param.dto);
    @test measurement_schedule[end] == Time(sim_horizon);
    @test minimum(target_trajectory)[1] == Time(0.0);
    @test minimum(target_trajectory)[2] == ego_pos_init_vec;
    @test maximum(target_trajectory)[1] == Time(14.8);
    @test maximum(target_trajectory)[2] == ego_pos_goal_vec;
    @test get_position(w_init.e_state) == ego_pos_init_vec;
    @test get_velocity(w_init.e_state) ==
        (ego_pos_goal_vec - ego_pos_init_vec)./
        norm(ego_pos_goal_vec - ego_pos_init_vec)*target_speed;
end


text:
using MatpowerCases, Test

println("- extract_case check")
for c in casenames()
    mpc = loadcase(c; describe=false)
    @assert isa(extract_case(mpc), MatpowerCases.Case)
    @assert isa(extract_case(c), MatpowerCases.Case)
end

c = joinpath(@__DIR__, "..", "data", "case118.mat")
loadcase(c; describe=false)
@assert isa(extract_case(c), MatpowerCases.Case)


text:
function clear_all!(window)
    wh = widths(window)
    glViewport(0,0, wh...)
    fb = framebuffer(window)
    glBindFramebuffer(GL_FRAMEBUFFER, fb.id)
    glDrawBuffers(2, [GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1])
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glBindFramebuffer(GL_FRAMEBUFFER, 0)
    glClear(GL_COLOR_BUFFER_BIT)
end

"""
Renders a single frame of a `window`
"""
function render_frame(window)
    fb = framebuffer(window)
    wh = widths(window)
    resize!(fb, wh)
    prepare(fb)
    glViewport(0,0, wh...)
    glClearColor(1,1,1,1)
    glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT)
    render(window)
    #Read all the selection queries
    push_selectionqueries!(window)
    display(fb, window)
end

function renderloop(window::Screen)
    while isopen(window)
        render_frame(window)
        swapbuffers(window)
        GLFW.PollEvents()
        yield()
    end
    destroy!(window)
end

function prepare(fb::GLFramebuffer)
    glDisable(GL_SCISSOR_TEST)
    glBindFramebuffer(GL_FRAMEBUFFER, fb.id)
    glDrawBuffers(2, [GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1])
end

function display(fb::GLFramebuffer, window)
    glDisable(GL_SCISSOR_TEST)
    glBindFramebuffer(GL_FRAMEBUFFER, 0)
    glViewport(0,0, widths(window)...)
    render(fb.postprocess)
end

function GLAbstraction.render(x::Screen, parent::Screen=x, context=x.area.value)
    if isopen(x) && !ishidden(x)
        sa    = value(x.area)
        sa    = SimpleRectangle(context.x+sa.x, context.y+sa.y, sa.w, sa.h) # bring back to absolute values
        pa    = context
        sa_pa = intersect(pa, sa) # intersection with parent
        if (
            sa_pa != SimpleRectangle(0,0,0,0) && # if it is in the parent area
            (sa_pa.w > 0 && sa_pa.h > 0)
        ) # if it is in the parent area
            glEnable(GL_SCISSOR_TEST)
            glScissor(sa_pa)
            glViewport(sa)
            colorbits = GL_DEPTH_BUFFER_BIT
            if alpha(x.color) > 0
                glClearColor(red(x.color), green(x.color), blue(x.color), alpha(x.color))
                colorbits = colorbits | GL_COLOR_BUFFER_BIT
            end
            glClear(colorbits)

            render(x.renderlist)
            for window in x.children
                render(window, x, sa)
            end
        end
    end
end


text:
using ITensors

examples_dir = joinpath(ITensors.examples_dir(), "src")
include(joinpath(examples_dir, "ctmrg_isotropic.jl"))
include(joinpath(examples_dir, "2d_classical_ising.jl"))

function run(; maxdim::Int,
               nsweeps::Int = 800,
               outputlevel::Int = 0,
               β::Float64 = 1.001 * βc)
  # Make Ising model MPO
  s = Index(2, "Site")
  sₕ = addtags(s, "horiz")
  sᵥ = addtags(s, "vert")

  T = ising_mpo(sₕ, sᵥ, β)

  χ0 = 1
  l = Index(χ0, "Link")
  lₕ = addtags(l, "horiz")
  lᵥ = addtags(l, "vert")

  # Initial CTM
  Cₗᵤ = ITensor(lᵥ, lₕ)
  Cₗᵤ[1, 1] = 1.0

  # Initial HRTM
  Aₗ = ITensor(lᵥ, lᵥ', sₕ)
  Aₗ[lᵥ => 1, lᵥ' => 1, sₕ => 1] = 1.0

  Cₗᵤ, Aₗ = ctmrg(T, Cₗᵤ, Aₗ;
                  χmax = maxdim,
                  cutoff = 0.0,
                  nsteps = nsweeps)

  lᵥ = commonind(Cₗᵤ, Aₗ)
  lₕ = noncommoninds(Cₗᵤ, Aₗ)[1]

  Aᵤ = replaceinds(Aₗ, lᵥ => lₕ, lᵥ' => lₕ', sₕ => sᵥ)

  ACₗ = Aₗ * Cₗᵤ * dag(Cₗᵤ')

  ACTₗ = prime(ACₗ * dag(Aᵤ') * T * Aᵤ, -1)
  κ = (ACTₗ * dag(ACₗ))[]

  Tsz = ising_mpo(sₕ, sᵥ, β; sz = true)
  ACTszₗ = prime(ACₗ * dag(Aᵤ') * Tsz * Aᵤ, -1)
  m = (ACTszₗ * dag(ACₗ))[] / κ

  if outputlevel > 0
    @show nsweeps
    @show κ, exp(-β * ising_free_energy(β))
    @show m, ising_magnetization(β)
  end

  return ITensors.maxdim(Cₗᵤ)
end


text:
using DashBootstrapComponents, DashHtmlComponents

textareas = html_div([
    dbc_textarea(className = "mb-3", placeholder = "A Textarea"),
    dbc_textarea(
        valid = true,
        size = "sm",
        className = "mb-3",
        placeholder = "A small, valid Textarea",
    ),
    dbc_textarea(invalid = true, size = "lg", placeholder = "A large, invalid Textarea"),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
# This file is part of Kpax3. License is MIT.

function selection(fitness::Array{Float64})
  (a1, a2, b1, b2) = StatsBase.sample(1:length(fitness), 4, replace=false)

  i1 = fitness[a1] > fitness[a2] ? a1 : a2
  i2 = fitness[b1] > fitness[b2] ? b1 : b2

  (i1, i2)
end


text:
# `CompressedFinder{T <: Integer}` is an optimized variant of UnionFinder{T}
# which does not support the addition of more edges. `CompressedFinder` also
# tracks the number of groups in the graph and garuantees that all group IDs
# will be between 1 and the total group number, inclusive.
mutable struct CompressedFinder{T <: Integer}
    ids :: Vector{T}
    groups :: T

    # `CompressedFinder(uf)` creates a `CompressedFinder` instance from the
    # groups within `uf`.
    function CompressedFinder{T}(uf :: UnionFinder{T}) where T
        groups = zero(T)
        ids = zeros(T, length(uf.parents))

        for i in one(T):convert(T, length(uf.parents))
            root = find!(uf, i)
            if ids[root] == 0
                groups += 1
                ids[root] = groups
            end
            ids[i] = ids[root]
        end

        return new(ids, convert(T, groups))
    end
end

CompressedFinder(uf :: UnionFinder) = CompressedFinder{eltype(uf.sizes)}(uf)

# `find(cf, node)` returns the group ID of `node`. `node` must be a valid
# index into `uf`.
function find(cf :: CompressedFinder{T}, node :: T) where T <: Integer
    if node <= 0 || node > length(cf.ids)
        throw(BoundsError())
    end

    return cf.ids[node]
end

# `length(cf)` returns the number of nodes in `cf`.
function Base.length(cf :: CompressedFinder)
    return length(cf.ids)
end

# `groups(cf)` returns the number of groups in `cf`.
function groups(cf :: CompressedFinder)
    return cf.groups
end


text:
using Makie
using MultipleViewGeometry, Test, Random
using MultipleViewGeometry.ModuleCostFunction
using MultipleViewGeometry.ModuleTypes
using MultipleViewGeometry.ModuleConstraints
using MultipleViewGeometry.ModuleConstruct
using MultipleViewGeometry.ModuleDraw
using MultipleViewGeometry.ModuleMove
using LinearAlgebra
using StaticArrays
using GeometryTypes
using Test

𝒳 = [Point3D(x,y,rand(50:100)) for x = -100:5:100 for y = -100:5:100]
𝒳 = 𝒳[1:50:end]
X = reshape(reinterpret(Float64,𝒳),(3,length(𝒳)))

# Specify the coordinate systems of the world, the camera frame and the picture
# plane.
world_basis = (Vec(1.0, 0.0, 0.0), Vec(0.0, 1.0, 0.0), Vec(0.0, 0.0, 1.0))
camera_basis = (Point(0.0, 0.0, 0.0), Vec(-1.0, 0.0, 0.0), Vec(0.0, -1.0, 0.0), Vec(0.0, 0.0, 1.0))
picture_basis = (Point(0.0, 0.0), Vec(-1.0, 0.0), Vec(0.0, -1.0))

# The focal length for both cameras is one.
f = 1
image_width = 640 / 10
image_height = 480 / 10

camera₁ = Pinhole(image_width, image_height, f, camera_basis..., picture_basis...)
camera₂ = Pinhole(image_width, image_height, f, camera_basis..., picture_basis...)

# Rotate and translate camera one.
𝐑₁ = Matrix{Float64}(I,3,3)
𝐭₁ = [-50.0, -2.0, 0.0]
relocate!(camera₁, 𝐑₁, 𝐭₁)

# Rotate and translate camera two.
𝐑₂ = Matrix{Float64}(I,3,3)
𝐭₂ = [50.0, 2.0, 0.0]
relocate!(camera₂, 𝐑₂, 𝐭₂)

scale = 20.0f0
x = Vec3f0(0); baselen = 0.2f0 * scale ; dirlen = 1f0 * scale
# create an array of differently colored boxes in the direction of the 3 axes
rectangles = [
    (HyperRectangle(Vec3f0(x), Vec3f0(dirlen, baselen, baselen)), RGBAf0(1,0,0,1)),
    (HyperRectangle(Vec3f0(x), Vec3f0(baselen, dirlen, baselen)), RGBAf0(0,1,0,1)),
    (HyperRectangle(Vec3f0(x), Vec3f0(baselen, baselen, dirlen)), RGBAf0(0,0,1,1))
]
meshes = map(GLNormalMesh, rectangles)

scene = mesh(merge(meshes))
scatter!(scene, X[1,:],X[2,:], X[3,:], markersize = 3, color = :red)
draw!(camera₁, scene)
draw!(camera₂, scene)

𝐑₁′, 𝐭₁′ = ascertain_pose(camera₁, world_basis... )
𝐊₁′ = obtain_intrinsics(camera₁, CartesianSystem())
𝐑₂′, 𝐭₂′ = ascertain_pose(camera₂, world_basis... )
𝐊₂′ = obtain_intrinsics(camera₂, CartesianSystem())

# Camera projection matrices.
𝐏₁ = construct(ProjectionMatrix(),𝐊₁′,𝐑₁′,𝐭₁′)
𝐏₂ = construct(ProjectionMatrix(),𝐊₂′,𝐑₂′,𝐭₂′)

# Set of corresponding points.
ℳ = project(camera₁,𝐏₁,𝒳)
ℳʹ = project(camera₂,𝐏₂,𝒳)

# Estimate of the fundamental matrix and the true fundamental matrix.
𝐅 = estimate(FundamentalMatrix(), DirectLinearTransform(), (ℳ, ℳʹ))

#𝐅ₜ = construct(FundamentalMatrix(), 𝐏₁, 𝐏₂)
𝐅ₜ = construct(FundamentalMatrix(), 𝐊₁′,𝐑₁′,-𝐭₁′,𝐊₂′,𝐑₂′,-𝐭₂′)

# Ensure the estimated and true matrix have the same scale and sign.
𝐅 = 𝐅 / norm(𝐅)
𝐅 = 𝐅 / sign(𝐅[3,1])
𝐅ₜ = 𝐅ₜ / norm(𝐅ₜ)
𝐅ₜ = 𝐅ₜ / sign(𝐅ₜ[3,1])

@test 𝐅 ≈ 𝐅ₜ

# Check that the fundamental matrix satisfies the corresponding point equation.
npts = length(ℳ)
residual = zeros(Float64,npts,1)
for correspondence in zip(1:length(ℳ),ℳ, ℳʹ)
    i, m , mʹ = correspondence
    𝐦  = hom(m)
    𝐦ʹ = hom(mʹ)
    residual[i] = (𝐦ʹ'*𝐅*𝐦)
end

@test isapprox(sum(residual), 0.0; atol = 1e-7)

# Test the Fundamental Numerical Scheme on the Fundamental matrix problem.
Λ₁ =  [SMatrix{3,3}(Matrix(Diagonal([1.0,1.0,0.0]))) for i = 1:length(ℳ)]
Λ₂ =  [SMatrix{3,3}(Matrix(Diagonal([1.0,1.0,0.0]))) for i = 1:length(ℳ)]
𝐅₀ = estimate(FundamentalMatrix(),DirectLinearTransform(),  (ℳ, ℳʹ))
𝐅 = estimate(FundamentalMatrix(),
              FundamentalNumericalScheme(vec(𝐅₀), 5, 1e-10),
               (Λ₁,Λ₂), (ℳ, ℳʹ))

# Ensure the estimated and true matrix have the same scale and sign.
𝐅 = 𝐅 / norm(𝐅)
𝐅 = 𝐅 / sign(𝐅[3,1])

@test 𝐅 ≈ 𝐅ₜ

# Test the Bundle Adjustment estimator on the Fundamental matrix problem.
𝐅, lsqFit = estimate(FundamentalMatrix(),
                      BundleAdjustment(vec(𝐅₀), 5, 1e-10),
                        (ℳ, ℳʹ))
𝐅 = 𝐅 / norm(𝐅)
𝐅 = 𝐅 / sign(𝐅[3,1])
@test 𝐅 ≈ 𝐅ₜ


text:
module DynamicHMCExamples

const src_path = @__DIR__

"Relative path using the DynamicHMCExamples src/ directory."
rel_path(parts...) = normpath(joinpath(src_path, parts...))

end # module


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/iup/include/iupgl.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

const IUP_BUFFER = "BUFFER"
const IUP_STEREO = "STEREO"
const IUP_BUFFER_SIZE = "BUFFER_SIZE"
const IUP_RED_SIZE = "RED_SIZE"
const IUP_GREEN_SIZE = "GREEN_SIZE"
const IUP_BLUE_SIZE = "BLUE_SIZE"
const IUP_ALPHA_SIZE = "ALPHA_SIZE"
const IUP_DEPTH_SIZE = "DEPTH_SIZE"
const IUP_STENCIL_SIZE = "STENCIL_SIZE"
const IUP_ACCUM_RED_SIZE = "ACCUM_RED_SIZE"
const IUP_ACCUM_GREEN_SIZE = "ACCUM_GREEN_SIZE"
const IUP_ACCUM_BLUE_SIZE = "ACCUM_BLUE_SIZE"
const IUP_ACCUM_ALPHA_SIZE = "ACCUM_ALPHA_SIZE"
const IUP_DOUBLE = "DOUBLE"
const IUP_SINGLE = "SINGLE"
const IUP_INDEX = "INDEX"
const IUP_RGBA = "RGBA"
const IUP_YES = "YES"
const IUP_NO = "NO"

function IupGLCanvasOpen()
  ccall( (:IupGLCanvasOpen, iup), Cvoid, (), )
end

function IupGLCanvas()
  ccall( (:IupGLCanvas, iup), Ptr{Cint}, (), )
end

function IupGLMakeCurrent(ih::Ptr{Cint})
  ccall( (:IupGLMakeCurrent, iup), Cvoid, (Ptr{Cint},), ih)
end

function IupGLIsCurrent(ih::Ptr{Cint})
  ccall( (:IupGLIsCurrent, iup), Cint, (Ptr{Cint},), ih)
end

function IupGLSwapBuffers(ih::Ptr{Cint})
  ccall( (:IupGLSwapBuffers, iup), Cvoid, (Ptr{Cint},), ih)
end

function IupGLPalette(ih::Ptr{Cint}, index::Cint, r::Cfloat, g::Cfloat, b::Cfloat)
  ccall( (:IupGLPalette, iup), Cvoid, (Ptr{Cint}, Cint, Cfloat, Cfloat, Cfloat), ih, index, r, g, b)
end

function IupGLUseFont(ih::Ptr{Cint}, first::Cint, count::Cint, list_base::Cint)
  ccall( (:IupGLUseFont, iup), Cvoid, (Ptr{Cint}, Cint, Cint, Cint), ih, first, count, list_base)
end

function IupGLWait(gl::Cint)
  ccall( (:IupGLWait, iup), Cvoid, (Cint,), gl)
end


text:
# Top level test script for Stan.jl
using CmdStan, Test, Statistics

if haskey(ENV, "JULIA_CMDSTAN_HOME")

  println("\nRunning tests for CmdStan-j1-v5:\n")

  code_tests = [
    "test_env.jl",
    "test_utilities.jl",
    "test_cmdtype.jl",
    "test_extract.jl"
  ]

  # Run execution_tests only if cmdstan is installed and CMDSTAN_HOME is set correctly.
  execution_tests = [
    "test_bernoulli.jl",
    "test_bernoulli_diagnostics.jl",
    "test_bernoulli_optimize.jl",
    "test_bernoulli_diagnose.jl",

    # Variational often fails on Travis.

    #"test_bernoulli_variational.jl",
    "test_bernoulliinittheta.jl",
    "test_bernoulliinittheta2.jl",
    "test_bernoulliinittheta3.jl",
    "test_bernoulliscalar.jl",
    "test_binomial.jl",
    "test_binormal.jl",
    "test_schools8.jl",
    "test_dyes.jl",
    "test_kidscore.jl",
    "test_fixed_param.jl",
    "test_zerolengtharray.jl",
    "test_parse_interpolate.jl",
    "test_init_in_stanmodel.jl",
    "test_mcmcchains.jl"
  ]

  if CMDSTAN_HOME != ""
    println("CMDSTAN_HOME set. Try to run tests.")
    @testset "CmdStan.jl" begin
      for my_test in code_tests
          println("\n\n\n  * $(my_test) *")
          include(my_test)
      end

      for my_test in execution_tests
          println("\n\n  * $(my_test) *\n")
          include(my_test)
      end

      println("\n")
    end
  else
    println("\n\nCMDSTAN_HOME not set or found.")
    println("Skipping all tests that depend on CmdStan!\n")
  end

  println("\n")

else
  println("\nJULIA_CMDSTAN_HOME not set. Skipping tests")
end


text:
mutable struct Workspace
    obj::PyObject
end

PyCall.PyObject(ws::Workspace) = getfield(ws, :obj)

function Base.show(io::IO, ws::Workspace)
    println(io, "MATLAB workspace with properties")
    print(io, PyObject(ws))
end

Base.getproperty(ws::Workspace, s::Symbol) = _convert_py_output(PyObject(ws).__getitem__(string(s)))
Base.setproperty!(ws::Workspace, s::Symbol, val) = PyObject(ws).__setitem__(string(s), val)
Base.propertynames(ws::Workspace) = (print(PyObject(ws)); ()) # A hack because reasons

Base.getindex(ws::Workspace, s) = _convert_py_output(PyObject(ws).__getitem__(string(s)))
Base.setindex!(ws::Workspace, val, s) = PyObject(ws).__setitem__(string(s), val)

"""
    eng = connect_matlab(; background=false)
    eng = connect_matlab(engine_name; background=false)

Connect to a named MATLAB session. This must be done before using the integrated MATLAB REPL.
If no name engine name is given, it will connect to the first named session returned from
`find_matlab()`.

## Examples
```julia
julia> using SharedMATLABEngine

julia> connect_matlab(:default_engine);
REPL mode MATLAB initialized. Press > to enter and backspace to exit.

>> a = magic(3)

a =

     8     1     6
     3     5     7
     4     9     2

>> \$a21 = a(2,1)
3.0

julia> a21 .+ mat"a"
3×3 Array{Float64,2}:
 11.0   4.0   9.0
  6.0   8.0  10.0
  7.0  12.0   5.0
```
"""
function connect_matlab(name=nothing; background=false)
    copy!(eng, matlab_engine.connect_matlab(name, background=background))
    setfield!(matlab_workspace, :obj, eng.workspace)
    start_repl()
    return eng
end

"""
    eng = start_matlab(; option="-nodesktop", background=false)

Start new MATLAB session
"""
function start_matlab(; option="-nodesktop", background=false)
    copy!(eng, matlab_engine.start_matlab(option=option, background=background))
    setfield!(matlab_workspace, :obj, eng.workspace)
    start_repl()
    return eng
end

"""
    find_matlab()

Find available MATLAB sessions
"""
find_matlab() = matlab_engine.find_matlab()

function start_repl()
    if isinteractive()
        initrepl(str -> Meta.parse("SharedMATLABEngine.matrepl\"$str\"");
            prompt_text = ">> ",
            start_key = ">",
            prompt_color = :default,
            mode_name = "MATLAB",
        )
    end
    return nothing
end


text:
using Gaston
# using PyCall
using DelimitedFiles

function loadcsv(inputcsv)
    contents = readdlm(inputcsv, ',', Float64, '\n'; header = true)
    return contents
end

function coldata(inputcsv, theset)
    contents = loadcsv(inputcsv)
    return contents[1][:, theset]
end

set(axis="loglog", plotstyle="linespoints", linewidth=2, pointsize = 2, color = "black", xlabel = "Number of nodes", ylabel = "Relative error", grid="on", title = "")

inputcsv = "LE1NAFEMS_MSH8_convergence.CSV"
x = coldata(inputcsv, 1)
y = coldata(inputcsv, 3)
plot(x, abs.(y), legend = "MSOE", marker = "edmd")

inputcsv = "LE1NAFEMS_MSH8_convergence.CSV"
x = coldata(inputcsv, 1)
y = coldata(inputcsv, 2)
plot!(x, abs.(y), legend="TBE", marker = "ecircle")

# inputcsv = "LE1NAFEMS_MST10_convergence.CSV"
# x = coldata(inputcsv, 1)
# y = coldata(inputcsv, 3)
# plot(x, abs.(y), legend = "MSOE", marker = "edmd")

# inputcsv = "LE1NAFEMS_MST10_convergence.CSV"
# x = coldata(inputcsv, 1)
# y = coldata(inputcsv, 2)
# plot!(x, abs.(y), legend="TBE", marker = "ecircle")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Pinchot_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Pinchot")
JLLWrappers.@generate_main_file("Pinchot", UUID("a0e007bb-30f3-5a40-a04b-58baf70e7e14"))
end  # module Pinchot_jll


text:
type Kronos
    data::Array{Number}
    labels::Array{String}
    Markers
    SamplingRate::Number
    Provenance::Dict
end


text:
using LightXML

function toXML(r::JTestResult)
    xroot = new_element("testcase")
    set_attribute(xroot, "id", r.expr)
    set_attribute(xroot, "name", r.expr)
    set_attribute(xroot, "time", string(elapsed_seconds(r)))
    if !r.result
        failure = new_element("failure")
        set_attribute(failure, "message", "Failure")
        set_attribute(failure, "type", "FAIL")
        set_content(failure, "test failed: $(string(r.expr))\nSource: $(r.file):$(r.line)")
        add_child(xroot, failure)
    end
    return xroot
end

function toXML(ts::JUnitTestSet)
    xroot = typeof(ts.results[1]) == JUnitTestSet ? new_element("testsuites") : new_element("testsuite")
    set_attribute(xroot, "id", ts.id)
    set_attribute(xroot, "name", ts.name)
    set_attribute(xroot, "time", ts.time.value / 1e9)

    for result in ts.results
        childnode = toXML(result)
        add_child(xroot, childnode)
    end

    set_attribute(xroot, "tests", ts.tests)
    set_attribute(xroot, "failures", ts.failures)
    return xroot
end

function write_report(ts::JUnitTestSet)
    el = toXML(ts)
    xroot = XMLDocument()
    set_root(xroot, el)
    return xroot
end