text:
@testset "ThermalGen data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    warn_message = "The data doesn't include devices of type ThermalStandard, consider changing the device models"
    model = DeviceModel(ThermalStandard, ThermalStandardUnitCommitment)
    c_sys5_re_only = build_system("c_sys5_re_only")
    op_problem =
        OperationsProblem(TestOpProblem, DCPPowerModel, build_system("c_sys5_re_only"))
    @test_logs (:warn, warn_message) construct_device!(op_problem, :Thermal, model)
end

################################### Unit Commitment tests ##################################
@testset "Thermal UC With DC - PF" begin
    bin_variable_names = [
        PSI.make_variable_name(PSI.ON, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.START, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.STOP, PSY.ThermalStandard),
    ]
    uc_constraint_names = [
        PSI.make_constraint_name(PSI.RAMP_UP, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.RAMP_DOWN, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.DURATION_UP, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.DURATION_DOWN, PSY.ThermalStandard),
    ]
    model = DeviceModel(ThermalStandard, ThermalStandardUnitCommitment)

    @info "5-Bus testing"
    c_sys5_uc = build_system("c_sys5_uc")
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_uc)
    construct_device!(op_problem, :Thermal, model)
    moi_tests(op_problem, false, 480, 0, 480, 120, 120, true)
    psi_constraint_test(op_problem, uc_constraint_names)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    op_problem =
        OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_uc; use_parameters = true)
    construct_device!(op_problem, :Thermal, model)
    moi_tests(op_problem, true, 480, 0, 480, 120, 120, true)
    psi_constraint_test(op_problem, uc_constraint_names)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    @info "14-Bus testing"
    c_sys14 = build_system("c_sys14")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, DCPPowerModel, c_sys14; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 480, 0, 240, 120, 120, true)
        psi_checkbinvar_test(op_problem, bin_variable_names)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal UC With AC - PF" begin
    bin_variable_names = [
        PSI.make_variable_name(PSI.ON, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.START, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.STOP, PSY.ThermalStandard),
    ]
    uc_constraint_names = [
        PSI.make_constraint_name(PSI.RAMP_UP, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.RAMP_DOWN, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.DURATION_UP, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.DURATION_DOWN, PSY.ThermalStandard),
    ]
    model = DeviceModel(ThermalStandard, ThermalStandardUnitCommitment)

    @info "5-Bus testing"
    c_sys5_uc = build_system("c_sys5_uc")
    op_problem = OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5_uc)
    construct_device!(op_problem, :Thermal, model)
    moi_tests(op_problem, false, 600, 0, 600, 240, 120, true)
    psi_constraint_test(op_problem, uc_constraint_names)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    op_problem =
        OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5_uc; use_parameters = true)
    construct_device!(op_problem, :Thermal, model)
    moi_tests(op_problem, true, 600, 0, 600, 240, 120, true)
    psi_constraint_test(op_problem, uc_constraint_names)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    @info "14-Bus testing"
    c_sys14 = build_system("c_sys14")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, ACPPowerModel, c_sys14; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 600, 0, 360, 240, 120, true)
        psi_checkbinvar_test(op_problem, bin_variable_names)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

################################### Basic Unit Commitment tests ############################
@testset "Thermal Basic UC With DC - PF" begin
    bin_variable_names = [
        PSI.make_variable_name(PSI.ON, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.START, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.STOP, PSY.ThermalStandard),
    ]
    model = DeviceModel(ThermalStandard, ThermalBasicUnitCommitment)

    @info "5-Bus testing"
    c_sys5_uc = build_system("c_sys5_uc")
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_uc)
    construct_device!(op_problem, :Thermal, model)
    moi_tests(op_problem, false, 480, 0, 240, 120, 120, true)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    op_problem =
        OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_uc; use_parameters = true)
    construct_device!(op_problem, :Thermal, model)
    moi_tests(op_problem, true, 480, 0, 240, 120, 120, true)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    @info "14-Bus testing"
    c_sys14 = build_system("c_sys14")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, DCPPowerModel, c_sys14; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 480, 0, 240, 120, 120, true)
        psi_checkbinvar_test(op_problem, bin_variable_names)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal Basic UC With AC - PF" begin
    bin_variable_names = [
        PSI.make_variable_name(PSI.ON, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.START, PSY.ThermalStandard),
        PSI.make_variable_name(PSI.STOP, PSY.ThermalStandard),
    ]
    model = DeviceModel(ThermalStandard, ThermalBasicUnitCommitment)

    @info "5-Bus testing"
    c_sys5_uc = build_system("c_sys5_uc")
    op_problem = OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5_uc)
    construct_device!(op_problem, :Thermal, model)
    moi_tests(op_problem, false, 600, 0, 360, 240, 120, true)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    op_problem =
        OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5_uc; use_parameters = true)
    construct_device!(op_problem, :Thermal, model)
    moi_tests(op_problem, true, 600, 0, 360, 240, 120, true)
    psi_checkbinvar_test(op_problem, bin_variable_names)
    psi_checkobjfun_test(op_problem, GAEVF)

    @info "14-Bus testing"
    c_sys14 = build_system("c_sys14")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, ACPPowerModel, c_sys14; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 600, 0, 360, 240, 120, true)
        psi_checkbinvar_test(op_problem, bin_variable_names)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

################################### Basic Dispatch tests ###################################
@testset "Thermal Dispatch With DC - PF" begin
    model = DeviceModel(ThermalStandard, ThermalDispatch)
    @info "5-Bus testing"
    c_sys5 = build_system("c_sys5")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 120, 0, 120, 120, 0, false)
        psi_checkobjfun_test(op_problem, GAEVF)
    end

    @info "14-Bus testing"
    c_sys14 = build_system("c_sys14")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, DCPPowerModel, c_sys14; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 120, 0, 120, 120, 0, false)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal Dispatch With AC - PF" begin
    model = DeviceModel(ThermalStandard, ThermalDispatch)
    @info "5-Bus testing"
    c_sys5 = build_system("c_sys5")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 240, 0, 240, 240, 0, false)
        psi_checkobjfun_test(op_problem, GAEVF)
    end

    @info "14-Bus testing"
    c_sys14 = build_system("c_sys14")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, ACPPowerModel, c_sys14; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 240, 0, 240, 240, 0, false)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

################################### No Minimum Dispatch tests ##############################

@testset "Thermal Dispatch NoMin With DC - PF" begin
    model = DeviceModel(ThermalStandard, ThermalDispatchNoMin)
    @info "5-Bus testing"
    c_sys5 = build_system("c_sys5")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 120, 0, 120, 120, 0, false)
        moi_lbvalue_test(op_problem, :P_lb__ThermalStandard__RangeConstraint, 0.0)
        psi_checkobjfun_test(op_problem, GAEVF)
    end

    @info "14-Bus testing"
    c_sys14 = build_system("c_sys14")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, DCPPowerModel, c_sys14; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 120, 0, 120, 120, 0, false)
        moi_lbvalue_test(op_problem, :P_lb__ThermalStandard__RangeConstraint, 0.0)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal Dispatch NoMin With AC - PF" begin
    model = DeviceModel(ThermalStandard, ThermalDispatchNoMin)
    @info "5-Bus testing"
    c_sys5 = build_system("c_sys5")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 240, 0, 240, 240, 0, false)
        moi_lbvalue_test(op_problem, :P_lb__ThermalStandard__RangeConstraint, 0.0)
        psi_checkobjfun_test(op_problem, GAEVF)
    end

    @info "14-Bus testing"
    c_sys14 = build_system("c_sys14")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, ACPPowerModel, c_sys14; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 240, 0, 240, 240, 0, false)
        moi_lbvalue_test(op_problem, :P_lb__ThermalStandard__RangeConstraint, 0.0)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

################################### Ramp Limited Testing ##################################
@testset "Thermal Ramp Limited Dispatch With DC - PF" begin
    constraint_names = [
        PSI.make_constraint_name(PSI.RAMP_UP, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.RAMP_DOWN, PSY.ThermalStandard),
    ]
    model = DeviceModel(ThermalStandard, ThermalRampLimited)
    @info "5-Bus testing"
    c_sys5_uc = build_system("c_sys5_uc")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_uc; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 120, 0, 216, 120, 0, false)
        psi_constraint_test(op_problem, constraint_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end

    @info "14-Bus testing"
    c_sys14 = build_system("c_sys14")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, DCPPowerModel, c_sys14; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 120, 0, 120, 120, 0, false)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

@testset "Thermal Ramp Limited Dispatch With AC - PF" begin
    constraint_names = [
        PSI.make_constraint_name(PSI.RAMP_UP, PSY.ThermalStandard),
        PSI.make_constraint_name(PSI.RAMP_DOWN, PSY.ThermalStandard),
    ]
    model = DeviceModel(ThermalStandard, ThermalRampLimited)
    @info "5-Bus testing"
    c_sys5_uc = build_system("c_sys5_uc")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, ACPPowerModel, c_sys5_uc; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 240, 0, 336, 240, 0, false)
        psi_constraint_test(op_problem, constraint_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end

    @info "14-Bus testing"
    c_sys14 = build_system("c_sys14")
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, ACPPowerModel, c_sys14; use_parameters = p)
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 240, 0, 240, 240, 0, false)
        psi_checkobjfun_test(op_problem, GQEVF)
    end
end

################################### ThermalMultiStart Testing ##################################

@testset "Thermal MultiStart Dispatch With DC - PF" begin
    constraint_names = [
        PSI.make_constraint_name(PSI.ACTIVE_RANGE_IC, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.START_TYPE, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.MUST_RUN_LB, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.STARTUP_TIMELIMIT_WARM, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.STARTUP_TIMELIMIT_HOT, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.STARTUP_INITIAL_CONDITION_LB, PSY.ThermalMultiStart),
        PSI.make_constraint_name(PSI.STARTUP_INITIAL_CONDITION_UB, PSY.ThermalMultiStart),
    ]
    model = DeviceModel(PSY.ThermalMultiStart, PSI.ThermalMultiStartUnitCommitment)
    no_less_than = Dict(true => 238, false => 234)
    @info "5-Bus testing"
    c_sys5_pglib = build_system("c_sys5_pglib")
    for p in [true, false]
        op_problem = OperationsProblem(
            TestOpProblem,
            DCPPowerModel,
            c_sys5_pglib;
            use_parameters = p,
        )
        construct_device!(op_problem, :Thermal, model)
        moi_tests(op_problem, p, 528, 0, no_less_than[p], 60, 192, true)
        psi_constraint_test(op_problem, constraint_names)
        psi_checkobjfun_test(op_problem, GAEVF)
    end
end

############################# UC validation tests ##########################################
branches = Dict{Symbol, DeviceModel}()
services = Dict{Symbol, ServiceModel}()
ED_devices = Dict{Symbol, DeviceModel}(
    :Generators => DeviceModel(ThermalStandard, ThermalRampLimited),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
)
UC_devices = Dict{Symbol, DeviceModel}(
    :Generators => DeviceModel(ThermalStandard, ThermalStandardUnitCommitment),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
)
# Testing Ramping Constraint
@testset "Solving ED with CopperPlate for testing Ramping Constraints" begin
    ramp_test_sys = build_system("c_ramp_test")
    template =
        OperationsProblemTemplate(CopperPlatePowerModel, ED_devices, branches, services)
    for p in [true, false]
        ED = OperationsProblem(
            TestOpProblem,
            template,
            ramp_test_sys;
            optimizer = Cbc_optimizer,
            use_parameters = p,
        )
        psi_checksolve_test(ED, [MOI.OPTIMAL], 11191.00)
        moi_tests(ED, p, 10, 0, 20, 10, 5, false)
    end
end

# Testing Duration Constraints
@testset "Solving UC with CopperPlate for testing Duration Constraints" begin
    node = Bus(1, "nodeA", "PV", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing)
    load = PowerLoad("Bus1", true, node, nothing, 0.4, 0.9861, 100.0, 1.0, 2.0)
    DA_dur = collect(
        DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):Hour(1):DateTime(
            "1/1/2024  6:00:00",
            "d/m/y  H:M:S",
        ),
    )
    gens_dur = [
        ThermalStandard(
            name = "Alta",
            available = true,
            status = true,
            bus = node,
            active_power = 0.40,
            reactive_power = 0.010,
            rating = 0.5,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.3, max = 0.9),
            reactive_power_limits = nothing,
            ramp_limits = nothing,
            time_limits = (up = 4, down = 2),
            operation_cost = ThreePartCost((0.0, 1400.0), 0.0, 4.0, 2.0),
            base_power = 100.0,
            time_at_status = 2.0,
        ),
        ThermalStandard(
            name = "Park City",
            available = true,
            status = false,
            bus = node,
            active_power = 1.70,
            reactive_power = 0.20,
            rating = 2.2125,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.7, max = 2.2),
            reactive_power_limits = nothing,
            ramp_limits = nothing,
            time_limits = (up = 6, down = 4),
            operation_cost = ThreePartCost((0.0, 1500.0), 0.0, 1.5, 0.75),
            base_power = 100.0,
            time_at_status = 3.0,
        ),
    ]

    duration_load = [0.3, 0.6, 0.8, 0.7, 1.7, 0.9, 0.7]
    load_data = SortedDict(DA_dur[1] => TimeArray(DA_dur, duration_load))
    load_forecast_dur = Deterministic("max_active_power", load_data)
    duration_test_sys = System(100.0)
    add_component!(duration_test_sys, node)
    add_component!(duration_test_sys, load)
    add_component!(duration_test_sys, gens_dur[1])
    add_component!(duration_test_sys, gens_dur[2])
    add_time_series!(duration_test_sys, load, load_forecast_dur)

    template =
        OperationsProblemTemplate(CopperPlatePowerModel, UC_devices, branches, services)
    UC = OperationsProblem(
        TestOpProblem,
        template,
        duration_test_sys;
        optimizer = Cbc_optimizer,
        use_parameters = true,
    )
    psi_checksolve_test(UC, [MOI.OPTIMAL], 8223.50)
    moi_tests(UC, true, 56, 0, 56, 14, 21, true)
end

## PWL linear Cost implementation test
@testset "Solving UC with CopperPlate testing Linear PWL" begin
    node = Bus(1, "nodeA", "PV", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing)
    load = PowerLoad("Bus1", true, node, nothing, 0.4, 0.9861, 100.0, 1.0, 2.0)
    gens_cost = [
        ThermalStandard(
            name = "Alta",
            available = true,
            status = true,
            bus = node,
            active_power = 0.52,
            reactive_power = 0.010,
            rating = 0.5,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.22, max = 0.55),
            reactive_power_limits = nothing,
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = ThreePartCost(
                [(589.99, 0.220), (884.99, 0.33), (1210.04, 0.44), (1543.44, 0.55)],
                532.44,
                5665.23,
                0.0,
            ),
            base_power = 100.0,
        ),
        ThermalStandard(
            name = "Park City",
            available = true,
            status = true,
            bus = node,
            active_power = 0.62,
            reactive_power = 0.20,
            rating = 221.25,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.62, max = 1.55),
            reactive_power_limits = nothing,
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = ThreePartCost(
                [(1264.80, 0.62), (1897.20, 0.93), (2594.4787, 1.24), (3433.04, 1.55)],
                235.397,
                5665.23,
                0.0,
            ),
            base_power = 100.0,
        ),
    ]

    DA_load_forecast = SortedDict{Dates.DateTime, TimeArray}()
    ini_time = DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S")
    cost_sos_load = [[1.3, 2.1], [1.3, 2.1]]
    for (ix, date) in enumerate(range(ini_time; length = 2, step = Hour(1)))
        DA_load_forecast[date] = TimeArray([date, date + Hour(1)], cost_sos_load[ix])
    end
    load_forecast_cost_sos = Deterministic("max_active_power", DA_load_forecast)
    cost_test_sys = System(100.0)
    add_component!(cost_test_sys, node)
    add_component!(cost_test_sys, load)
    add_component!(cost_test_sys, gens_cost[1])
    add_component!(cost_test_sys, gens_cost[2])
    add_time_series!(cost_test_sys, load, load_forecast_cost_sos)

    template =
        OperationsProblemTemplate(CopperPlatePowerModel, UC_devices, branches, services)
    UC = OperationsProblem(
        TestOpProblem,
        template,
        cost_test_sys;
        optimizer = Cbc_optimizer,
        use_parameters = true,
    )
    psi_checksolve_test(UC, [MOI.OPTIMAL], 9336.736919354838)
    moi_tests(UC, true, 32, 0, 8, 4, 10, true)
end

## PWL SOS-2 Cost implementation test
@testset "Solving UC with CopperPlate testing SOS2 implementation" begin
    node = Bus(1, "nodeA", "PV", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing)
    load = PowerLoad("Bus1", true, node, nothing, 0.4, 0.9861, 100.0, 1.0, 2.0)
    gens_cost_sos = [
        ThermalStandard(
            name = "Alta",
            available = true,
            status = true,
            bus = node,
            active_power = 0.52,
            reactive_power = 0.010,
            rating = 0.5,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.22, max = 0.55),
            reactive_power_limits = nothing,
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = ThreePartCost(
                [(1122.43, 0.22), (1617.43, 0.33), (1742.48, 0.44), (2075.88, 0.55)],
                0.0,
                5665.23,
                0.0,
            ),
            base_power = 100.0,
        ),
        ThermalStandard(
            name = "Park City",
            available = true,
            status = true,
            bus = node,
            active_power = 0.62,
            reactive_power = 0.20,
            rating = 2.2125,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.62, max = 1.55),
            reactive_power_limits = nothing,
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = ThreePartCost(
                [(1500.19, 0.62), (2132.59, 0.929), (2829.875, 1.24), (2831.444, 1.55)],
                0.0,
                5665.23,
                0.0,
            ),
            base_power = 100.0,
        ),
    ]
    DA_load_forecast = SortedDict{Dates.DateTime, TimeArray}()
    ini_time = DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S")
    cost_sos_load = [[1.3, 2.1], [1.3, 2.1]]
    for (ix, date) in enumerate(range(ini_time; length = 2, step = Hour(1)))
        DA_load_forecast[date] = TimeArray([date, date + Hour(1)], cost_sos_load[ix])
    end
    load_forecast_cost_sos = Deterministic("max_active_power", DA_load_forecast)
    cost_test_sos_sys = System(100.0)
    add_component!(cost_test_sos_sys, node)
    add_component!(cost_test_sos_sys, load)
    add_component!(cost_test_sos_sys, gens_cost_sos[1])
    add_component!(cost_test_sos_sys, gens_cost_sos[2])
    add_time_series!(cost_test_sos_sys, load, load_forecast_cost_sos)

    for g in gens_cost_sos
        @test PSI.pwlparamcheck(PSY.get_operation_cost(g).variable) == false
    end

    template =
        OperationsProblemTemplate(CopperPlatePowerModel, UC_devices, branches, services)
    UC = OperationsProblem(
        TestOpProblem,
        template,
        cost_test_sos_sys;
        optimizer = Cbc_optimizer,
        use_parameters = true,
    )
    psi_checksolve_test(UC, [MOI.OPTIMAL], 8500.89716, 10.0)
    moi_tests(UC, true, 32, 0, 8, 4, 14, true)
end


text:
module TestDNS

using LibCloud.DNS
using Base.Test
using PyCall

const TEST_DOMAIN = "julialibcloudtest.com."

function common(dns, upd=false)
    println("Regions:\n", list_regions(dns))

    println("Record Types:\n", list_record_types(dns))

    println("Creating test zone ($TEST_DOMAIN)...")
    zone = create_zone(dns, TEST_DOMAIN)
    @test isa(zone, Zone)

    if upd
        println("Updating test zone ($TEST_DOMAIN)...")
        zone = update_zone(dns, zone, TEST_DOMAIN, ttl=600)
        @test isa(zone, Zone)
    end

    println("Creating test records...")

    data = Dict{Any,Any}()
    data["rrdatas"] = ["127.0.0.1"]

    rec1 = create_record(dns, "www."*TEST_DOMAIN, zone, RecordType.A, data)
    @test isa(rec1, Record)
    rec2 = create_record(dns, "ftp."*TEST_DOMAIN, zone, RecordType.A, data)
    @test isa(rec2, Record)

    println("Listing zones...")
    for z in list_zones(dns)
        @test isa(z, Zone)
        println(z)
        for r in list_records(z)
            @test isa(r, Record)
            println("    - ", r)
        end
    end

    println("Deleting test records...")
    for rec in list_records(zone)
        if rec.typ == RecordType.A
            println("    - ", rec)
            @test delete_record(dns, rec)
        end
    end

    println("Deleting test zone...")
    @test delete_zone(dns, zone)
end

function dummy()
    dns = DNSDriver(DNSProvider.DUMMY, "key", "secret")
    common(dns)
end

function route53(key, secret)
    dns = DNSDriver(DNSProvider.ROUTE53, key, secret)
    common(dns)
end

function gce(service_account_id, key, project)
	dns = DNSDriver(DNSProvider.GOOGLE, service_account_id, key, project=project)
    common(dns)
end

end # module TestDNS

ext: jl


text:
module Clusters

using DataFrames, GLM, CovarianceMatrices, Distributions, GrowableArrays

abstract type WildWeights end
struct Rademacher <: WildWeights end
struct TwoPoint <: WildWeights end
struct Mammen1 <: WildWeights end
struct Mammen2 <: WildWeights end

abstract type MonteCarloModel end
abstract type MonteCarloModelOpt end

struct LinearRegressionCluster{G} <: MonteCarloModel
    ## Containers
    y::Array{Float64, 1}
    X::Array{Float64, 2}
    epsilon::Array{Float64, 1}
    eta::Array{Float64, 1}
    D::Array{Float64, 1}
    sigma_e::Array{Float64, 1}
    x̄::Array{Float64, 1}
    cl::Array{Int64, 1}
    iter::Array{Int64, 1}
    bstarts::Array{UnitRange{Int64},1}
    wts::Array{Float64, 1}
    sqwts::Array{Float64, 1}
    opt::G
end

struct LinearRegressionClusterOpt <: MonteCarloModelOpt
    G::Int64
    ng::Vector{Int64}
    ## Regressors std. dev.
    σ_z::Float64
    σ_ξ::Float64
    ## Errors std. dev
    σ_ϵ::Float64
    σ_α::Float64
    ## Third/Fourth design parms
    p::Float64
    γ::Float64
    δ::Float64
    σ_η::Float64
    b::Float64
    ## Null value
    β₀::Float64
    design::Int
end

LinearRegressionClusterOpt(G::Int64,
                           ng::Vector{Int64},
                           σ_z::Float64,
                           σ_ξ::Float64,
                           σ_ϵ::Float64,
                           σ_α::Float64,
                           p::Float64,
                           γ::Float64,
                           δ::Float64,
                           σ_η::Float64,
                           b::Float64,
                           β₀::Float64,
                           design::T) where T<:Union{Float16, Float32, Float64, Int} =
                           LinearRegressionClusterOpt(G, ng, σ_z, σ_ξ, σ_ϵ, σ_α, p,
                                                      γ, δ, σ_η, b, β₀, convert(Int, design))

function dispersion_ng(m::LinearRegressionClusterOpt)
    sum(m.ng^2)/sum(ng)^2
end

function icc_u(m::LinearRegressionClusterOpt, ::Type{Val{1}})
    σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, b = getparms(m)
    σ_α^2/(σ_α^2+σ_ϵ^2)
end

function icc_u(m::LinearRegressionClusterOpt, ::Type{Val{2}})
    σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, b = getparms(m)
    σ_α^2/(σ_α^2+σ_ϵ^2)
end

function icc_u(m::LinearRegressionClusterOpt, ::Type{Val{3}})
    σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, b = getparms(m)
    (σ_α^2 + (1-p)*p*γ^2*σ_ξ^2)/(σ_α^2 + σ_ϵ^2 + (1-p)*p*γ^2*(1+σ_ξ^2))
end

function icc_u(m::LinearRegressionClusterOpt, ::Type{Val{4}})
    σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, b = getparms(m)
    num = γ^2*σ_ξ^2*(δ^2*σ_z^2 + σ_η^2 + δ*(2*δ + 3*b)*σ_ξ^2)
    den = 1 + γ^2*(σ_z^2*(b^2 + σ_η^2) + (b^2 + δ^2*σ_z^2 + σ_η^2)*σ_ξ^2 + 2*δ^2*σ_ξ^4)
    num/den
end

function icc_xu(m::LinearRegressionClusterOpt)
    icc_xu(m, Val{convert(Int, m.design)})
end

function icc_xu(m::LinearRegressionClusterOpt, ::Type{Val{1}})
    icc_u(m)*icc_x(m)
end

function icc_xu(m::LinearRegressionClusterOpt, ::Type{Val{2}})
    σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, b = getparms(m)
    (σ_α^2*σ_ξ^2)/( (3*b^2*σ_ϵ^2)*(σ_ξ^4+2*σ_ξ^2*σ_z^2+σ_z^4)+σ_α^2*(σ_ξ^2+σ_z^2) )
end

function icc_xu(m::LinearRegressionClusterOpt, ::Type{Val{3}})
    σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, b = getparms(m)
    num = 3*γ^2*(p-1)*p*σ_ξ^4 + σ_ξ^2*(2*γ^2*(p-1)*p*σ_z^2-σ_α^2) + (p-1)*p*γ^2*σ_z^4
    den = (σ_ξ^2+σ_z^2)*(σ_α^2-3*γ^2*(p-1)*p*(σ_ξ^2+σ_z^2)+σ_ϵ^2)
    -num/den
end

function icc_xu(m::LinearRegressionClusterOpt, ::Type{Val{4}})
    σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, b = getparms(m)
    num = γ^2*(3*σ_ξ^4*(5*δ^2*σ_ξ^2+σ_η^2)+σ_z^4*(δ^2*σ_ξ^2+σ_η^2)+2σ_ξ^2*σ_z^2*(3*δ^2*σ_ξ^2+σ_η^2))
    den = (σ_ξ^2+σ_z^2)*(3*γ^2*(5*δ^2*σ_ξ^4+σ_ξ^2*(σ_η^2+δ^2*σ_z^2)+σ_η^2*σ_z^2)+σ_ϵ^2)
    num/den
end

function icc_x(m::LinearRegressionClusterOpt)
    σ_z = m.σ_z
    σ_ξ = m.σ_ξ
    σ_ξ^2/(σ_ξ^2+σ_z^2)
end

function icc_u(m::LinearRegressionClusterOpt)
    idx = convert(Int, m.design)
    icc_u(m, Val{idx})
end

icc_x(m::LinearRegressionCluster) = icc_x(m.opt)
icc_u(m::LinearRegressionCluster) = icc_u(m.opt)

struct ConfInterval
    lo::Array{Float64,1}
    hi::Array{Float64,1}
end

function confinterval(m::Real, s::Real, zl::Real, zh::Real)
    m = float(m)
    (m+zl*s, m+zh*s)
end

function ConfInterval(m::Array{Float64}, σ::Array{Float64}, z::Float64 = 1.96)
    zs = z.*σ
    ConfInterval(m - zs, m + zs)
end

function ConfInterval(m::Array{Float64}, σ::Array{Float64},
    zl::Array{Float64}, zh::Array{Float64})
    @assert length(m) == length(σ)
    @assert length(m) == length(zh)
    @assert length(zl) == length(zh)

    lo = m .+ zl.*σ
    hi = m .+ zh.*σ

    @assert all(hi.>=lo)
    ConfInterval(m .+ zl.*σ, m .+ zh.*σ)
end

Base.length(i::ConfInterval) = length(i.lo)

function average_length(i::Clusters.ConfInterval)
    l = 0.0
    n = length(i)
    @simd for j in 1:n
        @inbounds l += i.hi[j]-i.lo[j]
    end
    l/n
end

function Base.in(i::ConfInterval, a::Number)
    isin = BitArray(length(i))
    for j in 1:length(i)
        if (i.lo[j] <= a && i.hi[j] >= a)
            isin[j] = true
        else
            isin[j] = false
        end
    end
    isin
end

function coverage(i::Clusters.ConfInterval, a::Number)
    ifin = in(i, a)
    mean(ifin)
end

function initialize(::Type{LinearRegressionCluster}, opt)
    G = opt.G
    ng = opt.ng
    n = sum(ng)
    cl = vcat(map((i, n) -> repeat(i:i, outer = n), 1:G, ng)...)
    wts = vcat(map((i, n) -> repeat(1/n:1/n, outer = n), 1:G, ng)...)
    sqwts = sqrt.(wts)
    iter = mapreduce(x -> collect(1:x), vcat, ng)
    bstarts = [searchsorted(cl, j[2]) for j in enumerate(unique(cl))]
    y = Array{Float64}(n)
    X = Array{Float64}(n, 1)
    epsilon = similar(y)
    eta = Array{Float64}(G)
    D = similar(eta)
    sigma_e = similar(y)
    x̄ = similar(eta)
    LinearRegressionCluster(y, X, epsilon, eta, D, sigma_e, x̄, cl, iter, bstarts, wts, sqwts, opt)
end

function getcontainers(m::LinearRegressionCluster)
    (m.X, m.y, m.epsilon, m.eta, m.D, m.sigma_e, m.x̄)
end

getparms(m::LinearRegressionCluster) = getparms(m.opt)

function getparms(opt::LinearRegressionClusterOpt)
    (opt.σ_z,
    opt.σ_ξ,
    opt.σ_ϵ,
    opt.σ_α,
    opt.p,
    opt.γ,
    opt.δ,
    opt.σ_η,
    opt.b)
end

function gethyper(m::LinearRegressionCluster)
    (m.cl, m.iter, m.bstarts)
end

function simulate!(m::LinearRegressionCluster, ::Type{Val{1}})
    X, y, ϵ, η, D, σ_e, x̄ = getcontainers(m)
    σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, b = getparms(m)
    clus, iter, bstarts = gethyper(m)
    ## GENERATE X = z + ξ
    randn!(X)
    scale!(X, σ_z)
    randn!(η)
    scale!(η, σ_ξ)
    for (i, (n, g)) in enumerate(zip(iter, clus))
        X[i] += η[g]
    end
    ## GENERATE u = α + ξ
    randn!(ϵ)
    scale!(ϵ, σ_ϵ)
    randn!(η)
    scale!(η, σ_α)
    for (i, (n, g)) in enumerate(zip(iter, clus))
        ϵ[i] += η[g]
    end
end

function simulate!(m::LinearRegressionCluster, ::Type{Val{2}})
    X, y, ϵ, η, D, σ_e, x̄ = getcontainers(m)
    σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, b = getparms(m)
    clus, iter, bstarts = gethyper(m)

    ## GENERATE X = z + ξ
    randn!(X)
    scale!(X, σ_z)
    randn!(η)
    scale!(η, σ_ξ)
    for (i, (n, g)) in enumerate(zip(iter, clus))
        X[i] += η[g]
    end

    ## GENERATE u = α + ξ
    randn!(ϵ)
    for i in eachindex(ϵ)
        ϵ[i] *= abs(X[i])*b*σ_ϵ
    end
    randn!(η)
    scale!(η, σ_α)
    for (i, (n, g)) in enumerate(zip(iter, clus))
        ϵ[i] += η[g]
    end
end

function simulate!(m::LinearRegressionCluster, ::Type{Val{3}})
    X, y, ϵ, η, D, σ_e, x̄ = getcontainers(m)
    σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, b = getparms(m)
    clus, iter, bstarts = gethyper(m)

    ## GENERATE X = z + ξ
    randn!(X)
    scale!(X, σ_z)
    randn!(η)
    scale!(η, σ_ξ)
    for (i, (n, g)) in enumerate(zip(iter, clus))
        X[i] += η[g]
    end

    ## GENERATE u = α + ξ
    randn!(ϵ)
    scale!(ϵ, σ_ϵ)
    randn!(η)
    scale!(η, σ_α)
    for (i, (n, g)) in enumerate(zip(iter, clus))
        ϵ[i] += η[g]
    end

    rand!(η)
    for i in eachindex(D)
        D[i] = ifelse(η[i] .<= p, 1-p, -p)
    end
    for (i, (n, g)) in enumerate(zip(iter, clus))
        ϵ[i] += γ*X[i]*D[g]
    end
end

function simulate!(m::LinearRegressionCluster, ::Type{Val{4}})
    X, y, ϵ, η, D, σ_e, x̄ = getcontainers(m)
    σ_z, σ_ξ, σ_ϵ, σ_α, p, γ, δ, σ_η, b = getparms(m)
    clus, iter, bstarts = gethyper(m)
    ## GENERATE X = z + ξ
    randn!(X)
    scale!(X, σ_z)
    randn!(η)
    scale!(η, σ_ξ)
    for (i, (n, g)) in enumerate(zip(iter, clus))
        X[i] += η[g]
    end
    ## GENERATE u = α + ξ
    randn!(η)
    scale!(η, σ_η)
    for (i, j) in enumerate(bstarts)
        x̄[i] = mean(X[j])
        η[i] += δ*x̄[i]
    end

    randn!(ϵ)
    scale!(ϵ, σ_ϵ)
    for (i, (n, g)) in enumerate(zip(iter, clus))
        ϵ[i] += γ*X[i]*η[g]
    end

end

function simulate!(m::LinearRegressionCluster{T}) where T
    X, y, ϵ, _ = getcontainers(m)
    β₀ = m.opt.β₀
    design = m.opt.design
    simulate!(m, Val{design})
    ## Change this when the power is sought
    @simd for i in eachindex(y)
        @inbounds y[i] = X[i].*β₀ + ϵ[i]
    end
end

function montecarlo(m::Type{T} where T <: MonteCarloModel, opt::MonteCarloModelOpt; simulations::Int64 = 1000)
    model = initialize(m, opt)
    simulate!(model)
    out = estimatemodel(model)
    res = GrowableArray(out)
    sizehint!(res,simulations)
    @inbounds for j in 1:simulations
        simulate!(model)
        push!(res, estimatemodel(model))
    end
    ## Add information about the experiment

    res = convert(DataFrame, convert(Array,res))
    names!(res,
    [:theta_u, :V1_u, :V2_u, :V3_u, :V4_u, :V5_u,
    :G_u, :k_u,
    :qu_025, :qu_975,
    :qu_050, :qu_950,
    :qu0_025, :qu0_975,
    :qu0_050, :qu0_950,
    :theta_w, :V1_w, :V2_w, :V3_w, :V4_w, :V5_w,
    :G_w, :k_w,
    :qw0_025, :qw0_975, :qw0_050, :qw0_950])

end

function montecarlolight(m::Type{T} where T <: MonteCarloModel, opt::MonteCarloModelOpt; simulations::Int64 = 1000)
    model = initialize(m, opt)
    simulate!(model)
    out = estimatemodellight(model)
    res = GrowableArray(out)
    sizehint!(res,simulations)
    @inbounds for j in 1:simulations
        simulate!(model)
        push!(res, estimatemodellight(model))
    end
    ## Add information about the experiment

    res = convert(DataFrame, convert(Array,res))
    names!(res,
    [:theta_u, :V1_u, :V2_u, :V3_u, :V4_u, :V5_u,
    :G_u, :k_u,
    :qu_025, :qu_975,
    :qu_050, :qu_950,
    :qu0_025, :qu0_975,
    :qu0_050, :qu0_950,
    :theta_w, :V1_w, :V2_w, :V3_w, :V4_w, :V5_w,
    :G_w, :k_w,
    :qw0_025, :qw0_975, :qw0_050, :qw0_950])

end

function estimatemodel(m::LinearRegressionCluster)
    y = m.y
    X = m.X
    cl = m.cl
    iter = m.iter
    w  = Val{:weighted}
    uw = Val{:unweighted}
    try
        fitted_u = fit(GeneralizedLinearModel, X, y, Normal(), IdentityLink())

        theta_u = first(coef(fitted_u))
        V1_u = first(stderr(fitted_u))
        V2_u = first(stderr(fitted_u, HC1()))

        V3_u = faststderr(fitted_u, m, CRHC1(cl), uw)
        V4_u = faststderr(fitted_u, m, CRHC2(cl), uw)
        V5_u = faststderr(fitted_u, m, CRHC3(cl), uw)

        fitted_w = fit(GeneralizedLinearModel, X, y, Normal(), IdentityLink(), wts = m.wts)
        theta_w = first(coef(fitted_w))
        V1_w = fastiid(fitted_w, m)
        V2_w = first(stderr(fitted_w, HC1()))

        V3_w = faststderr(fitted_w, m, CRHC1(cl), w)
        V4_w = faststderr(fitted_w, m, CRHC2(cl), w)
        V5_w = faststderr(fitted_w, m, CRHC3(cl), w)

        G_u, k_u = kappastar(fitted_u, m, uw)
        G_w, k_w = kappastar(fitted_w, m, w)

        _, _, qw = fastwildboot_nonull(fitted_w, m, Rademacher(), rep = 499)
        _, _, qu = fastwildboot_nonull(fitted_u, m, Rademacher(), rep = 499)

        _, _, qw0 = fastwildboot_null(fitted_w, m, Rademacher(),  rep = 499)
        _, _, qu0 = fastwildboot_null(fitted_u, m, Clusters.Rademacher(),  rep = 499)

        [theta_u, V1_u, V2_u, V3_u, V4_u, V5_u, G_u, k_u,
        qu[1], qu[4], qu[2], qu[3],
        qu0[1], qu0[4], qu0[2], qu0[3],
        theta_w, V1_w, V2_w, V3_w, V4_w, V5_w, G_w, k_w,
        qw0[1], qw0[4], qw0[2], qw0[3]]
    catch
        fill(NaN, 28)
    end
end

function estimatemodellight(m::LinearRegressionCluster)
    y = m.y
    X = m.X
    cl = m.cl
    iter = m.iter
    w  = Val{:weighted}
    uw = Val{:unweighted}
    qu = zeros(4)
    qw = zeros(4)
    qu0 = zeros(4)
    qw0 = zeros(4)
    V4_w = 1.0
    V5_w = 1.0
    V4_u = 1.0
    V5_u = 1.
    G_u = 1.0
    k_u = 1.0
    G_w = 1.0
    k_w = 1.0
    try
        fitted_u = fit(GeneralizedLinearModel, X, y, Normal(), IdentityLink())

        theta_u = first(coef(fitted_u))
        V1_u = first(stderr(fitted_u))
        V2_u = first(stderr(fitted_u, HC1()))

        V3_u = faststderr(fitted_u, m, CRHC1(cl), uw)
        #V4_u = faststderr(fitted_u, m, CRHC2(cl), uw)
        #V5_u = faststderr(fitted_u, m, CRHC3(cl), uw)

        fitted_w = fit(GeneralizedLinearModel, X, y, Normal(), IdentityLink(), wts = m.wts)
        theta_w = first(coef(fitted_w))
        V1_w = fastiid(fitted_w, m)
        V2_w = first(stderr(fitted_w, HC1()))

        V3_w = faststderr(fitted_w, m, CRHC1(cl), w)
        #V4_w = faststderr(fitted_w, m, CRHC2(cl), w)
        #V5_w = faststderr(fitted_w, m, CRHC3(cl), w)

        #G_u, k_u = kappastar(fitted_u, m, uw)
        #G_w, k_w = kappastar(fitted_w, m, w)

        # _, _, qw = fastwildboot_nonull(fitted_w, m, Rademacher(), rep = 499)
        # _, _, qu = fastwildboot_nonull(fitted_u, m, Rademacher(), rep = 499)

        # _, _, qw0 = fastwildboot_null(fitted_w, m, Rademacher(),  rep = 499)
        # _, _, qu0 = fastwildboot_null(fitted_u, m, Clusters.Rademacher(),  rep = 499)

        [theta_u, V1_u, V2_u, V3_u, V4_u, V5_u, G_u, k_u,
        qu[1], qu[4], qu[2], qu[3],
        qu0[1], qu0[4], qu0[2], qu0[3],
        theta_w, V1_w, V2_w, V3_w, V4_w, V5_w, G_w, k_w,
        qw0[1], qw0[4], qw0[2], qw0[3]]
    catch
        fill(NaN, 28)
    end
end

##=
## Standard Errors
##=
function fastiid(f::GLM.AbstractGLM, m::LinearRegressionCluster)
    r = f.rr.wrkresid.*m.sqwts
    ichol = inv(cholfact(f.pp))
    sqrt(ichol[1]*mean(abs2.(r)))
end

function faststderr(f::GLM.AbstractGLM, m::LinearRegressionCluster, v::CovarianceMatrices.CRHC, w::T) where T
    B = fastmeat(f, m, v, w)
    A = first(inv(f.pp.chol))
    sqrt(B*A^2)
end

function fastmeat(f::GLM.AbstractGLM, m::LinearRegressionCluster, v::CovarianceMatrices.CRHC,
    w::Type{Val{:unweighted}})
    cl = m.cl
    bstarts = m.bstarts
    ichol  = inv(cholfact(f.pp))::Array{Float64, 2}
    e = copy(CovarianceMatrices.wrkresid(f.rr))
    CovarianceMatrices.adjresid!(v, m.X, e, ichol, bstarts)
    fastclusterize(m.X.*e, bstarts)
end

function fastmeat(f::GLM.AbstractGLM, m::LinearRegressionCluster, v::CovarianceMatrices.CRHC,
    w::Type{Val{:weighted}})
    cl = m.cl
    bstarts = m.bstarts
    ichol  = inv(cholfact(f.pp))::Array{Float64, 2}
    e = copy(CovarianceMatrices.wrkresid(f.rr))
    X = copy(m.X)
    broadcast!(*, X, X, m.sqwts)
    broadcast!(*, e, e, m.sqwts)
    CovarianceMatrices.adjresid!(v, X, e, ichol, bstarts)
    fastclusterize(X.*e, bstarts)
end

function fastmeat(ichol, e, X, m::LinearRegressionCluster, v::CovarianceMatrices.CRHC)
    bstarts = m.bstarts
    cl  = m.cl
    CovarianceMatrices.adjresid!(v, X, e, ichol, bstarts)
    fastclusterize(X.*e, bstarts)
end

function fastclusterize(U, bstarts)
    M = 0.0
    G = length(bstarts)
    for g = 1:G
        s = 0.0
        for i = bstarts[g]
            @inbounds s += U[i, 1]
        end
        M += s^2
    end
    M
end

####################################################################################
##
##
####################################################################################
function fastgamma(f::GLM.AbstractGLM, m::LinearRegressionCluster, ::Type{Val{:weighted}})
    A = first(inv(f.pp.chol))
    X = copy(m.X)
    X .= X.*m.sqwts
    fastgamma(A, X, m.bstarts)
end

function fastgamma(f::GLM.AbstractGLM, m::LinearRegressionCluster, ::Type{Val{:unweighted}})
    A = first(inv(f.pp.chol))
    X = copy(m.X)
    fastgamma(A, X, m.bstarts)
end

function kappastar(f::GLM.AbstractGLM, m::LinearRegressionCluster, w::T) where T
    G = m.opt.G
    Gamma = fastgamma(f, m, w)
    (Gamma, G/(1+Gamma))
end

function fastgamma(XXinv, X, bstarts)
    M = Array{Float64}(length(bstarts))
    kappaclusterize!(M, X, bstarts)
    gamma_g = XXinv^2.*M
    gamma_bar = mean(gamma_g)
    mean((gamma_g-gamma_bar).^2)/gamma_bar.^2
end

function kappaclusterize!(M, U, bstarts)
    ## Assumes that U is ng x 1
    G = length(bstarts)
    for m = 1:G
        s = 0.0
        for i = bstarts[m]
            @inbounds s += U[i, 1]
        end
        M[m] = s^2
    end
end

function fastwildboot_null(f::GLM.AbstractGLM, m::LinearRegressionCluster, WT::WildWeights; rep::Int = 499)
    ## Wildbootstrap with null imposed
    wts  = f.rr.wts
    is_weighted = !all(wts.==1.0)
    sqwts = m.sqwts
    betahat = first(coef(f))
    Hₙ = first(inv(f.pp.chol))
    if is_weighted
        Y = f.rr.y.*sqwts
        X = m.X.*sqwts
        uhat = copy(Y)
    else
        Y = f.rr.y
        X = copy(m.X)
        uhat = Y
    end
    β, σ = Clusters._wildboot_null(Y, X, uhat, Hₙ, WT, rep, m)
    q = quantile(β./σ, [.025, .05, .95,.975])
    (β, σ, q)
end

function fastwildboot_nonull(f::GLM.AbstractGLM, m::LinearRegressionCluster, WT::WildWeights; rep::Int = 999)
    wts  = f.rr.wts
    is_weighted = !all(wts.==1)
    sqwts = m.sqwts
    betahat = first(coef(f))
    Hₙ = first(inv(f.pp.chol))
    uhat = CovarianceMatrices.wrkresidwts(f.rr)
    if is_weighted
        Y = f.rr.y.*sqwts
        X = m.X.*sqwts
    else
        Y = f.rr.y
        X = copy(m.X)
    end
    β, σ = Clusters._wildboot_nonull(Y, X, uhat, betahat, Hₙ, WT, rep, m)
    q = quantile((β-betahat)./σ, [.025, .05, .95,.975])
    (β, σ, q)
end

function _wildboot_null(Y, X, uhat, Hₙ, WT, rep, m::LinearRegressionCluster)
    ustar = similar(uhat)
    bstarts = m.bstarts
    G = length(bstarts)
    ## Containers
    W  = Array{Float32}(G)
    β  = Array{Float64}(rep)
    σ  = Array{Float64}(rep)
    cr = CRHC1(m.cl)
    cX = similar(X)
    for h in 1:rep
        Clusters.wbweights!(WT, W)
        s = 0.0
        for r = 1:G
            @inbounds for i = bstarts[r]
                ustar[i] = uhat[i]*W[r]
                s += X[i]*ustar[i]
            end
        end
        ##
        β[h] = Hₙ*s
        ## Calculate the variance
        @simd for j in eachindex(ustar)
            @inbounds ustar[j] -= X[j]*Hₙ*s
        end
        B = Clusters.fastmeat(Hₙ, ustar, copy!(cX, X), m, cr)
        σ[h] = sqrt(B*Hₙ^2)
    end
    (β, σ)
end

function _wildboot_nonull(Y, X, uhat, betahat, Hₙ, WT, rep, m::LinearRegressionCluster)
    ustar = similar(uhat)
    bstarts = m.bstarts
    G = length(bstarts)
    ## Containers
    W = Array{Float32}(G)
    β = Array{Float64}(rep)
    σ = Array{Float64}(rep)
    cr = CRHC1(m.cl)
    cX = similar(X)
    for h in 1:rep
        Clusters.wbweights!(WT, W)
        s = 0.0
        for r = 1:G
            @inbounds for i = bstarts[r]
                ustar[i] = uhat[i]*W[r]
                s += X[i]*ustar[i]
            end
        end
        ##
        β[h] = betahat .+ Hₙ*s
        ## Calculate the variance
        @simd for j in eachindex(ustar)
            @inbounds ustar[j] -= X[j]*Hₙ*s
        end
        B = Clusters.fastmeat(Hₙ, ustar, copy!(cX, X), m, cr)
        σ[h] = sqrt(B*Hₙ^2)
    end
    (β, σ)
end

function wbweights!(x::Rademacher, W::Vector{Float32})
    G = length(W)
    for m = 1:G
        W[m] = rand() < 0.5 ? 1.0 : -1.0
    end
end

function wbweights!(x::TwoPoint, W::Vector{Float32})
    G = length(W)
    s5 = √5
    p = (s5 + 1.0)/(2*s5)
    for m = 1:G
        s = 0.0
        W[m] = rand() < p ? (1 + s5)/2 : -(s5-1)/2
    end
end

export LinearRegressionCluster, LinearRegressionClusterOpt,
initialize, simulate!, estimatemodel, montecarlo,
estimatemodellight, montecarlolight,
ConfInterval, coverage, average_length, icc_x, icc_u

end


text:
"""
$(TYPEDEF)

Hypograph of perspective function of real symmetric or complex Hermitian
log-determinant cone of dimension `dim` in svec format.

    $(FUNCTIONNAME){T, R}(dim::Int, use_dual::Bool = false)
"""
mutable struct HypoPerLogdetTri{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    d::Int
    is_complex::Bool
    rt2::T

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    ϕ::T
    ζ::T
    mat::Matrix{R}
    fact_W::Cholesky{R}
    Wi::Matrix{R}
    Wi_vec::Vector{T}
    mat2::Matrix{R}
    mat3::Matrix{R}
    mat4::Matrix{R}

    function HypoPerLogdetTri{T, R}(
        dim::Int;
        use_dual::Bool = false,
        ) where {T <: Real, R <: RealOrComplex{T}}
        @assert dim >= 3
        cone = new{T, R}()
        cone.use_dual_barrier = use_dual
        cone.dim = dim
        cone.rt2 = sqrt(T(2))
        cone.is_complex = (R <: Complex)
        cone.d = svec_side(R, dim - 2)
        return cone
    end
end

function setup_extra_data!(
    cone::HypoPerLogdetTri{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    dim = cone.dim
    d = cone.d
    cone.mat = zeros(R, d, d)
    cone.Wi = zeros(R, d, d)
    cone.Wi_vec = zeros(T, dim - 2)
    cone.mat2 = zeros(R, d, d)
    cone.mat3 = zeros(R, d, d)
    cone.mat4 = zeros(R, d, d)
    return cone
end

get_nu(cone::HypoPerLogdetTri) = 2 + cone.d

function set_initial_point!(
    arr::AbstractVector{T},
    cone::HypoPerLogdetTri{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    arr .= 0
    # central point data are the same as for hypoperlog
    (arr[1], arr[2], w) = get_central_ray_hypoperlog(cone.d)
    incr = (cone.is_complex ? 2 : 1)
    k = 3
    @inbounds for i in 1:cone.d
        arr[k] = w
        k += incr * i + 1
    end
    return arr
end

function update_feas(cone::HypoPerLogdetTri{T}) where {T <: Real}
    @assert !cone.feas_updated
    v = cone.point[2]

    if v > eps(T)
        u = cone.point[1]
        @views svec_to_smat!(cone.mat, cone.point[3:end], cone.rt2)
        fact = cone.fact_W = cholesky!(Hermitian(cone.mat, :U), check = false)
        if isposdef(fact)
            cone.ϕ = logdet(fact) - cone.d * log(v)
            cone.ζ = v * cone.ϕ - u
            cone.is_feas = (cone.ζ > eps(T))
        else
            cone.is_feas = false
        end
    else
        cone.is_feas = false
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::HypoPerLogdetTri{T}) where {T <: Real}
    u = cone.dual_point[1]

    if u < -eps(T)
        v = cone.dual_point[2]
        @views svec_to_smat!(cone.mat2, cone.dual_point[3:end], cone.rt2)
        fact = cholesky!(Hermitian(cone.mat2, :U), check = false)
        if isposdef(fact)
            return (v - u * (logdet(fact) + cone.d * (1 - log(-u))) > eps(T))
        end
    end

    return false
end

function update_grad(cone::HypoPerLogdetTri)
    @assert cone.is_feas
    v = cone.point[2]
    g = cone.grad
    ζ = cone.ζ
    ζi = inv(ζ)

    g[1] = ζi
    g[2] = -inv(v) - (cone.ϕ - cone.d) / ζ
    inv_fact!(cone.Wi, cone.fact_W)
    smat_to_svec!(cone.Wi_vec, cone.Wi, cone.rt2)
    # ∇ϕ = cone.Wi_vec * v
    vζi1 = -1 - v / ζ
    @. g[3:end] = vζi1 * cone.Wi_vec

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::HypoPerLogdetTri)
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    v = cone.point[2]
    H = cone.hess.data
    d = cone.d
    ζ = cone.ζ
    Wi_vec = cone.Wi_vec
    σζi = (cone.ϕ - d) / ζ
    vζi = v / ζ

    # u, v
    H[1, 1] = ζ^-2
    H[1, 2] = -σζi / ζ
    H[2, 2] = v^-2 + abs2(σζi) + d / (ζ * v)

    # u, v, w
    vζi2 = -vζi / ζ
    c1 = ((cone.ϕ - d) * vζi - 1) / ζ
    @. H[1, 3:end] = vζi2 * Wi_vec
    @. H[2, 3:end] = c1 * Wi_vec

    # w, w
    copytri!(cone.Wi, 'U', true)
    @views symm_kron!(H[3:end, 3:end], cone.Wi, cone.rt2)

    @inbounds for j in eachindex(Wi_vec)
        j2 = 2 + j
        c2 = vζi * Wi_vec[j]
        for i in 1:j
            i2 = 2 + i
            H[i2, j2] += vζi * (H[i2, j2] + c2 * Wi_vec[i])
        end
    end

    cone.hess_updated = true
    return cone.hess
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::HypoPerLogdetTri,
    )
    @assert cone.grad_updated
    v = cone.point[2]
    FU = cone.fact_W.U
    d = cone.d
    ζ = cone.ζ
    σ = cone.ϕ - d
    vζi1 = (v + ζ) / ζ
    w_aux = cone.mat3

    @inbounds for j in 1:size(arr, 2)
        p = arr[1, j]
        q = arr[2, j]
        @views svec_to_smat!(w_aux, arr[3:end, j], cone.rt2)
        copytri!(w_aux, 'U', true)
        rdiv!(w_aux, FU)
        ldiv!(FU', w_aux)

        qζi = q / ζ
        c1 = tr(Hermitian(w_aux, :U)) / ζ
        # ∇ϕ[r] = v * c1
        c2 = (v * c1 - p / ζ + σ * qζi) / ζ
        c3 = c2 * v - qζi

        prod[1, j] = -c2
        prod[2, j] = c2 * σ - c1 + (qζi * d + q / v) / v

        lmul!(vζi1, w_aux)
        for i in diagind(w_aux)
            w_aux[i] += c3
        end
        rdiv!(w_aux, FU')
        ldiv!(FU, w_aux)
        @views smat_to_svec!(prod[3:end, j], w_aux, cone.rt2)
    end

    return prod
end

function update_inv_hess(cone::HypoPerLogdetTri)
    @assert cone.grad_updated
    @assert !cone.inv_hess_updated
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    v = cone.point[2]
    @views w = cone.point[3:end]
    svec_to_smat!(cone.mat2, w, cone.rt2)
    W = Hermitian(cone.mat2, :U)
    Hi = cone.inv_hess.data
    d = cone.d
    ζ = cone.ζ
    ϕ = cone.ϕ
    ζv = ζ + v
    ζζvi = ζ / ζv
    c1 = v / (ζv + d * v) * v
    c2 = c1 / ζv

    # u, v
    Hi12 = Hi[1, 2] = c1 * (ζv * ϕ - d * ζ)
    Hi[1, 1] = ζζvi * d * v^2 + ζ^2 + (ϕ - ζζvi * d) * Hi12
    Hi[2, 2] = c1 * ζv

    # u, v, w
    c3 = (v * ζ + Hi12) / ζv
    @. Hi[1, 3:end] = c3 * w
    @. Hi[2, 3:end] = c1 * w

    # w, w
    @views Hiww = Hi[3:end, 3:end]
    symm_kron!(Hiww, W, cone.rt2)
    mul!(Hiww, w, w', c2, ζζvi)

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::HypoPerLogdetTri,
    )
    @assert cone.grad_updated
    v = cone.point[2]
    @views w = cone.point[3:end]
    svec_to_smat!(cone.mat4, w, cone.rt2)
    W = Hermitian(cone.mat4, :U)
    d = cone.d
    ζ = cone.ζ
    ϕ = cone.ϕ
    ζv = ζ + v
    ζζvi = ζ / ζv
    c1 = ζv / v + d
    w_aux = cone.mat2
    w_aux2 = cone.mat3

    @inbounds for j in 1:size(arr, 2)
        p = arr[1, j]
        q = arr[2, j]
        @views r = arr[3:end, j]
        svec_to_smat!(w_aux, r, cone.rt2)
        copytri!(w_aux, 'U', true)

        trrw = dot(w, r)
        c2 = v * (ζv * (ϕ * p + q) - d * ζ * p + trrw) / c1
        c3 = (c2 + ζ * v * p) / ζv

        prod[1, j] = ζ * ((v * (d * p * v + trrw) - d * c2) / ζv + ζ * p) + ϕ * c2
        prod[2, j] = c2

        mul!(w_aux2, w_aux, W)
        mul!(w_aux, W, w_aux2)
        @views prod_w = prod[3:end, j]
        smat_to_svec!(prod_w, w_aux, cone.rt2)
        axpby!(c3, w, ζζvi, prod_w)
    end

    return prod
end

function dder3(cone::HypoPerLogdetTri, dir::AbstractVector)
    @assert cone.grad_updated
    v = cone.point[2]
    dder3 = cone.dder3
    p = dir[1]
    q = dir[2]
    d = cone.d
    ζ = cone.ζ
    FU = cone.fact_W.U
    σ = cone.ϕ - d
    viq = q / v
    vζi1 = (v + ζ) / ζ
    rwi = cone.mat2
    w_aux = cone.mat3
    w_aux2 = cone.mat4

    @views svec_to_smat!(rwi, dir[3:end], cone.rt2)
    copytri!(rwi, 'U', true)
    rdiv!(rwi, FU)
    ldiv!(FU', rwi)

    tr1 = tr(Hermitian(rwi, :U))
    tr2 = sum(abs2, rwi)

    χ = (-p + σ * q + tr1 * v) / ζ
    c1 = (viq * (2 * tr1 - viq * d) - tr2) / (2 * ζ)
    c2 = (abs2(χ) - v * c1) / ζ
    c3 = -q * χ / ζ
    c4 = (χ * v - q) / ζ
    c5 = c3 + c2 * v

    dder3[1] = -c2
    dder3[2] = c2 * σ + (abs2(viq) - d * c3 - c4 * tr1) / v - tr2 / ζ - c1

    copyto!(w_aux2, I)
    axpby!(vζi1, rwi, c4, w_aux2)
    mul!(w_aux, Hermitian(rwi, :U), w_aux2)
    @inbounds for i in diagind(w_aux)
        w_aux[i] += c5
    end
    rdiv!(w_aux, FU')
    ldiv!(FU, w_aux)
    @views smat_to_svec!(dder3[3:end], w_aux, cone.rt2)

    return dder3
end


text:
using FromFile
@from "Core.jl" import CONST_TYPE, MAX_DEGREE, Node, Options, Dataset
@from "EquationUtils.jl" import stringTree
@from "PopMember.jl" import PopMember
@from "LossFunctions.jl" import EvalLoss
using Printf: @sprintf

""" List of the best members seen all time in `.members` """
mutable struct HallOfFame
    members::Array{PopMember, 1}
    exists::Array{Bool, 1} #Whether it has been set

    # Arranged by complexity - store one at each.
end

"""
    HallOfFame(options::Options)

Create empty HallOfFame. The HallOfFame stores a list
of `PopMember` objects in `.members`, which is enumerated
by size (i.e., `.members[1]` is the constant solution).
`.exists` is used to determine whether the particular member
has been instantiated or not.
"""
function HallOfFame(options::Options)
    actualMaxsize = options.maxsize + MAX_DEGREE
    HallOfFame([PopMember(Node(convert(CONST_TYPE, 1)), 1f9) for i=1:actualMaxsize], [false for i=1:actualMaxsize])
end

"""
    calculateParetoFrontier(dataset::Dataset{T}, hallOfFame::HallOfFame,
                            options::Options) where {T<:Real}
"""
function calculateParetoFrontier(dataset::Dataset{T},
                                 hallOfFame::HallOfFame,
                                 options::Options)::Array{PopMember, 1} where {T<:Real}
    # Dominating pareto curve - must be better than all simpler equations
    dominating = PopMember[]
    actualMaxsize = options.maxsize + MAX_DEGREE
    for size=1:actualMaxsize
        if hallOfFame.exists[size]
            member = hallOfFame.members[size]
            curMSE = EvalLoss(member.tree, dataset, options)
            member.score = curMSE
            betterThanAllSmaller = all([
                (!(hallOfFame.exists[i])
                 || curMSE < EvalLoss(hallOfFame.members[i].tree, dataset, options))
                for i=1:(size-1)
            ])
            if betterThanAllSmaller
                push!(dominating, member)
            end
        end
    end
    return dominating
end

"""
    calculateParetoFrontier(X::AbstractMatrix{T}, y::AbstractVector{T},
                            hallOfFame::HallOfFame, options::Options;
                            weights=nothing, varMap=nothing) where {T<:Real}

Compute the dominating Pareto frontier for a given hallOfFame. This
is the list of equations where each equation has a better loss than all
simpler equations.
"""
function calculateParetoFrontier(X::AbstractMatrix{T},
                                 y::AbstractVector{T},
                                 hallOfFame::HallOfFame,
                                 options::Options;
                                 weights=nothing,
                                 varMap=nothing) where {T<:Real}
    calculateParetoFrontier(Dataset(X, y, weights=weights, varMap=varMap), hallOfFame, options)
end

function string_dominating_pareto_curve(hallOfFame, baselineMSE,
                                        dataset, options,
                                        avgy)
    output = ""
    curMSE = baselineMSE
    lastMSE = curMSE
    lastComplexity = 0
    output *= "Hall of Fame:\n"
    output *= "-----------------------------------------\n"
    output *= @sprintf("%-10s  %-8s   %-8s  %-8s\n", "Complexity", "Loss", "Score", "Equation")

    #TODO - call pareto function!
    actualMaxsize = options.maxsize + MAX_DEGREE
    for size=1:actualMaxsize
        if hallOfFame.exists[size]
            member = hallOfFame.members[size]
            curMSE = EvalLoss(member.tree, dataset, options)
            betterThanAllSmaller = all([
                    (
                         !(hallOfFame.exists[i])
                         || curMSE < EvalLoss(hallOfFame.members[i].tree, dataset, options)
                    ) for i=1:(size-1)
               ])
            if betterThanAllSmaller
                delta_c = size - lastComplexity
                delta_l_mse = log(curMSE/lastMSE)
                score = convert(Float32, -delta_l_mse/delta_c)
                output *= @sprintf("%-10d  %-8.3e  %-8.3e  %-s\n" , size, curMSE, score, stringTree(member.tree, options, varMap=dataset.varMap))
                lastMSE = curMSE
                lastComplexity = size
            end
        end
    end
    output *= "\n"
    return output
end


text:
include("common.jl")
import ProfileView
import Profile

sim = init_sim("Bomex")
tendencies = copy(sim.state.prog)

Profile.@profile update_n(sim, tendencies, 100)
Profile.print()
# Profile.print(; format = :flat, sortedby = :count)
ProfileView.@profview update_n(sim, tendencies, 1000) # compile first
ProfileView.@profview update_n(sim, tendencies, 1000)
Profile.clear_malloc_data()


text:
#=
const so_suffix = string(sys_suffix, ".so")

function call_make_one(c::String,
                    path_bin::String, path_code::String,
                    so_suffix::String=so_suffix,
            )
    so = string(c, so_suffix)
    so_to = string(path_bin, "/", so)
    if isfile(so_to) return end
    @show so_to
    path_source = string(path_code, "/", c)
    run(`make -C $path_source all`)
    @show string(path_source, "/", so)
    cp(string(path_source, "/", so), so_to, )
    ncf = readlines(string(path_source, "/cpfiles"))
    for i in ncf
        isempty(i) && continue
        cp(string(path_source, "/", i), string(path_bin, "/", i), force=true)
    end

end

function call_make(path_bin::String, path_code::String)
    sod = ["apex", "chapman", "hwm14", "nrlmsise00", "spa"]
    for i in sod
        call_make_one(i, path_bin, path_code)
    end

end

# call_make(bin_dir, string(bin_dir, "/code"),) =#

include(string(so_dir, "/apex/", wrapperjlfile))
include(string(so_dir, "/chapman/", wrapperjlfile))
include(string(so_dir, "/hwm14/", wrapperjlfile))
include(string(so_dir, "/iri2016/", wrapperjlfile))
include(string(so_dir, "/nrlmsise00/", wrapperjlfile))
include(string(so_dir, "/spa/", wrapperjlfile))


text:
export ALL_solvers_stopping

# Valid combinations

ALL_solvers_stopping = Function[]

include("ARCSpectral.jl")
push!(ALL_solvers_stopping, ARCSpectral)

include("ARCSpectral_abs.jl")
push!(ALL_solvers_stopping, ARCSpectral_abs)

include("TRSpectral.jl")
push!(ALL_solvers_stopping, TRSpectral)

include("TRSpectral_abs.jl")
push!(ALL_solvers_stopping, TRSpectral_abs)

include("ARCLDLt.jl")
push!(ALL_solvers_stopping, ARCLDLt)

include("ARCLDLt_abs.jl")
push!(ALL_solvers_stopping, ARCLDLt_abs)

include("TRLDLt.jl")
push!(ALL_solvers_stopping, TRLDLt)

include("TRLDLt_HO.jl")
push!(ALL_solvers_stopping, TRLDLt_HO)

include("TRLDLt_BBK.jl")
push!(ALL_solvers_stopping, TRLDLt_BK)

include("TRLDLt_HO_Sham.jl")
push!(ALL_solvers_stopping, TRLDLt_HO_Sham)

include("TRLDLt_HO_Sham_lambda.jl")
push!(ALL_solvers_stopping, TRLDLt_HO_Sham_λ)

include("TRLDLt_HO_Sham_BFGS.jl")
push!(ALL_solvers_stopping, TRLDLt_HO_Sham_BFGS)

# include("TRLDLt_HO_MP.jl")                # Still a work in progress
# push!(ALL_solvers_stopping, TRLDLt_HO_MP)

include("TRLDLt_abs.jl")
push!(ALL_solvers_stopping, TRLDLt_abs)

include("ARCqKOp.jl")
push!(ALL_solvers_stopping, ARCqKOp)

include("ARCqKOp1.jl")
push!(ALL_solvers_stopping, ARCqKOp1)

include("ARCqKOp1b.jl")
push!(ALL_solvers_stopping, ARCqKOp1b)

include("ARCqKOp2.jl")
push!(ALL_solvers_stopping, ARCqKOp2)

include("ARCqKOp0_5.jl")
push!(ALL_solvers_stopping, ARCqKOp05)

include("ARCqKSparse.jl")
push!(ALL_solvers_stopping, ARCqKsparse)

include("ARCqKdense.jl")
push!(ALL_solvers_stopping, ARCqKdense)

include("TRKOp.jl")
push!(ALL_solvers_stopping, TRKOp)

include("TRKsparse.jl")
push!(ALL_solvers_stopping, TRKsparse)

include("TRKdense.jl")
push!(ALL_solvers_stopping, TRKdense)

include("ST_TROp.jl")
push!(ALL_solvers_stopping, ST_TROp)

include("ST_TROp_Sham.jl")
push!(ALL_solvers_stopping, ST_TROp_Sham)

include("ST_TROp_Sham_BFGS.jl")
push!(ALL_solvers_stopping, ST_TROp_Sham_BFGS)

include("ST_TRsparse.jl")
push!(ALL_solvers_stopping, ST_TRsparse)

include("ST_TRdense.jl")
push!(ALL_solvers_stopping, ST_TRdense)

include("ST_ARCOp.jl")
push!(ALL_solvers_stopping, ST_ARCOp)

include("ST_ARCsparse.jl")
push!(ALL_solvers_stopping, ST_ARCsparse)

include("ST_ARCdense.jl")
push!(ALL_solvers_stopping, ST_ARCdense)

include("ARCMA97.jl")
push!(ALL_solvers_stopping, ARCMA97)

include("ARCMA97_abs.jl")
push!(ALL_solvers_stopping, ARCMA97_abs)

include("TRMA97.jl")
push!(ALL_solvers_stopping, TRMA97)

include("TRMA97_abs.jl")
push!(ALL_solvers_stopping, TRMA97_abs)

include("ARCMA57.jl")
push!(ALL_solvers_stopping, ARCMA57)

include("ARCMA57_abs.jl")
push!(ALL_solvers_stopping, ARCMA57_abs)

include("TRMA57.jl")
push!(ALL_solvers_stopping, TRMA57)

include("TRMA57_Sham.jl")
push!(ALL_solvers_stopping, TRMA57_Sham)

include("TRMA57-2.jl")
push!(ALL_solvers_stopping, TRMA57_2)

include("TRMA57_Sham_2.jl")
push!(ALL_solvers_stopping, TRMA57_Sham_2)

include("TRMA57_Sham_BFGS.jl")
push!(ALL_solvers_stopping, TRMA57_Sham_BFGS)

include("TRMA57_abs.jl")
push!(ALL_solvers_stopping, TRMA57_abs)


text:
module NBTest_Water_Analysis
using NBTesting
is_testing = true

NBTesting.clear_test_times!()

# -----------------------------------
println("# Top title"); flush(STDOUT)

# Cell #1
using NBTesting
using Base.Test

using PyPlot

# This defines
#     plot(args...; kwargs...) = nothing
# but only during testing.
@testing_noop plot

# -----------------------------------
println("### Basic Arithmetic"); flush(STDOUT)

# Cell #2
x = 10 / 2
@test x+x == 10
@test 6*10 == 60

# Cell #3
N = is_testing ? 5 : 100

srand(1)
M1 = ones(N, N)
M2 = rand(N, N);

if is_testing
    @test 30 < sum(M1 * M2) < 100
end

sum(M1 * M2)

# -----------------------------------
println("#### Some plot"); flush(STDOUT)

# Cell #4
plot([1,2,1]);

end  # module


text:
using OutlierDetectionNeighbors
using OutlierDetectionTest
using Distances: Cityblock

# Test the metadata of all exported detectors
test_meta.(eval.(OutlierDetectionNeighbors.MODELS))

data = TestData()
run_test(detector) = test_detector(detector, data)

# ABOD
run_test(ABODDetector())
run_test(ABODDetector(static=false))
run_test(ABODDetector(parallel=true, enhanced=true))

# COF
run_test(COFDetector())
run_test(COFDetector(metric=Cityblock()))
run_test(COFDetector(static=false))
run_test(COFDetector(parallel=true))

# DNN
run_test(DNNDetector(d=1))
run_test(DNNDetector(d=1, static=false))
run_test(DNNDetector(d=1, parallel=true))

# KNN
run_test(KNNDetector(reduction=:maximum))
run_test(KNNDetector(static=false))
run_test(KNNDetector(parallel=true, reduction=:mean, algorithm=:balltree))
run_test(KNNDetector(parallel=true, reduction=:median))

# LOF
run_test(LOFDetector())
run_test(LOFDetector(static=false))
run_test(LOFDetector(parallel=true))


text:
using BandedMatrices, LinearAlgebra, LazyArrays, Test
import LazyArrays: MemoryLayout, MulAdd, DenseColumnMajor, ConjLayout
import Base: BroadcastStyle
import BandedMatrices: BandedStyle, BandedRows

@testset "broadcasting" begin
    @testset "general" begin
        n = 1000
        A = brand(n,n,1,1)
        B = Matrix{Float64}(undef, n,n)
        B .= exp.(A)
        @test B == exp.(Matrix(A)) == exp.(A)

        @test exp.(A) isa BandedMatrix
        @test bandwidths(exp.(A)) == (n-1,n-1)

        C = similar(A)
        @test_throws BandError C .= exp.(A)

        @test identity.(A) isa BandedMatrix
        @test bandwidths(identity.(A)) == bandwidths(A)

        @test (z -> exp(z)-1).(A) isa BandedMatrix
        @test bandwidths((z -> exp(z)-1).(A)) == bandwidths(A)

        @test A .+ 1 isa BandedMatrix
        @test (A .+ 1) == Matrix(A) .+ 1

        @test A ./ 1 isa BandedMatrix
        @test bandwidths(A ./ 1) == bandwidths(A)

        @test 1 .+ A isa BandedMatrix
        @test (1 .+ A) == 1 .+ Matrix(A)

        @test 1 .\ A isa BandedMatrix
        @test bandwidths(1 .\ A) == bandwidths(A)

        A = brand(10,1,1,1)
        @test A[:,1] isa Vector
        @test norm(A .- A[:,1]) == 0
        @test A ≈ A[:,1]
    end

    @testset "identity" begin
        n = 100
        A = brand(n,n,2,2)
        A.data[1,1] = NaN
        B = brand(n,n,2,2)
        B.data[1,1] = NaN
        A .= B
        @test A == B
        B = brand(n,n,1,1)
        B.data[1,1] = NaN
        A .= B
        B = brand(n,n,3,3)
        B[band(3)] .= 0
        B[band(-3)] .= 0
        B.data[1,1] = NaN
        B.data[end,end] = NaN
        A .= B
        @test A == B

        B = brand(n,n,0,3)
        B[band(3)] .= 0
        B.data[1,1] = NaN
        A .= B
        @test A == B

        B = brand(n,n,3,0)
        B[band(-3)] .= 0
        B.data[end,end] = NaN
        A .= B
        @test A == B

        B = brand(n,n,-1,1)
        A .= B
        @test A == B

        B = brand(n,n,1,-1)
        A .= B
        @test A == B
    end

    @testset "lmul!/rmul!" begin
        n = 1000
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B .= (-).(A)
        @test -A isa BandedMatrix
        @test (-).(A) isa BandedMatrix
        @test bandwidths(A) == bandwidths(-A) == bandwidths((-).(A))
        @test B == -A == (-).(A)
        @test A-I isa BandedMatrix
        @test I-A isa BandedMatrix
        @test bandwidths(A) == bandwidths(A-I) == bandwidths(I-A)

        n = 1000
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B .= 2.0.*A

        @test B ==  2A == 2.0.*A
        @test 2A isa BandedMatrix
        @test 2.0.*A isa BandedMatrix
        @test bandwidths(2A) == bandwidths(2.0.*A) == bandwidths(A)

        A .= 2.0.*A
        @test A == B

        n = 1000
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B .= A.*2.0

        @test B ==  A*2 == A.*2.0
        @test A*2 isa BandedMatrix
        @test A .* 2.0 isa BandedMatrix
        @test bandwidths(A*2) == bandwidths(A.*2.0) == bandwidths(A)
        A .= A.*2.0
        @test A == B

        n = 1000
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B .= A ./ 2.0

        @test B == A/2 == A ./ 2.0
        @test A/2 isa BandedMatrix
        @test A ./ 2.0 isa BandedMatrix
        @test bandwidths(A/2) == bandwidths(A ./ 2.0) == bandwidths(A)

        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B .= 2.0 .\ A

        @test B == A/2 == A ./ 2.0
        @test 2\A isa BandedMatrix
        @test 2.0 .\ A isa BandedMatrix
        @test bandwidths(2\A) == bandwidths(2.0 .\ A) == bandwidths(A)

        A = brand(5,5,1,1)
        A.data .= NaN
        lmul!(0.0,A)
        @test isnan(norm(A)) == isnan(norm(lmul!(0.0,[NaN])))

        lmul!(false,A)
        @test norm(A) == 0.0

        A = brand(5,5,1,1)
        A.data .= NaN
        rmul!(A,0.0)
        @test isnan(norm(A)) == isnan(norm(rmul!([NaN],0.0)))

        rmul!(A,false)
        @test norm(A) == 0.0

        n = 100
        A = brand(n,n,2,2)
        B = brand(n,n,1,1)
        A[band(2)] .= A[band(-2)] .= 0
        B .= A ./ 2.0
        @test B == A / 2.0 == Matrix(A)/2.0

        B .= 2.0 .\ A
        @test B == 2.0 \ A == 2.0 \ Matrix(A)

        n = 100
        A = brand(n,n,2,2)
        B = brand(n,n,1,3)
        A[band(-2)] .= 0
        B .= A ./ 2.0
        @test B == A / 2.0 == Matrix(A)/2.0

        B .= 2.0 .\ A
        @test B == 2.0 \ A == 2.0 \ Matrix(A)

        A = brand(n,n,2,2)
        B = brand(n,n,3,1)
        A[band(2)] .= 0
        B .= A ./ 2.0
        @test B == A / 2.0 == Matrix(A)/2.0

        B .= 2.0 .\ A
        @test B == 2.0 \ A == 2.0 \ Matrix(A)
    end

    @testset "axpy!" begin
        n = 100
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        C = brand(n,n,3,3)
        C .= A .+ B
        @test C == A + B == A .+ B  == Matrix(A) + Matrix(B)
        @test A + B isa BandedMatrix
        @test A .+ B isa BandedMatrix
        @test bandwidths(A+B) == bandwidths(A.+B) == (2,2)
        B .= A .+ B
        @test B == C

        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        C = brand(n,n,3,3)
        C .= A .* B
        @test C == A .* B  == Matrix(A) .* Matrix(B)
        @test A .* B isa BandedMatrix
        @test bandwidths(A.*B) == (2,2)
        @time B .= A .* B
        @test B == C

        n = 100
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B[band(2)] .= B[band(-2)] .= 0
        C = brand(n,n,1,1)
        C .= A .+ B
        @test C == A + B == A .+ B == Matrix(A) + Matrix(B)

        n = 100
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B[band(-2)] .= 0
        C = brand(n,n,1,2)
        C .= A .+ B
        @test C == A + B == A .+ B == Matrix(A) + Matrix(B)

        n = 100
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B[band(2)] .= 0
        C = brand(n,n,2,1)
        C .= A .+ B
        @test C == A + B == A .+ B == Matrix(A) + Matrix(B)

        n = 100
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        C = brand(n,n,3,3)

        C .= 2.0 .* A .+ B
        @test C == 2A+B == 2.0.*A .+ B
        @test all(BLAS.axpy!(2.0, A, copy(B)) .=== C)

        @test 2A + B isa BandedMatrix
        @test 2.0.*A .+ B isa BandedMatrix
        @test bandwidths(2A+B) == bandwidths(2.0.*A .+ B) == (2,2)
        B .= 2.0 .* A .+ B
        @test B == C
    end

    @testset "gbmv!" begin
        n = 100
        x = randn(n)
        A = brand(n,n,1,1)
        y = similar(x)
        y .= Mul(A,x)
        @test Matrix(A)*x ≈ y
        @test all(BLAS.gbmv!('N', n, 1, 1, 1.0, A.data, x, 0.0, copy(y)) .===
                    Broadcast.materialize!(MulAdd(1.0,A,x,0.0,similar(y))) .=== y)
        z = similar(y)
        z .= 2.0.*Mul(A,x) .+ 3.0.*y
        @test 2Matrix(A)*x + 3y ≈ z
        @test all(BLAS.gbmv!('N', n, 1, 1, 2.0, A.data, x, 3.0, copy(y)) .=== z)
        @test MemoryLayout(A') == BandedRows(DenseColumnMajor())
        y .= Mul(A',x)
        @test Matrix(A')*x ≈ Matrix(A)'*x ≈ y
        @test all(BLAS.gbmv!('T', n, 1, 1, 1.0, A.data, x, 0.0, copy(y)) .=== y)
        z = similar(y)
        z .= 2.0.*Mul(A',x) .+ 3.0.*y
        @test 2Matrix(A')*x + 3y ≈ z
        @test all(BLAS.gbmv!('T', n, 1, 1, 2.0, A.data, x, 3.0, copy(y)) .=== z)
        @test MemoryLayout(transpose(A)) == BandedRows(DenseColumnMajor())
        y .= Mul(transpose(A),x)
        @test Matrix(A')*x ≈ Matrix(A)'*x ≈ y
        @test all(BLAS.gbmv!('T', n, 1, 1, 1.0, A.data, x, 0.0, copy(y)) .=== y)

        n = 100
        x = randn(ComplexF64,n)
        A = brand(ComplexF64,n,n,1,1)
        y = similar(x)
        y .= Mul(A,x)
        @test Matrix(A)*x ≈ y
        @test all(BLAS.gbmv!('N', n, 1, 1, 1.0+0.0im, A.data, x, 0.0+0.0im, copy(y)) .=== y)
        @test LazyArrays.MemoryLayout(A') == ConjLayout(BandedRows(DenseColumnMajor()))
        z = similar(y)
        z .= (2.0+0.0im).*Mul(A,x) .+ (3.0+0.0im).*y
        @test 2Matrix(A)*x + 3y ≈ z
        @test all(BLAS.gbmv!('N', n, 1, 1, 2.0+0.0im, A.data, x, 3.0+0.0im, copy(y)) .=== z)
        y .= Mul(A',x)
        @test Matrix(A')*x ≈ Matrix(A)'*x ≈ y
        @test all(BLAS.gbmv!('C', n, 1, 1, 1.0+0.0im, A.data, x, 0.0+0.0im, copy(y)) .=== y)
        @test MemoryLayout(transpose(A)) == BandedRows(DenseColumnMajor())
        y .= Mul(transpose(A),x)
        @test Matrix(transpose(A))*x ≈ transpose(Matrix(A))*x ≈ y
        @test all(BLAS.gbmv!('T', n, 1, 1, 1.0+0.0im, A.data, x, 0.0+0.0im, copy(y)) .=== y)
    end

    @testset "non-blas matvec" begin
        n = 10
        A = brand(Int,n,n,1,1)
        x = rand(Int,n)
        y = similar(x)
        y .= Mul(A,x)
        @test y == Matrix(A)*x
        y .= Mul(A',x)
        @test y == Matrix(A')*x
        y .= Mul(transpose(A),x)
        @test y == Matrix(transpose(A))*x

        A = brand(Int,n,n,1,1) + im*brand(Int,n,n,1,1)
        x = rand(eltype(A),n)
        y = similar(x)
        y .= Mul(A,x)
        @test y == Matrix(A)*x
        y .= Mul(A',x)
        @test y == Matrix(A')*x
        y .= Mul(transpose(A),x)
        @test y == Matrix(transpose(A))*x
    end

    @testset "gbmm!" begin
        n = 10
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        C = brand(n,n,3,3)
        C .= Mul(A,B)
        @test Matrix(C) ≈ Matrix(A)*Matrix(B)
        C .= Mul(A',B)
        @test Matrix(C) ≈ Matrix(A')*Matrix(B)
        C .= Mul(A,B')
        @test Matrix(C) ≈ Matrix(A)*Matrix(B')
        C = brand(n,n,4,4)
        C .= Mul(A,B)
        @test Matrix(C) ≈ Matrix(A)*Matrix(B)
        B = randn(n,n)
        C = similar(B)

        C .= Mul(A,B)
        @test C ≈ Matrix(A)*Matrix(B)
        C .= Mul(B,A)
        @test C ≈ Matrix(B)*Matrix(A)
    end

    @testset "Subarray" begin
        n = 10
        A = brand(n,n,1,1)
        B = view(brand(2n,2n,2,0),Base.OneTo(n),1:n)
        @test BroadcastStyle(typeof(B)) == BandedStyle()

        @test A .+ B isa BandedMatrix
        @test A + B isa BandedMatrix
        @test bandwidths(A .+ B) == (2,1)
        @test A .+ B == Matrix(A) + Matrix(B)
    end
end


text:
@testset "Conv-lx" begin
    cd(td)
    # Exception instead of ArgumentError as may fail with system error
    @test_throws Exception F.check_input_rpath("aldjfk")
end

@testset "Conv-html" begin
    @test_throws F.HTMLFunctionError F.convert_html("{{fill bb cc}}")
    @test_throws F.HTMLFunctionError F.convert_html("{{insert bb cc}}")
    @test_throws F.HTMLFunctionError F.convert_html("{{href aa}}")
    @test (@test_logs (:warn, "Unknown dictionary name aa in {{href ...}}. Ignoring") F.convert_html("{{href aa bb}}")) == "??"
    @test_throws F.HTMLBlockError F.convert_html("{{if asdf}}{{end}}")
    @test_throws F.HTMLBlockError F.convert_html("{{if asdf}}")
    @test_throws F.HTMLBlockError F.convert_html("{{isdef asdf}}")
    @test_throws F.HTMLBlockError F.convert_html("{{ispage asdf}}")
end

@testset "Conv-md" begin
    s = """
        @def blah
        """
    @test (@test_logs (:warn, "Found delimiters for an @def environment but it didn't have the right @def var = ... format. Verify (ignoring for now).") (s |> fd2html_td)) == ""

    s = """
        Blah
        [^1]: hello
        """ |> fd2html_td
    @test isapproxstr(s, "
Blah
")
end

@testset "Franklin" begin
    cd(td); mkpath("foo"); cd("foo"); write("config.md","")
    @test_throws ArgumentError serve(single=true)
    cd(td)
end

@testset "RSS" begin
    F.set_var!(F.GLOBAL_VARS, "website_descr", "")
    F.RSS_DICT["hello"] = F.RSSItem("","","","","","","",Date(1))
    @test (@test_logs (:warn, """
              I found RSS items but the RSS feed is not properly described:
              at least one of the following variables has not been defined in
              your config.md: `website_title`, `website_descr`, `website_url`.
              The feed will not be (re)generated.""") F.rss_generator()) === nothing
end

@testset "parser-lx" begin
    s = raw"""
        \newcommand{hello}{hello}
        """
    @test_throws F.LxDefError (s |> fd2html)
    s = raw"""
        \foo
        """
    @test_throws F.LxComError (s |> fd2html)
    s = raw"""
        \newcommand{\foo}[2]{hello #1 #2}
        \foo{a} {}
        """
    @test_throws F.LxComError (s |> fd2html)
end

@testset "ocblocks" begin
    s = raw"""
        @@foo
        """
    @test_throws F.OCBlockError (s |> fd2html)
end

@testset "tofrom" begin
    s = "jμΛΙα"
    @test F.from(s) == 1
    @test F.to(s) == lastindex(s)
end

@testset "{{toc}}" begin
    s = """
        ~~~
        {{toc 1 2}}
        ~~~
        # Hello
        ## Goodbye
        """ |> fd2html_td
    @test isapproxstr(s, raw"""

Hello

Goodbye

Hello

Goodbye
""")
end

@testset "franklin" begin
    fs1()
    gotd()
    write(joinpath(td, "config.md"), "foo")
    @test_throws ArgumentError serve(single=true)
end


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 781a3590-d1ed-11eb-17ae-fdf1ffd5bdb0
begin
	using Plots
	using PlutoUI
end

# ╔═╡ f2bc9edf-548c-4958-82f4-e65d6fe99f3d
begin
	using Futbolista: Plotter, DataLoaders
end

# ╔═╡ af773761-1267-4dcb-9e60-a5358f6bb319
begin
	const pl = Plotter.StatsbombPitch
	pl.plot_pitch(line_color=:silver, line_width=1)
end

# ╔═╡ e61162bf-c488-4c11-a446-11f65e3b8512
begin
	const sb = DataLoaders.Statsbomb
	ldr = sb.SBLoader("local", raw"C:/repository/open-data/data")
	events_df = sb.load_events(ldr, 22912)

	players = unique(events_df.player)
	players[2]
end

# ╔═╡ f934ef47-7558-4960-98f4-c264772e23b2
@bind player_index Slider(2:29)

# ╔═╡ 82e4a033-347b-423e-a8a3-6329f50023c0
player_name = players[player_index]["name"]

# ╔═╡ f4e11cf7-b68e-4881-b328-2640ef562240
begin
	pass_df = filter(row -> row.type["name"] == "Pass", events_df)
	player_df = filter(row -> row.player["name"] == player_name, pass_df);

	player_df.location = eval.(Meta.parse.(player_df.location))
    player_df[!, "x"] = map(ar -> ar[1], player_df."location")
    player_df[!, "y"] = map(ar -> ar[2], player_df."location")

	player_df[!, "end_x"] = map(d -> d["end_location"][1], player_df."pass")
	player_df[!, "end_y"] = map(d -> d["end_location"][2], player_df."pass")

	p = pl.plot_pitch(;line_color=:silver, line_width=1)
	for row in eachrow(player_df[:, ["x", "y", "end_x", "end_y"]])
		plot!(p, [row.x, row.end_x], [row.y,row.end_y], color=:white,                          linewidth=1.1, label="")
	end
	scatter!(p, player_df.x, player_df.y)
	p
end

# ╔═╡ Cell order:
# ╠═781a3590-d1ed-11eb-17ae-fdf1ffd5bdb0
# ╠═f2bc9edf-548c-4958-82f4-e65d6fe99f3d
# ╠═af773761-1267-4dcb-9e60-a5358f6bb319
# ╠═e61162bf-c488-4c11-a446-11f65e3b8512
# ╠═f934ef47-7558-4960-98f4-c264772e23b2
# ╠═82e4a033-347b-423e-a8a3-6329f50023c0
# ╠═f4e11cf7-b68e-4881-b328-2640ef562240


text:
using StaticArrays

const BitBoard = UInt64
const BOARD_WIDTH = 8
const BOARD_HEIGHT = 8

@assert BOARD_WIDTH >= 1
@assert BOARD_HEIGHT >= 1

const BOARD_SIZE = BOARD_WIDTH * BOARD_HEIGHT
@assert BOARD_SIZE <= 8 * sizeof(BitBoard)

const IntPair = SVector{2,Int}
const ALL_COORDS = [SA[i, j] for i = 1 : BOARD_HEIGHT for j = 1 : BOARD_WIDTH]
const ROOK_DISPLACEMENTS = SA[SA[+1, 0], SA[-1, 0], SA[0, +1], SA[0, -1]]
const BISHOP_DISPLACEMENTS = SA[SA[+1, +1], SA[+1, -1], SA[-1, +1], SA[-1, -1]]

to_square(coords::IntPair) =
    @inbounds ((coords[1] - 1) * BOARD_WIDTH + (coords[2] - 1))

to_board(square::Int) = BitBoard(1) << square
to_board(coords::IntPair) = to_board(to_square(coords))

in_bounds(coords::IntPair) =
    @inbounds ((1 <= coords[1] <= BOARD_HEIGHT) &&
               (1 <= coords[2] <= BOARD_WIDTH))

function slider_mask(origin::IntPair, displs)
    result = BitBoard(0)
    for displ in displs
        temp = origin + displ
        while in_bounds(temp + displ)
            result |= to_board(temp)
            temp += displ
        end
    end
    return result
end

function slider_moves(origin::IntPair, displs, blockers::BitBoard)
    result = BitBoard(0)
    for displ in displs
        temp = origin + displ
        while in_bounds(temp)
            result |= to_board(temp)
            if (to_board(temp) & blockers) != 0
                break
            end
            temp += displ
        end
    end
    return result
end

function bitwise_subsets(x::BitBoard)
    if x == 0
        return BitBoard[0]
    else
        index = trailing_zeros(x)
        bit = BitBoard(1) << index
        y = x & (x - 1)
        rec = bitwise_subsets(y)
        result = BitBoard[]
        for z in rec
            push!(result, z)
            push!(result, z | bit)
        end
        return result
    end
end

move_table(origin::IntPair, displs) = [
    (blockers, slider_moves(origin, displs, blockers))
    for blockers in bitwise_subsets(slider_mask(origin, displs))
]

const BitBoardTable = Vector{Tuple{BitBoard,BitBoard}}
const ROOK_TABLES = [move_table(origin, ROOK_DISPLACEMENTS)
                     for origin in ALL_COORDS]
const BISHOP_TABLES = [move_table(origin, BISHOP_DISPLACEMENTS)
                       for origin in ALL_COORDS]

function is_valid_bit_count(table::BitBoardTable, bit_count)
    dict = Dict{BitBoard,BitBoard}()
    for (long_key, value) in table
        key = long_key >> bit_count
        if haskey(dict, key)
            if dict[key] != value
                return false
            end
        else
            dict[key] = value
        end
    end
    return true
end

function find_bit_count(table::BitBoardTable, magic_number::BitBoard)
    magic_table = [(magic_number * key, value) for (key, value) in table]
    lo = 0
    hi = 64
    for i = 1 : 6
        mid = lo + div(hi - lo, 2)
        if is_valid_bit_count(magic_table, mid)
            lo = mid
        else
            hi = mid
        end
    end
    return lo
end

view_board(board::BitBoard) = transpose(reshape(
    [((board >> i) & 1) != 0 for i = 0 : BOARD_WIDTH * BOARD_HEIGHT - 1],
    (BOARD_WIDTH, BOARD_HEIGHT)
))

function find_magic_numbers()
    magic_numbers = [(0, UInt64(0)) for _ in ALL_COORDS]
    while true
        magic_number = rand(BitBoard) & rand(BitBoard)
        for i = 1 : BOARD_SIZE
            table = ROOK_TABLES[i]
            bit_count = find_bit_count(table, magic_number)
            if magic_numbers[i][1] < bit_count
                magic_numbers[i] = (bit_count, magic_number)
                for j = 1 : BOARD_SIZE
                    if i == j
                        println('*', lpad(j, 2),
                                "*: ", lpad(magic_numbers[j][1], 2),
                                " : 0x", string(magic_numbers[j][2], base=16))
                    else
                        println(' ', lpad(j, 2),
                                " : ", lpad(magic_numbers[j][1], 2),
                                " : 0x", string(magic_numbers[j][2], base=16))
                    end
                end
                println()
            end
        end
    end
end

const ROOK_MAGIC_NUMBERS = UInt64[
    0x6580001040008120, 0x0cc0004410082001, 0x03800ad000832000, 0x31000804500100a0,
    0x02001830604e0044, 0x4600101402004b48, 0x02002804a2002409, 0x0200040240802506,
    0x060b8004804008a8, 0x2356400850016000, 0x00090040e0049101, 0x0202002150406a00,
    0x822a002a0a00b022, 0x8827001900040012, 0x6086000883042200, 0x4202002400c48601,
    0x134560800880c005, 0xb80c830040010428, 0x4143050040a00012, 0x0140990020b00100,
    0x0800808018001402, 0xc312008004008086, 0x34280c0010080d12, 0x648422000984004b,
    0x03c0028080002642, 0x4910084040012002, 0x160010c100200107, 0x8504633900100100,
    0x8071002700080090, 0x0a260106000c5028, 0x3a485d1c00a81006, 0x000188820005034c,
    0x880a401088800024, 0x0051104005002181, 0x4281801000802000, 0xd520801000800800,
    0x834200208e000870, 0x235200487200100c, 0x484630080c001691, 0x0012006492000403,
    0x4000882040018008, 0x7c20004530024000, 0x9189002000c30012, 0x24ab690090010020,
    0x10282cb801010010, 0x8002e40006008080, 0x1a225001085c0086, 0x0b12418420420001,
    0xcc44832200c90600, 0xcc44832200c90600, 0x004201b482604200, 0x0202002150406a00,
    0x0340280104008080, 0x8002e40006008080, 0x6086000883042200, 0x098d4ac08c110a00,
    0x00890512c1220082, 0x00890512c1220082, 0x00004d4011002001, 0xc0816900a00c1001,
    0x012100141800b003, 0x284200107c050812, 0x818761930210080c, 0x11800c010020da82
]

const ROOK_BIT_COUNT = 52

const BISHOP_MAGIC_NUMBERS = UInt64[
    0x5341280611c64100, 0xa660142d36042025, 0x45d00d5041018703, 0x0418061144002080,
    0xaa06021000113410, 0x2182c4206004a8c4, 0x1175413008220208, 0x1004124210042000,
    0x39820858b0740a41, 0xa08a2204030c0104, 0x4040410301190a09, 0x03000c0402888c30,
    0x8200360211016041, 0x100a02084c8418d8, 0x03000c0402888c30, 0x0043452103482028,
    0x4340300488080b30, 0x0022811812081600, 0x1110081807401260, 0x42920a4401220002,
    0xd124001205210640, 0xd032006501860140, 0x0043452103482028, 0x2401c0048a181102,
    0x617620434108ca10, 0x0031a03810040301, 0x00866203f0040442, 0x00400c0010410260,
    0x80160020020080c5, 0x12c802000d405200, 0x94010e0994025904, 0x0000c25121820801,
    0x8c0410159440ec10, 0x200c3a20818a0407, 0x28cc004110081208, 0x00b7601800110105,
    0x012c420400120090, 0x24020800404a004a, 0x818808088a804248, 0x2cc1104486120600,
    0xc42611106806c0a1, 0x500200a208082004, 0x0e01028041081001, 0x0026212051018800,
    0x0114a0061c008082, 0x0908090806010460, 0x82700442a0851c08, 0x8022008a22026081,
    0x1175413008220208, 0x04428a080b244459, 0x2004020462081680, 0x8494825184141064,
    0x4f0822d08a020020, 0x08344124980280a4, 0xa660142d36042025, 0x90bf900cc2c65aac,
    0x1004124210042000, 0x076c16494808036e, 0x0012160900824112, 0x02c11810ece08800,
    0x2c2a008820a0c100, 0x70042418b0311a04, 0x39820858b0740a41, 0x5341280611c64100
]

const BISHOP_BIT_COUNT = 55

function print_magic_tables()
    for i = 1 : BOARD_SIZE
        magic_table = zeros(UInt64, 1 << (64 - ROOK_BIT_COUNT))
        for (blockers, value) in ROOK_TABLES[i]
            key = ((blockers * ROOK_MAGIC_NUMBERS[i]) >> ROOK_BIT_COUNT) + 1
            @assert (magic_table[key] == 0) || (magic_table[key] == value)
            magic_table[key] = value
        end
        for i = 1 : 8 : length(magic_table)
            println("    ", string(magic_table[i:i+7])[8:end-1], ',')
        end
    end
end

print_magic_tables()